Partial inhibition and bilevel optimization in flux balance analysis
 Giuseppe Facchetti^{1} and
 Claudio Altafini^{2}Email author
DOI: 10.1186/1471210514344
© Facchetti and Altafini; licensee BioMed Central Ltd. 2013
Received: 13 February 2013
Accepted: 19 November 2013
Published: 29 November 2013
Abstract
Motivation
Within Flux Balance Analysis, the investigation of complex subtasks, such as finding the optimal perturbation of the network or finding an optimal combination of drugs, often requires to set up a bilevel optimization problem. In order to keep the linearity and convexity of these nested optimization problems, an ON/OFF description of the effect of the perturbation (i.e. Boolean variable) is normally used. This restriction may not be realistic when one wants, for instance, to describe the partial inhibition of a reaction induced by a drug.
Results
In this paper we present a formulation of the bilevel optimization which overcomes the oversimplified ON/OFF modeling while preserving the linear nature of the problem. A case study is considered: the search of the best multidrug treatment which modulates an objective reaction and has the minimal perturbation on the whole network. The drug inhibition is described and modulated through a convex combination of a fixed number of Boolean variables. The results obtained from the application of the algorithm to the core metabolism of E.coli highlight the possibility of finding a broader spectrum of drug combinations compared to a simple ON/OFF modeling.
Conclusions
The method we have presented is capable of treating partial inhibition inside a bilevel optimization, without loosing the linearity property, and with reasonable computational performances also on large metabolic networks. The more finegraded representation of the perturbation allows to enlarge the repertoire of synergistic combination of drugs for tasks such as selective perturbation of cellular metabolism. This may encourage the use of the approach also for other cases in which a more realistic modeling is required.
Background
In recent years, genomescale metabolic networks have represented an important paradigm of systems biology, well describing how interesting (and relevant) biological features can be deduced in spite of the complexity of the model [14]. Thanks to the use of genomic techniques, metabolic networks have been reconstructed for many organisms, ranging from small bacteria up to the human cells. In parallel, the development of quantitative descriptions of these large and complex systems based on simple computational framework such as Flux Balance Analysis (FBA) [5, 6] has increased both their characterization [79] and the spectrum of applications. Two important examples are (i) strain improvement [10, 11], i.e. the identification of the best knockout or gene manipulation maximizing the biosynthesis of a key metabolite, (ii) support to drug discovery through the identification of new inhibition targets [1215] or of new drug therapies for various medical purposes [16, 17]. All the studies just mentioned are based on the FBA formalism.
The constraints (1) and (2) generate a convex and bounded set $W\subset {\mathbb{R}}^{r}$ to which the vector v has to belong. To obtain the reaction fluxes (a point in W) which describe the metabolic state of the organism, one has to perform the maximization of a function Φ(v) (or equivalently the minimization of Φ(v)). The choice of Φ(v) depends on the context and on the application: common examples are biomass production [18], ATP production [2] and minimal metabolic adjustment [19].
where arg min_{x∈Q}f(x) stands for the set of x for which the function f attains its minimum value in Q (or equivalently, its maximum value for "arg max"). Therefore, (3) says that the output of the bilevel optimization is the vector h such that the corresponding vector v(h), which minimizes Φ on W(h), maximizes the function Ψ. The reformulation of a bilevel optimization as a single optimization is commonly obtained through duality theory [21, 22]. When performing this reformulation, in order to save the linear nature of the optimization procedure, the variables h (the real output of the algorithm) have to be Boolean, rather than continuous [2, 10, 13]. While this approach is correct for gene knockout, in the case of drug treatment (where the enzymes are inhibited by drug) or gene modification (on which one changes the enzymes activity) it represents only a rough approximation which may not constitute a realistic description of the biological effect. It is in fact more plausible to assume that a drug acting on an enzyme leads to a partial loss of functionality of the latter, and hence to a partial inhibition of the corresponding reaction(s). In order to treat the inhibition as a variable to be optimized and to avoid the ON/OFF oversimplification it is necessary to reformulate the bilevel optimization as a nonlinear (nonconvex) single optimization problem [23] but this leads to a more complicated situation from a numerical point of view.
Although a complete inhibition of a diseasecausing target may not represent the right therapeutic solution (in healthy cells, the level of each metabolite must be in a finite range) and although expected to be a potential strategy in a multitarget approach [24], partial inhibition has been considered only in a few computational works; moreover, studies like [25, 26] dealt with a small part of the network, modeling the kinetic reactions explicitly and solving then numerically. The partial inhibition then amounts, for instance, to a modulation of one or more kinetic parameters. Because of the complexity of the metabolic networks and because of the impossibility of knowing the kinetic parameters of all biochemical reactions, this approach cannot be applied at a genomewide level. On the other hand, the authors of [27] consider the whole network, but the (partial) inhibition is given as initial fixed parameter of the model and only the effect of the perturbation is quantified. A different approach is presented in [28] in the context of the prediction of new drug targets; these targets are identified though a twostages FBA (which differs from a bilevel formulation because the two optimizations are not nested). However, the potential targets obtained with this method must be verified exhaustively, which may represent a problem for networks with more than the 26 reactions of the human hyperuricemia metabolic pathway considered in [28].
Therefore, the aim of the present paper is to describe a novel algorithm which allows to provide a more realistic description of the partial inhibition induced by the drugs on large networks while still remaining within the framework of Linear Programming (LP). In order to introduce the algorithm, we refer to a realistic case where a bilevel minimization is used. Namely we consider the search for the optimal combination of drugs capable, through a synergistic effect, to inhibit (or enhance) an objective reaction (i.e. a putative target for a disease) while inducing the minimal perturbation on the rest of the network. Indeed, the selectivity of the therapy is one of the most important aspects of any drug discovery project. Replacing a single Boolean variable by a convex combination of a fixed number of Boolean variables, we are able to model the inhibition as any number belonging to a discretized representation of the interval [0,1]. This approach preserves the linear nature of the final problem. Notice that the method we propose can be extended to any bilevel optimization which needs to deviate from the simple ON/OFF description.
The paper is organized as follows. We first formalize our example about drug combinations; then, within this case study, we describe why Boolean variables are necessary in the reformulation of a bilevel optimization problem via the strong duality theorem of LP. The presentation of the basic idea of the proposed algorithm and the discussion of its limits conclude the Methods Section. Results and Discussion sections present and comment the outcomes obtained on a benchmark application to the E.coli core metabolism and to some other larger networks. Final considerations are then reported in the Conclusion.
Methods
Optimal drug combination: a guiding example
Symbol, value range, meaning and type of all quantities used in the algorithm description
Network and drugs  
Symbol  Range  Description  Type 
r 
 Number of reactions  Parameter 
