Analysis of the impact degree distribution in metabolic networks using branching process approximation
Kazuhiro Takemoto, Takeyuki Tamura, Yang Cong, Wai-Ki Ching, Jean-Philippe Vert, Tatsuya Akutsu
aa r X i v : . [ q - b i o . M N ] A ug Analysis of the impact degree distribution inmetabolic networks using branching processapproximation
Kazuhiro Takemoto a , b Takeyuki Tamura c Yang Cong d Wai-Ki Ching d Jean-Philippe Vert e , f , g Tatsuya Akutsu c , ∗ a PRESTO, Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012,Japan b Department of Biophysics and Biochemistry, University of Tokyo, Hongo 7-3-1,Bunkyo-ku, Tokyo 113-0033, Japan c Bioinformatics Center, Institute for Chemical Research, Kyoto University,Gokasho, Uji, Kyoto, 611-0011, Japan d Department of Mathematics, The University of Hong Kong, Pokfulam Road,Hong Kong e Centre for Computational Biology, Mines ParisTech, 35 rue Saint-Honor´e, 77305Fontainebleau cedex, France f Institut Curie, 75248 Paris, France g INSERM, U900, 75248 Paris, France
Abstract
Theoretical frameworks to estimate the tolerance of metabolic networks to variousfailures are important to evaluate the robustness of biological complex systems insystems biology. In this paper, we focus on a measure for robustness in metabolicnetworks, namely, the impact degree , and propose an approximation method to pre-dict the probability distribution of impact degrees from metabolic network struc-tures using the theory of branching process. We demonstrate the relevance of thismethod by testing it on real-world metabolic networks. Although the approxima-tion method possesses a few limitations, it may be a powerful tool for evaluatingmetabolic robustness.
Key words:
Metabolic network, Branching process, Power law, Cascading failure
PACS: ∗ Corresponding author.
Email address: [email protected] (Tatsuya Akutsu).
Preprint submitted to Elsevier 19 August 2018
Introduction
Robustness is a key feature in the analysis of complex systems, especiallyfor complex biological systems. Many organisms have strong adaptability toenvironmental changes or failures in some of their components, and can liveeven if some of their genes are mutated. In particular, it is known that cancercells are very robust [1]. Therefore, understanding the origin of robustness ofliving cells has become an important research topic.In particular, extensive studies have focused on the analysis of structural ro-bustness of metabolic networks. Structural robustness refers to the toleranceof the system’s behavior to changes in the structure of networks, and mostexisting studies focus on changes caused by knockout of gene(s) or enzyme(s).One of the reasons why extensive studies have been done on structural robust-ness of metabolic networks is that rather accurate data of metabolic networksare available via such databases as the Kyoto Encyclopedia of Genes andGenomes (KEGG) [2] and the Encyclopedia of
Escherichia coli
K-12 Genesand Metabolism (EcoCyc) [3], and kinetic parameters, which are not neces-sarily available, are not required.In order to analyze the structural robustness of metabolic networks, the fluxbalance analysis (FBA) methods have been widely used. In many of theseapproaches, elementary flux modes (EFMs) play a key role, where an EFMis a minimal set of reactions that can operate at steady state [4]. Based onFBA and/or EFM, several studies have focused on finding a minimum reac-tion cut [5,6,7,8], that is, a minimum set of reaction (or enzyme) removalswhich prevent the production of a specified set of compounds. Other FBA-based measures of robustness have also been proposed. Behre et al. proposeda measure based on the number of remaining EFMs after knockout versus thenumber of EFMs in the unperturbed situation [9]. Deutscher et al. proposedanother measure using the Shaply value from game theory [10].Other approaches have been proposed based on Boolean models of metabolicnetworks in which reactions and compounds are modeled as AND and ORnodes, respectively. Handorf et al. analyzed robustness of metabolic networksby introducing the concept of scope [11]. Li et al., Sridhar et al., and Tamuraet al. developed integer programming-based methods for finding a minimumreaction cut in Boolean models of metabolic networks [12,13,14,15]. Smart etal. defined the topological flux balance (TFB) criterion based on a Booleanmodel of metabolic networks and analyzed the damage (number of reactions)caused by knockout of a single reaction under TFB [16]. Jiang et al. definedand analyzed the impact degree , which is the number of reactions inactivatedby knockout of a specific reaction [17]. Although there are some differences inthe treatment of reversible reactions, the damage and the impact degree are2ery similar concepts. Cong et al. extended the impact degree for knockout ofmultiple reactions [18].In this paper, we study the distribution of the impact degree caused by randomknockout of a single reaction using the theory of branching process [19,20]. Inorder to analyze earthquakes, Saichev et al. proposed a branching processwith power-law distributions of offspring d : P ( d ) ∝ /d γ +1 , where γ is someconstant, and approximately derived the distribution of the total number ofoffsprings [20]. We regard propagation of the impact of knockout of a reac-tion as a branching process, and apply their method to estimate the impactdegree distribution, where the impact degree in our problem corresponds tothe total number of offsprings in the branching process. In order to applythis method, we develop a simple method for estimating the offspring distri-bution in a metabolic network. Although Smart et al. have already appliedpercolation theory and branching process to analysis of the size distributionof rigid clusters, defined as clusters of contagion nodes that do not containany branched metabolite nodes (see Fig. 3A in [16]), they did not explicitlyestimate the damage distribution (i.e., the impact degree distribution). Wefinally show an estimation method for the damage distributions. The pro-posed method is applied to analysis of metabolic networks of four species: Escherichia coli , Bacillus subtilis , Saccharomyces cerevisiae and
Homo sapi-ens . The results show good agreement of impact degree distributions betweenempirical results and theoretical estimates.
Jiang et al. proposed the impact degree as a measure of the importance ofeach reaction in a metabolic network [17]. The impact degree is defined asthe number of inactivated reactions caused by knockout of a single reaction.However, they did not consider the effect of cycles in metabolic networks. Sincecycles play an important role in metabolic networks, Cong et al. extended theimpact degree so that the effect of cycles is taken into account [18] by usinga concept of the maximal valid assignment [15]. Here, we briefly review thedefinition of this extended impact degree [18].As in other works, we regard each metabolic network as a bipartite directedgraph. Let V c = { C , . . . , C m } and V r = { R , . . . , R n } be a set of compoundnodes and a set of reaction nodes respectively, where V c ∩ V r = {} . A metabolicnetwork is defined as a directed graph G ( V c ∪ V r , E ) in which either ( u ∈ V c ) ∧ ( v ∈ V r ) or ( u ∈ V r ) ∧ ( v ∈ V c ) holds for each edge ( u, v ) ∈ E . Each reaction and compound takes one of two states: 0 or 1, where 0 and 1 A ∧ B means logical AND of A and B . R i is knocked out. Then, we start with the global statewhere all compounds are active (i.e., C k = 1 for all C k ∈ V c ) and all reactionsbut R i are active (i.e., R j = 1 for all R j ∈ V r \{ R i } and R i = 0). Then, wealternatively update the states of reactions and compounds by the followingrules.(1) For each reaction, there are three different compounds: consumed com-pounds (i.e., substrates), produced compounds (i.e., products), and di-rectly unrelated compounds.(2) A reaction is inactivated if any of its consumed or produced compoundis inactivated.(3) For each compound, there are three different reactions: consuming reac-tions, producing reactions, and directly unrelated reaction.(4) A compound is inactivated if all its consuming reactions or all its pro-ducing reactions are inactivated.Since no activation is possible in this process, the procedure converges to astable state in a finite number of iterations. The impact degree of reaction R i is defined as the number of inactivated reactions in the stable state. Thisprocedure simultaneously gives the definition of the impact degree and analgorithm to compute it. R R R R C R C C R C R C C C Fig. 1. Example of a metabolic network. Boxes and circles correspond to reactionsand compounds, respectively.
Let us illustrate the above process with the metabolic network shown in Fig. 1.Suppose that reaction R is knocked out. Then, the states of nodes change asshown in Table 1. Since four reactions (including R ) are inactivated in thestable state, the resulting impact degree is four. Next, suppose that reaction R is knocked out. In this case, the states of nodes change as shown in Table 2 andthe resulting impact degree is two. It is to be noted that C is not inactivatedbecause R is still active in this case. 4 able 1Impact degree calculation when R in Fig. 1 is knocked out. R R R R R R R C C C C C C C R in Fig. 1 is knocked out. R R R R R R R C C C C C C C We here explain the branching process approximation for estimating the im-pact degree distributions of metabolic networks.The branching process is a stochastic process in which each progenitor gener-ates offsprings according to a fixed probability distribution called the offspringdistribution. We propose that the branching process approximation is usefulfor estimating the impact degree distribution because the propagation of animpact on a network is essentially similar to cascading failures, which is asequence of failures caused by an accident. The branching process model ofcascading failure is a standard Galton-Watson branching process [21]. Theapproximation method using branching process has been already applied toloading-dependent cascading failure, and its relevance has been shown [22,23].However, the branching process approximation needs the assumption of treestructure of networks; thus, this is a limitation of the branching process ap-proximation because metabolic networks generally have cycles.5 .1 The number of offsprings in metabolic networks
To estimate the impact degree distributions using the branching process ap-proximation, we need to define the notion of offsprings for each reaction nodein metabolic networks, and the distribution of the number of offsprings.For that purpose, we consider the reaction network obtained as the unipartiteprojection of the metabolic network, where we draw an edge from reaction Ato reaction B when at least one product of A is a substrate of B [25,26,27,28].Fig. 2 shows the reaction network obtained from the metabolic network of Fig.1 by this procedure.As an easy example, we consider 2 metabolic reactions, A and B. In thiscase, the edge is drawn from A to B (i.e., A → B) if at least 1 product ofreaction A corresponds to at least 1 substrate of reaction B (e.g., the caseof a → A → b → B → c, where a, b, and c are chemical compounds). Througha similar procedure, we obtain a reaction network from a given metabolicnetwork (see Fig. 2). R R R R R R R Fig. 2. The reaction network transformed from the metabolic network in Fig. 1.
We now observe that the number of reactions inactivated by the failure of agiven reaction, which we want to model in the branching process, does not cor-respond to its outdegree (i.e., the number of out-going edges from a reactionnode) in the reaction network because the spreading of an impact is depressedwhen there are alternative synthetic pathways. For example, assuming thatthe reaction R is inactive (or disrupted), the cascading of the impact does notoccur because the reaction R remains active due to the chemical compoundgenerated through the reaction R . In contrast, the cascading of the impactcontinues when the reaction R becomes inactive for instance because the re-action R is dependent on this reaction only. In short, the impact spreadsthough reactions whose substrates are synthesized via unique metabolic reac-tions (i.e., reaction nodes with the indegree of 1) when assuming tree structuresof networks.Based on this analysis, we define the number of offsprings for each reaction6ode in metabolic networks as follows: d i = k out i (if k in i = 1)0 (otherwise) , (1)where k out i and k in i are the outdegree and indegree of reaction node i in thereaction network, respectively. To analytically estimate the impact degree distributions, it is useful to assumethat the number of offsprings for each node follows a probability distribution.We here consider two types of distributions, which are frequently observed inreal worlds.
The simplest case of branching processes is a Poisson branching process (here-after called
Poisson model ) in which the number of offsprings d for each pro-genitor follows the Poisson distribution: µ d e − µ /d !, where µ corresponds to themean of this distribution. In this model, the total number of offsprings (i.e.,impact degree) r is distributed according to the Borel distribution [24]: P ( r ) = ( µr ) r − e − µr r ! . (2)Using Stirling’s formula (i.e., r ! ≈ √ πrr r e − r ), the above equation leads tothe approximation P ( r ) ∝ r − / e − r (ln µ − µ +1) . (3)In particular, when µ = 1 (i.e., the critical case), the impact degree follows apower-law distribution with exponent − / In addition to the Poisson case, we consider the case where the number ofoffsprings is determined based on a power-law distribution (hereafter called
Power-law model ). Indeed, the number of offsprings for each reaction nodeis based on the outdegree in metabolic networks. Since real-world complex7etworks including metabolic networks have power-law degree distributions[25,29,30], the number of offsprings for each reaction node may also obey apower-law distribution.Saichev et al. [20,33] showed analytical asymptotic approximations for thedistribution P ( r ) of the total number of offsprings (i.e., impact degree) r inthe case where the number of offsprings d for each progenitor follows asymp-totically a power-law distribution P ( d ) ∝ /d γ +1 and the mean number ofoffsprings has a given value µ . In particular they derive the following approx-imation for large impact degree r , in the case 1 < γ < P ( r ) ≃ µνr /γ ϕ γ (1 − µ ) r − µ − νr /γ ! , (4)where ϕ γ ( x ) = ∞ Z exp (cid:20) u γ cos (cid:18) πγ (cid:19)(cid:21) cos (cid:20) u γ sin (cid:18) πγ (cid:19) + ux (cid:21) du , (5)and ν = µ ( γ − γ /γ − Γ( − γ ) /γ , where Γ( x ) is the Gamma function.When µ = 1 (i.e., the critical case), in particular, the distribution of the totalnumber of offsprings obeys the power-law distribution: P ( r ) ∝ /r /γ .In addition, P ( r ) in the case of γ > The above probability distributions may be unsuitable to approximate real-world offspring distributions. In addition, we also consider a branching pro-cess model using empirical offspring distributions (hereafter called
Empiricalmodel ). This way, we can estimate the distribution P ( r ) of the total numberof offsprings (i.e., impact degree distributions) without the approximation ofoffspring distributions, although the model is analytically intractable. More-over, we can evaluate whether the prediction accuracy of the proposed methodis influenced by the approximation of offspring distributions or the fidelity ofbranching processes.Let F ( s ) be the probability generating function of the impact degree r (i.e., thetotal number of offsprings), the function F ( s ) satisfies the recursive relation820,21]: F ( s ) = f ( sF ( s )) , (6)where f ( s ) denotes the probability generating function of the number d ofoffsprings of each node.Using the Lagrange expansion and the relation of P ( r ) = (1 /r !) d r F ( s ) /ds r | s =0 ,the distribution P ( r ) (i.e., impact degree distribution) is derived from theabove implicit equation as the following explicit equation [20,21]: P ( r ) = 1 r ! d r − ds r − " f r ( s ) df ( s ) ds s =0 ( r > , (7)where f ( s ) = P d max d =0 P ( d ) s d . The value d max indicates the maximum of d , andthe function P ( d ) corresponds to the probability density function of empirical d (i.e., empirical offspring distribution). In addition, P ( r ) = f (0) when r = 0. To apply the Poisson and the power-law models, we need to estimate themodel parameters, namely the mean µ and the exponent γ , from real metabolicnetworks.We estimate the mean of the number of offsprings for each reaction node bythe empirical average: µ = 1 N N X i =1 d i , (8)where N is the total number of reaction nodes in a metabolic network.We estimate the exponent of a power-law offspring distribution using the max-imum likelihood estimation method [34]: γ = | N ∗ | " X i ∈ N ∗ ln d i d min − , (9)where N ∗ is the set of reaction nodes with d i >
0, and | N ∗ | indicates the totalnumber of such reaction nodes. d min is the minimum of d i in the set of N ∗ .9 Evaluation of the branching process approximation
We evaluated the above estimation methods for the impact degree distribu-tions on several real metabolic networks.We selected two bacteria [
Escherichia coli (eco) and
Bacillus subtilis (bsu)]and two eukaryotes [
Saccharomyces cerevisiae (sce) and
Homo sapiens (hsa)]whose metabolic pathways have been well-identified. We downloaded the dataof their metabolic networks, represented as bipartite networks as shown in Fig.1, from the KEGG database (version 0.7.1) [2,31]. The parenthetic three-lettercodes correspond to KEGG organism identifiers [32].Based on the KEGG metabolic network data, the impact degree distribu-tions in the metabolic networks were calculated using the method explainedin Sec. 2. Moreover, we constructed the reaction networks of these species,and obtained the offspring distributions. Using Eqs. (8) and (9), the modelparameters µ and γ were extracted (see Table 3). All metabolic networks show1 < γ <
2, implying that the assumption of power-law model is suitable.
Table 3Model parameters extracted from real metabolic networks. The character µ Exponent γ Escherichia coli
Bacillus subtilis
Saccharomyces cerevisiae
Homo sapiens
To test the validity of the offspring distribution models, we compared thefitting results on empirical offspring distributions in real metabolic networksbetween the power-law distribution and the Poisson distribution. Fig. 3 showsthe cumulative offspring distributions from the real metabolic networks andthe cumulative representation of theoretical distributions. The figure clearlyindicates that the power-law distributions are more appropriate for modelingoffspring distributions than the Poisson distributions. However, the power-lawdistributions may be not the best model because of the poor fittings for thelarger d in the case of H. sapiens (Fig. 3D)Using Eqs. (2), (4), and (7), we obtained the estimated impact degree distribu-tions using the branching process approximation. Fig. 4 shows the comparisonbetween the observed cumulative impact degree distributions and estimatedones. The theoretical predictions (lines) are in good agreement with the real10 ower-lawPoisson 10 -3 -2 -1 P c u m ( d ) -3 -2 -1 P c u m ( d ) d10 d 10 -4 -3 -2 -1 P c u m ( d ) d10 -4 -3 -2 -1 P c u m ( d ) d (B)(A)(C) (D) Fig. 3. Cumulative offspring distributions of
Escherichia coli (A),
Bacillus subtilis (B),
Saccharomyces cerevisiae (C), and
Homo sapiens (D). P cum ( x ) is defined as P ( X ≥ x ). Note that P cum (1) < d = 0. P cum (0) is not shown due to the logarithmic display. Thesymbols indicate observed data. The black solid lines and dashed lines correspondto the cumulative representations of the power-law distribution with the exponentestimated by Eq. (9) and the Poisson distribution with the mean obtained from Eq.(8), respectively. impact degree distributions (symbols), suggesting the relevance of branchingprocess approximations. Note that the impact degree distributions does notfollow a clear power law and show an exponential cut-off for larger impactdegrees because µ < ower-lawPoissonEmpirical 10 -3 -2 -1 P c u m (r) -3 -2 -1 P c u m (r) r10 r 10 -3 -2 -1 P c u m (r) r10 -3 -2 -1 P c u m (r) r (B)(A)(C) (D) Fig. 4. Cumulative impact degree r distributions of Escherichia coli (A),
Bacillussubtilis (B),
Saccharomyces cerevisiae (C), and
Homo sapiens (D). The symbolsindicate observed data. The black solid lines and dashed lines correspond to thecumulative representations of theoretical distributions of the power-law model andthe Poisson model, respectively. The gray solid lines are the the cumulative repre-sentations of theoretical distributions from the empirical model. Note that P cum ( x )is defined as P ( X ≥ x ). tions and theoretical distributions (Table 4). The power-law model is betterthan the Poisson model on all networks. The empirical model outperformsthe power law model on both bacterial networks. Surprisingly, the power lawmodel outperforms the empirical model on both eukaryote’s networks.12 able 4Prediction accuracy for the impact degree distributions: Kolmogorov-Smirnov (KS)distance, defined as sup x | R ( x ) − M ( x ) | , where R ( x ) and M ( x ) are empirical distri-butions and theoretical distributions, respectively. The parenthetic values indicatethe logarithmic P -values p from the KS test, defined as − log ( p ). The emphasizedvalues correspond to the best accuracy.Species Poisson model Power-law model Empirical model Escherichia coli
Bacillus subtilis
Saccharomyces cerevisiae
Homo sapiens
We proposed a model to estimate the impact degree distributions in metabolicnetworks, using a branching process approximation, and demonstrated its va-lidity on real data.The power-law model could more accurately estimate the impact degree distri-butions in real metabolic networks than the Poisson model because the numberof offsprings for each reaction node is assumed to follow the power-law dis-tribution. Especially, the power-law model showed the significant agreementsbetween the predicted distribution and the observed distribution although thecase of
H. sapiens represented the small P -value for the KS test (i.e., the lowprobability that the distribution is similar between models and observed data).However, there is no great difference of the prediction accuracy for estimatingthe impact degree distributions between the power-law model and the Poissonmodel; thus, the Poisson model may be useful for a rough estimate of theimpact degree distributions.Intrinsically, the distribution of the total number of offsprings [i.e., P ( r )] is notsignificantly different between the power-law model and the Poisson model inthe case of smaller r . As a simple example, we here consider the case of µ = 1(i.e., the critical case). In this case, the impact degree distributions P ( r ) ofthe power-law model and the Poisson model correspond to ∝ r − (1+1 /γ ) and ∝ r − / , respectively. Especially, P ( r ) of the power-law model is the power-law distribution with the exponent ranging between − − . < γ <
2; thus it is not critically different from r − / of the Possion model forsmaller r . Since the impact degrees of real metabolic networks (i.e., r ) wererelatively small ( r < P ( r ) of the total number of offsprings has the universal property ofa power-law tail with the exponent − / r → ∞ under mild conditions onoffspring distributions, and it implies that the types of offspring distributionshardly influence the distribution of the total number of offsprings (i.e., impactdegree distributions). Note the Otter’s theorem does not contradict with theanalytical distribution P ( r ) of the power-law model because this theorem isnot directly applicable to the power-law model due to the different assump-tions in the derivation of P ( r ) between the power-law model and the Otter’stheorem.The prediction performance is influenced by the assumption of offspring distri-butions and the fidelity of branching processes. To purely evaluate the validityof branching process approximation, the empirical model is useful because ofno assumption of offspring distributions. It is expected that the empiricalmodel show the best prediction accuracy because of using empirical offspringdistributions. In the case of bacteria (i.e., E. coli and
B. subtilis ), this expec-tation is true, suggesting that the branching process approximation is usefulfor estimating the impact degree distributions. In the case of eukaryotes (i.e.,
S. cerevisiae and
H. sapiens ), on the other hand, we observed the unexpectedresults: the empirical model shows the relatively-low prediction accuracy. Es-pecially, the prediction accuracy of the empirical model is lowest in the caseof
H. sapiens (human). This result implies limitations to the estimation ofimpact degree distributions based on the branching process approximation inthe case of metabolic networks of eukaryotes (i.e., higher organisms).A limitation of the branching process approximation is that we need to assumetree structures of networks (i.e., no cycles). The presence of cycles may leadto an overestimation of the number of offsprings d i for each reaction, becausesome offsprings of a progenitor (reaction node) may have already been inacti-vated due to cycle structures. From this reason, the models may overestimatethe impact degree distributions. On the other hand, however, some reactionswith more than one incoming edge may be inactivated in the presence of cy-cles, if all their parents are inactivated through different paths. In this case,the number of offsprings is underestimated. The estimation of the numberof offsprings depends on the relative importance of these two effects, and itmay be not simple. The difficulty in this estimation is also a limitation of themodel.The branching process model needs to be improved by considering additionalassumptions other than metabolic network structures in order to obtain bet-ter predictions. For example, the assumption of variable propagations in thebranching process [35], in which the mean of offspring distributions differs14t each propagation stage, may be useful because the degree of propagationmay depend on metabolic dynamics such as gene expressions and metaboliteconcentrations. To apply this modified branching process model, time-seriesdata on metabolic dynamics after gene disruptions, which are obtained bymetabolomic analysis, are necessary for estimating the mean of offspring distri-bution at each propagation stage (i.e., time after gene disruptions). Since suchdata are unavailable at present, however, it is difficult to apply this modifiedbranching process model. Similarly, other biologically-suitable assumptions arehardly determined because of few observed data on metabolic dynamics. Al-though there are above constraints on observed data on metabolic dynamics,we believe that the consideration of the additional information improves theprediction of impact degree distributions. In the future, the improvement ofthe prediction of the impact degree distributions using the branching processapproximation may be possible with available data on metabolic dynamics.The branching process approximation is useful for estimating the impact de-gree distributions in metabolic networks although it has the above limitations;thus, it may be a powerful tool for evaluating a robustness of biological sys-tems. Acknowledgment
This work was partially supported by the JST PRESTO program. TA, TT andJPV were partially supported by a JSPS/INSERM grant. TA was partiallysupported by Grant-in-Aid
References [1] H. Kitano, Cancer as a robust system: Implications for anticancer therapy, Nat.Rev. Cancer 4 (2004) 227–235.[2] M. Kanehisa, S. Goto, M. Furumichi, M. Tanabe, M. Hirakawa, KEGG forrepresentation and analysis of molecular networks involving diseases and drugs,Nucleic Acids Res. 38 (2010) D355–D360.[3] I. M. Keseler et al., EcoCyc: a comprehensive database of
Escherichia coli biology,Nucleic Acids Res., 39 (2011) D583–D590.[4] J. A. Papin, J. Stelling, N. D. Price, S. Klamt, S. Schuster, B. O. Palsson,Comparison of network-based pathway analysis methods, Trends Biotechnol. 22(2004) 400–405.
5] V. Acu˜na, F. Chierichetti, V. Lacroix, A. Marchetti-Spaccamela, M. F. Sagot,L. Stougie, Modes and cuts in metabolic networks: Complexity and algorithms,Biosystems 95 (2009) 51–60.[6] A. P. Burgard, P. Pharkya, C. D. Maranas, OptKnock: a bilevel programmingframework for identifying gene knockout strategies for microbial strainoptimization, Biotechnol. Bioeng. 84 (2003) 647–657.[7] U. U. Haus, S. Klamt, T. Stephen, Computing knock-out strategies in metabolicnetworks, J. Comput. Biol. 15 (2008) 259–268.[8] S. Klamt, E. D. Gilles, Minimal cut sets in biochemical reaction networks,Bioinformatics 20 (2004) 226–234.[9] J. Behre, T. Wilhelm, A. von Kamp, E. Ruppin, S. Schuster, Structuralrobustness of metabolic networks with respect to multiple knockouts, J. Theor.Biol. 252 (2008) 433–441.[10] D. Deutscher, I. Meilijson, S. Schuster, E. Ruppin, Can single knockoutsaccurately single out gene functions? BMC Syst. Biol. 2 (2008) 50.[11] T. Handorf, O. Ebenh¨oh, R. Heinrich, Expanding metabolic networks: scopesof compounds, robustness, and evolution, J. Mol. Evol. 61 (2005) 498–512.[12] N. Lemke, F. Her´edia, C. K. Barcellos, A. N. dos Reis, J. C. M. Mombach,Essentiality and damage in metabolic networks, Bioinformatics 20 (2004) 115–119.[13] Z. Li, R.-S. Wang, X.-S, Zhang, L. Chen, Detecting drug targets with minimumside effects in metabolic networks, IET Syst. Biol. 3 (2009) 523–533.[14] P. Sridhar, B. Song, T. Kahveci, S. Ranka, Mining metabolic networks foroptimal drug targets, In: Proc. Pacific Symposium on Biocomputing 2008 (2008)291–302.[15] T. Tamura, K. Takemoto, T. Akutsu, Finding minimum reaction cutsof metabolic networks under a Boolean model using integer programmingand feedback vertex sets, International Journal of Knowledge Discovery inBioinformatics 1 (2009) 14–31.[16] A. G. Smart, L. A. N. Amaral and J. M. Ottino, Cascading failure androbustness in metabolic networks, Proc. Nat. Acad. Sci. USA 105 (2008) 13223–13228.[17] D. Jiang, S. Zhou, Y-P. P. Chen, Compensatory ability to null mutation inmetabolic networks, Biotechnol. Bioeng. 103 (2009) 361–369.[18] Y. Cong, T. Tamura, T. Akutsu, W-K. Ching, Efficient computation of impactdegrees for multiple reactions in metabolic networks with cycles, In: Proc. ACMThird International Workshop on Data and Text Mining in Bioinformatics (2009)67.