m 
 Number of metabolites  Parameter 
d 
 Number of drugs  Parameter 
S  ${\mathbb{R}}^{m\times r}$  Stoichiometric matrix of the network  Parameter 
${v}_{1}^{\text{ut}},\dots ,{v}_{r}^{\text{ut}}$  ${\mathbb{R}}_{+}$  Reaction fluxes of the untreated network  Parameter 
${v}_{1}^{\text{tr}},\dots ,{v}_{r}^{\text{tr}}$  ${\mathbb{R}}_{+}$  Reaction fluxes after drug treatment h  Variable 
mod  1,…,r  Index of the reaction which must be modulated  Parameter 
W    Set of feasible solutions for unperturbed network  Parameter 
W(h)    Set of feasible solutions for network perturbed by drug inhibition h  Parameter(*) 
   Set of drugs ($\left\mathcal{D}\right=d$)  Parameter 
${\mathcal{T}}_{k}$    Set of targets of drug k (k = 1,…,d)  Parameter 
Outer problem  
Symbol  Range  Description  Type 
τ  (0,1)  Threshold for the reaction that must be inhibited (${v}_{\text{mod}}^{\text{tr}}\le \tau {v}_{\text{mod}}^{\text{ut}}$)  Parameter 
>1  Threshold for the reaction that must be activated (${v}_{\text{mod}}^{\text{tr}}\ge \tau {v}_{\text{mod}}^{\text{ut}}$)  
P  ${\mathbb{N}}_{0}$  Precision parameter  Parameter 
h_{1},…,h_{ d }  [0,1]  Drug inhibitions  Variable 
x_{1,0},x_{1,1},…,x_{d,P}  {0,1}  Drugdosage variables  Variable 
b  10^{3}  Correction against an "overselection" of drugs  Parameter 
Inner problem  
Symbol  Range  Description  Type 
v_{1},…,v_{ r }  ${\mathbb{R}}_{+}$  Reaction fluxes (primal problem and outer problem)  Variable 
L_{1},…,L_{ r }  ${\mathbb{R}}_{+}$  Lowerbounds for the reaction fluxes (of the primal problem)  Parameter 
U_{1},…,U_{ r }  ${\mathbb{R}}_{+}$  Upperbounds for the reaction fluxes (of the primal problem)  Parameter 
a_{1},…,a_{ r }  ${\mathbb{R}}_{+}$  Absolute fluxes differences ${a}_{i}={v}_{i}{v}_{i}^{\text{ut}}$ (primal problem)  Variable 
μ_{1},…,μ_{ m } 
 Dual variables associated to the FBA steady state constraint  Variable 
λ_{1},…,λ_{ r }  ${\mathbb{R}}_{+}$  Dual variables associated to the unperturbed upperbounds  Variable 
δ_{1},…,δ_{ t }  ${\mathbb{R}}_{+}$  Dual variables associated to the drug inhibition; $t:={\sum}_{k=1}^{d}{\mathcal{T}}_{k}$  Variable 
${\delta}_{1}^{\text{max}},\dots ,{\delta}_{t}^{\text{max}}$  ${\mathbb{R}}_{+}$  Upperbounds of the dual variables δ’s  Variable 
α_{1},…,α_{ r }  ${\mathbb{R}}_{+}$  Dual variables associated to the first absolute value inequalities  Variable 
β_{1},…,β_{ r }  ${\mathbb{R}}_{+}$  Dual variables associated to the second absolute value inequalities  Variable 
In the following the side effect of a drug treatment is quantified in terms of the distance ∥v^{tr}(h)  v^{ut}∥_{1} used in (5): the greater the distance, the bigger the impact of the drugs on the whole network.
The problem can be stated as follows:
Given:

a metabolic network, which means a stoichiometric matrix$\mathbf{\text{S}}\in {\mathbb{R}}^{m\times r}$ and the upperbounds$\mathbf{U}\in {\mathbb{R}}^{r}$of the reaction fluxes v;

the unperturbed fluxes v^{ut};

the set of drugs together with their inhibition targets${\left\{{\mathcal{T}}_{k}\right\}}_{k=1,\dots ,\left\mathcal{D}\right}$,

the index (denoted by "mod") of the objective reaction whose flux (v_{mod}) must be modulated;

a threshold τ ∈ [0,1) for the modulation constraint on v_{mod};
we want to find the inhibition h ∈ [0,1]^{ d } such that${v}_{\text{mod}}^{\text{tr}}(\mathbf{h})\le \tau {v}_{\text{mod}}^{\text{ut}}$and such that it causes the minimal side effect (i.e. the minimal distance ∥v^{tr}(h)v^{ut}∥_{1}). Of course, a different definition of side effect as well as a different constraint on ${v}_{\text{mod}}^{\text{tr}}$ (perhaps on its maximal value [17]) can be used if needed by the problem.
The bilevel optimization (6) is a minmin linear program. The inner problem adjusts the fluxes so as to achieve the minimal metabolic adjustment, subject to the drug inhibitions imposed by the outer problem and to the stoichiometric constraints. The outer problem selects the combination of drugs which has the minimum side effect and guarantees a modulated flux lower than the desired threshold.
The strong duality theorem and the need of Boolean variables
providing that both sets are not empty"[22], where x are the primal variables and θ are the dual variables (note that, because of the use of the transpose of the matrix A in the dual problem, there is a dual variable for each constraint of the primal problem). Therefore, the application of this theorem to the inner problem consists in appending a list of constraints, A^{ T } θ ≥ c, corresponding to the dual form of the constraints of the inner problem and setting the inner objective function equal to its dual c^{ T }x = b^{ T }θ (see Additional file 2 for a depiction of the structure of the matrix eventually obtained). Since, from the strong duality theorem, this equality holds only at the optimal points of both primal and dual problems, the resulting set of constraints is equivalent to selecting only the solutions of the inner problem.
where ${\delta}_{i}^{\mathit{\text{max}}}$ is the upper bound for the dual variable δ_{ i }.
The restriction to Boolean variables saves the linear nature of the problem (which however requires now Mixed Integer Linear Programming) but it implies the assumption that drugs can only act as switches on the reactions, or equivalently, that we are considering only an ON/OFF model.
Partial inhibition
Inhibitions and activation of the objective reaction
The evaluation of the effect on the fluxes induced by the drugs is performed through the MOMA formalism. It is known that this approach describes well the spreading across the network of the effect of the perturbation: many processes are downregulated or upregulated in order to adjust and compensate the effect of the perturbation (see for example [4]). For the same reason, a second intervention may amplify the deactivation (recovery) of a certain metabolic function that was downregulated (activated) after the first perturbation [30]. In terms of multiple drug effect, this means that a drug synergism may reinforce both the inhibition and the activation of the reaction fluxes.
(for τ > 1) the algorithm generates drug interactions that upregulates the objective reaction.
In the following, both versions are applied.
Cases of multiple equivalent solutions (nonuniqueness)
where the parameter b is chosen small enough (10^{3} in our computations, much smaller than 10^{3}, the common upperbounds of the fluxes) in order to keep this term smaller than the difference in the side effects and therefore not to change the order between nonequivalent solutions.
Since the problem is not strictly convex and since equivalent drug combinations may always appear (for instance, the combination of drug A which targets reaction 1 plus drug B which targets reactions 2 and 3 is equivalent to the combination of drug C which targets reaction 2 plus drug D which targets reactions 1 and 3), the problem of multiple equivalent solutions needs to be considered. In all these cases, because of the numerical implementation of the algorithm, a random choice of one of these optimal solutions is taken. Since, by definition, all these equivalent solutions fulfill the conditions on the objective reaction and on the minimization of the side effect, their differences are irrelevant: indeed they concern only some other fluxes on which we do not have any specific requirements. For this reason any solution chosen by the implementation of the algorithm can be considered acceptable.
Results
Computational performaces
The proposed algorithm has been implemented in MATLAB (2012R) and the optimization has been performed using ILOGIBM CPLEX 12.1 ( http://www01.ibm.com/software/commerce/optimization/cplexoptimizer/) under academic license.
Main properties of the metabolic networks used in this paper (from BIGG http://bigg.ucsd.edu/ )
Moreover, we run the algorithm on the metabolic network of the six microorganisms listed in Table 2. Our scope is to evaluate the impact of the size of the network (parametrized by the number of reactions r) on the computational performaces. In order to limit the interference of other parameters, these calculations are carried out with the same objective reaction (in particular we still keep ribose5phosphate isomerase since it appears on all networks we have considered) at constant precision (P = 2) and threshold (τ = 0.35), with the same number of drugs (d = 8), and choosing their inhibition targets in a random manner (unfeasible problems are ignored). However, since on very large networks it is quite unlikely to induce the sought modulation on the objective reaction when only 8 targets are inhibited, the number of targets of each drug is proportionally increased (on average the total number of inhibitions is approximately 6% of the total number of reactions). Because of the randomness in the choice of the targets, the computational times may present a significant variation. Therefore, Figure 3 (b) shows the whole distributions of the computational time over 100 runs for each one of the six metabolic networks we considered. Finally, Figure 3 (c) reports the mean and the standard deviation of these distributions as a function of the size of the network. On average, also for very large networks, the computational time is approximately one hour (on a 2.3 GHz CPU). All these characterizations show the good performaces of the algorithm.
Screening for optimal drug combinations
For each pair, we perform a screening that considers each metabolic reaction as objective process to be modulated (down or upregulated depending on the value of τ) and finds the most selective drug combination. The following characterization of the solutions is performed. For a given P and for a given objective reaction v_{mod}, we consider the solutions h at different values of τ. When the same drug combination is found for two values of τ (for example τ_{1} = 0.1 and τ_{2} = 0.5), the solution is considered valid only for the most stringent constraint (τ_{1} = 0.1, in the example; similarly, if τ_{1} = 1.5 and τ_{2} = 2.0 then the solution is associated to τ_{2} = 2.0 only). This procedure allows to considered only cases when passing to a weaker constraint on v_{ mod } the severity (for instance the dosage) of the corresponding optimal drug treatment is reduced too. We analyze the results by looking at the following four indices.
Averaged cardinality of the drug combinations from the screening
τ = 0.0  τ = 0.1  τ = 0.5  τ = 1.5  τ = 2.0  

P = 0  2.5  2.9  2.1  2.1  1.4 
P = 1  2.7  3.0  2.3  3.5  3.1 
P = 2  2.9  3.3  2.4  3.7  3.2 
Cardinality of the solutions: More details are presented in panel (b) of Figure 5, which reports the histogram of the cardinality of the solutions and their mean values for the case of τ = 0.5 (averages for each value of τ are reported in Table 3). When the precision increases, the distribution of the cardinality shifts slightly to higher values, meaning that multiple drug treatments are (slightly) preferred.
Perturbation induced by the solutions: For each solution that we have identified during this screening, also the corresponding perturbation (i.e. the side effect ∥v^{tr}v^{ut}∥_{1}) can be evaluated. We calculate the frequency of these perturbation values (regardless of the value of τ). The result is shown in Figure 5 (c). We notice that at higher precision, smaller perturbations become slightly more probable: as expected, for high values of P, the algorithm can modulate the inhibition more accurately and therefore reduce the impact on the network, while still satisfying the request on the flux of the objective reaction.
since it holds ${v}_{\text{mod}}^{\text{ut}}={v}_{\text{mod}}^{\text{tr}}(0,0,\dots ,0)$. From this definition, η = 0 means linear behavior and η > 0 nonlinear. Therefore, for each solutions of the screening, we calculate the corresponding η(h) and we analyzed the distribution of its values (still ignoring the parameter τ): the result is shown in panel (d) of Figure 5. It is clear that increasing the value of P the nonlinearity index tends to be higher. It seems that, thanks to the higher precision, the algorithm may exploit more efficiently the nonlinearity property and, by consequence, it can limit the dosage of the drug and consequently reduce the perturbation.
Drug interaction surfaces: three case studies
These calculations has been performed with MOMA based on the L^{1}norm because, as already mentioned, it allows a definition of a linear function in the optimization problem. When compared to the surfaces obtained from the original quadratic formulation of MOMA (last column: panels (c), (f) and (j) of Figure 6), the results from L^{1}norm show some irregularity of the surfaces (which makes the original L^{2} version more reliable) but the main features of the drugdrug interaction are still well described.
Discussion
Optimization is a concept widely used in many scientific fields; for instance, in systems biology, FBA makes use of it for discriminating reaction fluxes in large metabolic networks. Following the same philosophy, in order to cope with more complex situations, multiple optimization criteria can be needed simultaneously leading in some situations, like the one discussed in this paper, to a bilevel optimization problem. The bilevel approach is promising for studying several features and applications of metabolic networks, for instance for identifying metabolic objective functions [23] or for studying perturbations around a nominal optimum [10, 11]. In the context of drug combinatorics, in order to efficiently solve the bilevel optimization, Boolean variables are commonly used in the outer problem. However, this ON/OFF description of the corresponding biological quantities may represent a very rough approximation, as it is the case for the (partial) inhibition induced by drugs acting on the enzymes of a metabolic network.
In order to overcome this limitation, we propose an improvement on the formulation of the bilevel optimization in which a single Boolean variable is replaced by a convex combination of several Boolean quantities: in this manner the convex and linear nature of the problem is preserved and the description of the inhibitory effects becomes more realistic. Since the problem contains Boolean variables, the optimization falls in the Mixed Interger Linear Programming (MILP) class: compared to LP, the NPhard complexity of MILP [22] makes the new algorithm more expensive from a computational point of view. For the tasks at hand (see Figure 3), the algorithm behaves well also for large metabolic networks. The logarithm of the computational time scales linearly with the number of reactions, but with a small slope, so that on average the solution is found in a reasonable computational time, also for networks with around 2500 reactions and for P = 2.
For testing purposes, we run the algorithm on the central carbon metabolism of E.coli screening all reactions. We have found that increasing the number of Boolean variables used in the convex combination (the precision parameter P), it is more likely to find a solution which succeeds on the modulation of the objective reaction (see Figure 5 (a)). In particular, partial inhibitions (i.e. modulations of the dosage of the drugs) are more frequent for multicomponent solutions (panel (b) of Figure 5): this result may be interpreted as a wider possibility, offered by the synergism, to calibrate a drug treatment according to the specific needs. Moreover our computations represent a confirmation on large networks of the expected, but still not verified, higher efficiency of multiple targets drug treatments in presence of partial inhibition [24]. In this perspective, the results show that this approach may also lead to treatments which are more selective (panel (c) of Figure 5 and Table 3).
A possible explanation can be found in the unexpected or hardly predictable drug synergism which are typical of complex systems such as metabolic networks, even in a simplified framework like FBA. In particular, concerning the synergistic interactions between drugs, the analysis done through the drugdrug interaction surface (Figure 6) reveals that nonlinear effects, not explained by superposition of the single drug perturbation, are significant and can be captured and exploited by the method proposed, unlike with a more coarsegrained ON/OFF description. We should mention, that the three case studies presented in Figure 6 do not pretend to have any clinical value: they have been selected only for the purpose of illustrating the method and the advantages it may give in the context of drug synergism and drug reprofiling for reconstructed metabolic networks.
Conclusion
The purpose of this paper is to present a novel algorithm, able to relax the assumption on the variables of a bilevel optimization problem from ON/OFF type to more finegraded description. This setting is of interest in the context of FBA of metabolic networks and in particular in the modulation of the fluxes by means of drugs, capable of reducing (but not suppressing completely) the activity of the metabolic enzymes on which they are acting. With our algorithm, the problem can be formulated as a MILP problem of moderate practical complexity. Indeed we have shown that the algorithm performs well also on large metabolic networks at a more finegraded level of resolution than an ON/OFF representation. Furthermore, it is capable of exploiting with higher efficiency the peculiar nonlinearity which originates from the topology of the network, of finding more selective solutions and, therefore, of offering a larger spectrum of drug combinations. These features become more evident when the modulation required for the objective reaction is itself a fraction (τ ≠ 0) of the nominal flux, rather than a simple complete switch off (τ = 0). Indeed, if a disease corresponds to an anomalous biosynthesis of certain compounds, most often the cure consists in regulating back those fluxes to an healthy range, not to a complete inhibition.
It is worth noting that the problem of drug synergism we presented in this paper must be read as a guiding example for a more general class of situations: indeed, the idea we have proposed for treating bilevel optimization can be applied to any other case which requires a more realistic modeling with respect to the oversimplified ON/OFF description, in biology as well as in all the other fields where LP is already used.
Abbreviations
 FBA:

Flux balance analysis
 MOMA:

Minimization of metabolic adjustment
 LP:

Linear programming
 MILP:

Mixed integer linear programming.
Declarations
Acknowledgments
We would like to thank the Reviewers for their valuable and costructive comments; we thank also Mattia Zampieri (ETH Zurich) for usefull discussion on the method. This work is supported by a grant from MIUR (Ministero dell’Istruzione, Università e Ricerca).
Authors’ Affiliations
References
 Bordbar A, Feist A, UsaiteBlack R, Woodcock J, Palsson B, Famili I: A multitissue type genomescale metabolic network for analysis of wholebody systems physiology. BMC Syst Biol. 2011, 5: 18010.1186/175205095180.PubMed CentralView ArticlePubMedGoogle Scholar
 Folger O, Jerby L, Frezza C, Gottlieb E, Ruppin E, Shlomi T: Predicting selective drug targets in cancer through metabolic networks. Mol Syst Biol. 2011, 7: 501PubMed CentralView ArticlePubMedGoogle Scholar
 Snitkin E, Segrè D, Mackay T: Epistatic interaction maps relative to multiple metabolic phenotypes. PLoS Genetics. 2011, 7 (2): e100129410.1371/journal.pgen.1001294.PubMed CentralView ArticlePubMedGoogle Scholar
 Cornelius S, Lee J, Motter A: Dispensability of Escherichia coli’s latent pathways. Proc Natl Acad Sci. 2011, 108 (8): 312410.1073/pnas.1009772108.PubMed CentralView ArticlePubMedGoogle Scholar
 Raman K, Chandra N: Flux balance analysis of biological systems: applications and challenges. Brief Bioinform. 2009, 10 (4): 435449. 10.1093/bib/bbp011.View ArticlePubMedGoogle Scholar
 Lee J, Gianchandani E, Papin J: Flux balance analysis in the era of metabolomics. Brief Bioinform. 2006, 7 (2): 140150. 10.1093/bib/bbl007.View ArticlePubMedGoogle Scholar
 Jeong H, Tombor B, Albert R, Oltvai Z, Barabási A: The largescale organization of metabolic networks. Nature. 2000, 407 (6804): 651654. 10.1038/35036627.View ArticlePubMedGoogle Scholar
 Deutscher D, Meilijson I, Kupiec M, Ruppin E: Multiple knockout analysis of genetic robustness in the yeast metabolic network. Nat Gen. 2006, 38 (9): 993998. 10.1038/ng1856.View ArticleGoogle Scholar
 Klamt S, Gilles E: Minimal cut sets in biochemical reaction networks. Bioinformatics. 2004, 20 (2): 226234. 10.1093/bioinformatics/btg395.View ArticlePubMedGoogle Scholar
 Burgard A, Pharkya P, Maranas C: OptKnock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng. 2003, 84 (6): 647657. 10.1002/bit.10803.View ArticlePubMedGoogle Scholar
 Kim J, Reed J: OptORF: Optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains. BMC Syst Biol. 2010, 4: 5310.1186/17520509453.PubMed CentralView ArticlePubMedGoogle Scholar
 Jamshidi N, Palsson B: Investigating the metabolic capabilities of Mycobacterium tuberculosis H37Rv using the in silico strain iNJ661 and proposing alternative drug targets. BMC Syst Biol. 2007, 1: 2610.1186/17520509126.PubMed CentralView ArticlePubMedGoogle Scholar
 Huthmacher C, Hoppe A, Bulik S, Holzhutter H: Antimalarial drug targets in Plasmodium falciparum predicted by stagespecific metabolic network analysis. BMC Syst Biol. 2010, 4: 12010.1186/175205094120.PubMed CentralView ArticlePubMedGoogle Scholar
 Bazzani S, Hoppe A, Holzhütter H: Networkbased assessment of the selectivity of metabolic drug targets in Plasmodium falciparum with respect to human liver metabolism. BMC Syst Biol. 2012, 6: 11810.1186/175205096118.PubMed CentralView ArticlePubMedGoogle Scholar
 Chavali A, Hewlett E, Pearson R, Papin J, et al: A metabolic network approach for the identification and prioritization of antimicrobial drug targets. Trends Microbiol. 2012, 20: 113123. 10.1016/j.tim.2011.12.004.PubMed CentralView ArticlePubMedGoogle Scholar
 Suthers P, Zomorrodi A, Maranas C: Genomescale gene/reaction essentiality and synthetic lethality analysis. Mol Syst Biol. 2005, 5: 301Google Scholar
 Facchetti G, Zampieri M, Altafini C: Predicting and characterizing selective multiple drug treatments for metabolic diseases and cancer. BMC Syst Biol. 2012, doi:10.1186/175205096115,Google Scholar
 Palsson B, Varma A: Metabolic capabilities of Escherichia coli II: optimal growth pattern. J Theor Biol. 1993, 165: 503522. 10.1006/jtbi.1993.1203.View ArticleGoogle Scholar
 Segré D, Vitkup D, Church G: Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. 2002, 99 (23): 1511215117. 10.1073/pnas.232349399.PubMed CentralView ArticlePubMedGoogle Scholar
 Dempe S: Foundations of bilevel programming. 2010, New York: Kluwer Academic PublisherGoogle Scholar
 Matoušek J, Gärtner B: Understanding and Using Linear Programming. 2000, Berlin: SpringerGoogle Scholar
 Schrijver A: Theory of Linear and Integer Programming. 1986, New York: John Wiley & SonsGoogle Scholar
 Burgard A, Maranas C: Optimizationbased framework for inferring and testing hypothesized metabolic objective functions. Biotechnol Bioeng. 2003, 82 (6): 670677. 10.1002/bit.10617.View ArticlePubMedGoogle Scholar
 Csermely P, Agoston V, Pongor S: The efficiency of multitarget drugs: the network approach might help drug design. Trends Pharmacol Sci. 2005, 26: 178182. 10.1016/j.tips.2005.02.007.View ArticlePubMedGoogle Scholar
 Peng H, Wen J, Li H, Chang J, Zhou X: Drug inhibition profile prediction for NFκB pathway in multiple myeloma. PloS one. 2011, 6 (3): e1475010.1371/journal.pone.0014750.PubMed CentralView ArticlePubMedGoogle Scholar
 Salvador A: Synergism analysis of biochemical systems. I. Conceptual framework. Math Biosci. 2000, 163 (2): 105129. 10.1016/S00255564(99)000565.View ArticlePubMedGoogle Scholar
 Holzhütter H: The generalized fluxminimization method and its application to metabolic networks affected by enzyme deficiencies. Biosystems. 2006, 83 (2): 98107.View ArticlePubMedGoogle Scholar
 Li Z, Wang R, Zhang X: Twostage flux balance analysis of metabolic networks for drug target identification. BMC Syst Biol. 2011, 5 (Suppl 1): S1110.1186/175205095S1S11.View ArticleGoogle Scholar
 Mahadevan R, Schilling C, et al: The effects of alternate optimal solutions in constraintbased genomescale metabolic models. Metab Eng. 2003, 5 (4): 26410.1016/j.ymben.2003.09.002.View ArticlePubMedGoogle Scholar
 Motter A, Gulbahce N, Almaas E, Barabási A: Predicting synthetic rescues in metabolic networks. Mol Syst Biol. 2008, 4: 168PubMed CentralView ArticlePubMedGoogle Scholar
 Boghigian BA, Armando J, Salas D, Pfeifer BA: Computational identification of gene overexpression targets for metabolic engineering of taxadiene production. Appl Microbiol Biotechnol. 2012, 93 (5): 20632073. 10.1007/s0025301137251.View ArticlePubMedGoogle Scholar
 Snitkin ES, Segrè D: Epistatic interaction maps relative to multiple metabolic phenotypes. PLoS Genet. 2011, 7 (2): e100129410.1371/journal.pgen.1001294.PubMed CentralView ArticlePubMedGoogle Scholar
 Becker S, Feist A, Mo M, Hannum G, Palsson B, Herrgard M: Quantitative prediction of cellular metabolism with constraintbased models: the COBRA Toolbox. Nat Protoc. 2007, 2 (3): 727738. 10.1038/nprot.2007.99.View ArticlePubMedGoogle Scholar
 Ranganathan S, Suthers P, Maranas C: OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions. PLoS Comput Biol. 2010, 6 (4): e100074410.1371/journal.pcbi.1000744.PubMed CentralView ArticlePubMedGoogle Scholar
 Orth J, Fleming R, Palsson B: Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide. EcoSal. 2009,, doi: 10.1128/ecosal.10.2.1,Google Scholar
 Becker S, Palsson B: Genomescale reconstruction of the metabolic network in Staphylococcus aureus N315: an initial draft to the twodimensional annotation. BMC Microbiol. 2005, 5 (1): 810.1186/1471218058.PubMed CentralView ArticlePubMedGoogle Scholar
 Feist A, Scholten JCM, Palsson B, Brockman F, Ideker T: Modeling methanogenesis with a genomescale metabolic reconstruction of Methanosarcina barkeri. Molecular Syst Biol. 2006, 2: doi:10.1038/msb4100046Google Scholar
 Pinchuk G, Hill E, De Ingeniis J, Zhang X, Osterman A, et al: Constraintbased model of Shewanella oneidensis MR1 Metabolism: a tool for data analysis and hypothesis generation. PLoS Comput Biol. 2010, 6: e100082210.1371/journal.pcbi.1000822.PubMed CentralView ArticlePubMedGoogle Scholar
 Edwards JS, Palsson BO: The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proc Natl Acad Sci USA. 2000, 97 (10): 55285533,. 10.1073/pnas.97.10.5528. [http://www.pnas.org/content/97/10/5528.abstract],PubMed CentralView ArticlePubMedGoogle Scholar
 Raghunathan A, Reed J, Shin S, Palsson B, Daefler S: Constraintbased analysis of metabolic capacity of Salmonella typhimurium during hostpathogen interaction. BMC Syst Biol. 2009, 3 (1): 3810.1186/17520509338.PubMed CentralView ArticlePubMedGoogle Scholar
 Knox C, Law V, Jewison T, Liu P, Ly S, Frolkis A, Pon A, Banco K, Mak C, Neveu V, et al: DrugBank 3.0: a comprehensive resource for 'omics’ research on drugs. Nucleic Acids Res. 2011, 39 (suppl 1): D1035—D1041PubMed CentralPubMedGoogle Scholar
 Segré D, De Luna A, Church G, Kishony R: Modular epistasis in yeast metabolism. Nat Genet. 2005, 37 (1): 7783.PubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.