¨ t Berlin Freie Universita Fachbereich Mathematik und Informatik

Master’s Thesis

Thermodynamic Constraints for Metabolic Networks

Author: ¨ ller Arne Mu

Supervisor: Prof. Dr. Alexander Bockmayr

February 27, 2012

Abstract Background For scientific and biotechnology purposes biologists study how cells work and how to manipulate them to derive more information or desired industrial products. To this end, system biologist try to model and simulate cells on the computer. One such modelling approach are genome-scale metabolic networks, with which growth rates can be predicted and knock-out targets can be found. In practice, genome-scale metabolic networks consist of thousands of reactions in one cell (the iJO1366 model of E.coli has 2583 reactions [93]). Thus, it is very hard to understand the role of one single reaction in the context of the whole cell. A large research area established that solely focuses on the interplay of all reactions in the cell, on so called genome-scale metabolic networks. Flux Balance Analysis (FBA) is one of the most often used methods to analyze the capabilities of a metabolic network. It is based on the steady-state assumption. This means that we assume that every internal metabolite that is produced must also be consumed at the same rate. In particular, FBA does not require knowledge on kinetic parameters of the reactions. It is an optimization based approach that can be performed efficiently by solving a linear program. FBA can also give very precise results. For example, it is plausible that species evolved to grow as fast as possible and indeed, FBA predicts growth rates very well [35]. In addition to the growth rate, FBA also tells us an optimal flux (pathway) through the reactions, i.e., it tells us how fast each reaction must transform its reactants into its products. However, the flux that achieves optimal growth rate is usually not unique. This is a major drawback, because the production rate of byproducts can depend on the flux that is used by the cell in vivo [71]. Since there usually is not only one such optimal flux, the computed flux may represent the flux in vivo only poorly. Hence, we are not only interested in one optimal flux, but in the whole space of optimal fluxes. Among other methods, flux variability analysis (FVA) is one such method that can be used to understand the flux space [81]. If flux variability analysis is performed without thermodynamic constraints, obviously unrealistic effects can be witnessed, like that the range of flux through a reaction is unbounded. This effect is caused by internal cycles in the network. These cycles can carry unbounded flux, since they do not take up any nutrients. The most basic thermodynamic constraints already prohibit such internal cycles [9].

i

Results In this thesis, I will show that already these most basic thermodynamic constraints make thermodynamic constrained FBA NP-hard and thus computationally difficult for large networks. Because of this complexity result, it is important to develop a good understanding of the mathematical properties of thermodynamic constraints and properties of common genome-scale networks that appear in practice. In the literature many different (mostly physically motivated) notions of thermodynamically feasible fluxes are used. Not all of these definitions are equivalent. Hence, we will analyze definitions and properties of thermodynamically feasible fluxes. This will lead to computational methods for checking feasibility and existence theorems. FBA is only one possible analysis method for metabolic networks, hence we will analyze the effect of thermodynamic constraints for sampling methods, elementary modes, and optimization based approaches. Since optimization based approaches like FBA and FVA work also well on genomescale networks, the integration of thermodynamic constraints into these methods is challenging and most of this thesis will be devoted to this topic. Currently, thermodynamic constrained FBA is performed using mixed integer linear programming (MILP). Although the networks consist of thousands of reactions and thus the MILPformulation has thousands of boolean variables, current MILP solvers are able to solve these networks quickly. We will analyze this MILP formulation regarding its LP-relaxation, possible cuts and reformulations. We will see that most metabolic networks that are used in practice have few internal circuits. This property is used to derive a parametrized tractability result that is also applicable in practice. This theoretical approach is realized by a constraint programming approach that allows us to perform thermodynamically constrained FVA (on E. coli iAF1260 ) in less than 5 minutes compared to 5 hours of the MILP-approach. The implementation can be downloaded from http://page.mi.fu-berlin.de/arnem/fast-tfva.html.

ii

Zusammenfassung Hintergrund Um neue wissenschaftliche Erkenntnisse zu gewinnen und neue Medikamente, Bio-Treibstoffe etc. zu entwickeln, sind Biologen daran interessiert die Funktionsweise von Zellen zu verstehen und gezielt zu manipulieren. Zu diesem Zweck haben sich Computermodelle der Zellen etabliert. Eine solche Art von Modell sind metabolische Netzwerke. Diese erlauben die Vorhersage von Zellwachstum und den Effekt von Gen Knock-Outs. Um m¨ oglichst genaue Ergebnisse zu erzielen, werden metabolische Netzwerke erstellt, die m¨oglichst alle Reaktionen, die in der Zelle stattfinden, enthalten. Selbst f¨ ur Bakterien, wie E. coli, erh¨ alt man Netzwerke mit mehreren tausend Reaktionen. Eine manuelle Auswertung der Bedeutung von einzelnen Reaktionen im Kontext der ganzen Zelle ist somit fast unm¨oglich. Somit wurden und werden Computer gest¨ utzte Analyse-Methoden entwickelt. Fluss Balance Analyse (FBA) ist eine der meist genutzten Methoden um die F¨ahigkeiten eines metabolischen Netzwerkes zu analysieren. Diese Methode basiert auf der Annahme von Flusserhaltung, d.h. jeder Stoff innerhalb der Zelle muss in der gleichen Menge produziert wie konsumiert werden. Dies hat den Vorteil, dass kinetische Parameter, die meist unbekannt sind, nicht ben¨otigt werden. Mittels linearer Programmierung kann FBA effizient gel¨ost werden und produziert sogar sehr gute Ergebnisse, wie z.B. die Vorhersage von Wachstumsraten von E. coli gezeigt hat [35]. Außerdem sagt uns FBA in welchen Geschwindigkeiten die einzelnen Reaktionen ablaufen m¨ ussen. Nun ist es aber so, dass die Flussgeschwindigkeiten in der Regel nicht eindeutig bestimmt sind. Das ist ein großer Nachteil, da somit die Produktionsrate von Nebenprodukten nur schwerlich vorhergesagt werden kann [71]. Somit ist man nicht nur an einem Fluss interessiert, sondern an Eigenschaften aller Fl¨ usse. Fluss Variabilit¨ ats Analyse (FVA) ist eine Methode, die daf¨ ur genutzt werden kann [81]. Damit werden die maximalen und minimalen Fl¨ usse durch die Reaktionen bestimmt. Wenn man dies aber tut ohne thermodynamische Nebenbedingungen hinzuzuf¨ ugen, wie das u ¨blicher Weise der Fall ist, dann kommt es vor, dass die Ergebnisse der Analyse offensichtlich nicht realistisch sind. Durch innere Kreise in dem Netzwerk kann es n¨amlich zu unbeschr¨ankten Zirkulationen kommen. Schon die einfachsten aller thermodynamischen Nebenbedingungen jedoch verbieten diese Kreise [9].

iii

Ergebnisse In dieser Arbeit werde ich zeigen, dass allerdings schon die einfachsten dieser thermodynamischen Nebenbedingungen daf¨ ur sorgen, dass FBA NP-schwer wird. Dadurch werden diese, aus biologischer Sicht einfachen Nebenbedingungen, berechnungstechnisch kompliziert. Um dennoch effiziente Algorithmen zu entwickeln, ist es um so wichtiger diese Nebenbedinungen genau zu verstehen und auch die Eigenschaften von metabolischen Netzwerken, wie sie denn in der Praxis auftauchen, zu kennen. Die Problematik beginnt schon dabei, dass in der Literatur unterschiedliche Definitionen von thermodynamisch g¨ ultigen Fl¨ ussen verwendet werden. Nicht alle diese Definitionen sind ¨aquivalent. Deswegen werden wir diese Definitionen uns genauer anschauen und Eigenschaften beweisen, die dann auch dazu genutzt werden k¨ onnen um effizient zu pr¨ ufen, ob ein gegebener Fluss thermodynamisch g¨ ultig ist. FBA ist nur eine von mehreren Analyse-Methoden f¨ ur metabolische Netzwerke. Deswegen werden wir neben Methoden wie FBA und FVA auch Sampling und Elementarmodi Analyse betrachten. Da FBA und FVA aber die einzigen Methoden sind, die auch gut auf große Netzwerke anwendbar sind, haben ich mich darauf konzentriert, wie man thermodynamischen Nebenbedingungen in FBA und FVA integrieren kann. Die meisten heutzutage g¨ angigen Methoden f¨ ur FBA mit thermodynamischen Nebenbedinungen benutzen gemischt ganzahlige Programmierung (MILP). Trotz der Gr¨oße k¨onnen metabolische Netzwerke mit mehreren tausend Reaktionen (und fast entsprechend vielen ganzzahligen Variablen im MILP) schnell von g¨ angigen MILP-L¨osern gel¨ost werden. Wir werden diese MILPFormulierung auf ihre LP-Relaxierung, Schnittebenen-Verfahren und Reformulierungen hin analysieren. Viele metabolische Netzwerke, die in der Praxis benutzt werden, haben nur sehr wenige innere Kreise. Ich werde zeigen, dass wenn die Anzahl der Kreise beschr¨ankt ist, man das Problem in polynomieller Zeit l¨ osen kann. Dieses theoretische Ergebnis setze ich dann mittels Constraint Programmierung um. Mit dem daraus entstandenen Algorithmus ist es m¨oglich FVA mit thermodynamischen Nebenbedinungen (auf E. coli iAF1260 ) anstelle von 5 Stunden in weniger als 5 Minuten zu l¨ osen. Die Implementierung steht auf http://page.mi.fu-berlin.de/arnem/ fast-tfva.html zum Download bereit.

iv

Acknowledgments When I finished school, I wasn’t sure what I liked best: Mathematics or Computer Science. I liked to design algorithms, in particular because these things could be used to compute actual stuff and thus get results that would otherwise be inobtainable. On the other hand, I considered simple programming without much theoretic background dull and boring. So this lead me to combinatorial optimization and mixed integer programming. Since an internship in grade 10 at Max Planck Institute for Molecular Biology, I also developed a fascination on genetics and molecular biology. So, first of all I want to thank my supervisor Alexander Bockmayr for giving me the opportunity to work in the highly fascinating and interesting intersection of mathematics, computer science and biology. I also want to thank him for pointing me to the topic of thermodynamic constraints in metabolic networks with its rich mathematical background, which I will also continue studying after finishing my master’s thesis. In that context I also want to mention Daniel A. Beard and Hong Qian, whose work was the starting point of my research. In particular the article where they used oriented matroids to describe thermodynamic feasibility was one big source of fascination and motivation. I also want to thank the members of my workgroup. In particular I want to thank Heike Siebert and Xingyi Yang with whom I shared an office for a study-friendly environment; Yaron Goldstein with whom I had several discussions on thermodynamic constraints, because they do not satisfy one of the lattice axioms he is working with to describe flows in metabolic networks. I also want to thank Shahrad Jamshidi, Hannes Klarner and Aljoscha Palinkas who accompanied me at lunchtime very often, where we had loads of fun and interesting discussions. Next to my workgroup the influence of the Berlin Mathematical School must not be forgotten. There I met many new international students, several of them I now count to my friends. It also motivated me to take many classes at TU Berlin and thus continue the path that I started with a lecture by Andreas Paffenholz on “Geometry and Optimization”. There, Rolf M¨ohring, Andreas Bley, Martin Skutella and Martin Gr¨otschel showed me the interesting facets of integer and mixed integer programming in lectures, seminars and a highly interesting summer school on network flows organized by Martin Skutella. ´ Furthermore, I want to thank my friends and fellow students Agnes Cseh, Giovanni Conforti, Julie Meißner, Han Cheng Lie, Christoph Hansknecht, Wilhelm Neubert, Matthias Frey, Paul Podlech, and Gerrit Gruben together with whom I attended several lectures and seminars, and who accompanied me (and probably also had some influence) on my specialization path. In this list of friends I also do not want to miss Atul Shekhar, Ahmad Afuni, Andrian Gonzales Casanova-Soberon with whom I unfortunately never attended classes together but spend many nice evenings with.

v

I also want to thank the other two members of Piratenfraktion Steglitz-Zehlendorf, Eric L¨ uders and Georg Boroviczeny, with whom I spend many long evenings learning how to do communal politics as elected Bezirksverordnete. I also want to thank all the other members of the Piratenpartei from Steglitz-Zehlendorf in their encouragement and support for this new and compared to mathematics entirely different kind of work that I have gotten myself into. Last but not least, I want to thank my parents and my brother for supporting me in my path and for many interesting mathematical discussions during supper. Berlin, February 27, 2012 Arne M¨ uller

vi

Contents Abstract

i

Zusammenfassung

iii

Acknowledgments

v

1 Introduction

1

1.1

Metabolic Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Steady-State Assumption and Thermodynamic Constraints . . . . . . . . . . . .

4

1.3

Mathematical Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4

Structure of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2 Thermodynamically Feasible Pathways 2.1

2.2

Without Information on Metabolite Potentials

9 . . . . . . . . . . . . . . . . . . .

11

2.1.1

Oriented Matroids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.2

Internal Cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

With Information on Metabolite Potentials . . . . . . . . . . . . . . . . . . . . .

19

3 Flux Space

25

3.1

Structural Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.2

Thermodynamically Feasible Flux through a Given (Internal) Reaction . . . . . .

27

3.3

Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

4 Generating Vectors, Extreme Rays and Elementary Modes

31

4.1

Enumeration of Elementary Modes . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.2

Elementary Modes and Oriented Matroids . . . . . . . . . . . . . . . . . . . . . .

33

4.3

Thermodynamically Infeasible Elementary Modes . . . . . . . . . . . . . . . . . .

36

5 Optimization Based Approaches

39 vii

5.1

Thermodynamically Infeasible Optimal Solutions . . . . . . . . . . . . . . . . . .

40

5.2

General Difficulty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.3

Nonlinear Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

5.4

Mixed Integer Linear Programming . . . . . . . . . . . . . . . . . . . . . . . . . .

45

5.4.1

Mixed Integer Linear Programming (MILP) formulations . . . . . . . . .

45

5.4.2

Quality of the LP Relaxation . . . . . . . . . . . . . . . . . . . . . . . . .

48

5.4.3

Cutting Planes for MILP formulations . . . . . . . . . . . . . . . . . . . .

50

Constraint Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

5.5.1

A Heuristic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

5.5.2

Parametrized Complexity . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

Performance Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

5.6.1

Optimization on Single Reactions . . . . . . . . . . . . . . . . . . . . . . .

66

5.6.2

Flux Variability Analysis on Networks of Different Sizes . . . . . . . . . .

70

5.5

5.6

6 Future Work

75

6.1

Optimization over Fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

6.2

Optimization over Potentials, Integrated Models . . . . . . . . . . . . . . . . . .

77

6.3

Conclusion

78

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Appendix

79

Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Index

81

Bibliography

83

viii

Chapter 1

Introduction Living organisms are very complex systems. In the past biologists tried to deal with the problem by looking at single enzymes or very special subsystems only. These methods were ignoring the interactions between the subsystems and were never able to give a whole-cell system wide view. With the increasing amount of sequenced genomes, experimental high-throughput analyses and computational technology, the field of systems biology emerged. System biologists do not try to understand single enzymes, but the functionality of the whole cell. Many system biologists work on the reconstruction of metabolic networks. Metabolic networks help to understand physiological processes, which in turn help to bring enzymes in a whole cell context and annotate them with functional properties [72]. But not only deeper understanding of the functionality of the cell can be obtained [122, 76], there are als many applications in medicine [11] and bioengineering [121]. With the advent of wholly sequenced genomes and annotated genome-databases, the reconstruction of genome-scale metabolic networks becomes easier and easier [34, 47, 88, 129]. Genome-scale models give us the possibility to analyze the behavior of the cell as a whole on the computer. This is why biologists also talk in terms of in silico experiments in contrast to in vitro and in vivo experiments that are carried out in the wet lab or on living organisms respectively. Since many genes can be associated to enzymes and the enzymes to the biochemical reactions they catalyze, the effects of gene knock-outs can be analyzed without even performing mutations on real cells [23, 118]. Also symbiotic and parasite metabolism can be analyzed, which is very difficult in vitro [128]. But, genome scale models are huge, for example the iAF1260 model of the bacterium Escherichia coli contains 2077 reactions and 1668 metabolites [36]. Hence, mathematical and computational methods are necessary to analyze those networks [125].

1.1

Metabolic Networks

Before we can talk about methods for metabolic models, we need to understand what a metabolic model is. A metabolic network consists of a set of metabolites M and a set of reactions R. But this is already everything all metabolic models have in common. The amount of information on the reactions and metabolites may vary depending on purpose of the network. Information you may encounter in metabolic networks are 1

CHAPTER 1. INTRODUCTION

• Stoichiometries (stoichiometric matrix S ∈ RM×R ): Every metabolic model contains the stoichiometries of the reactions, but the quality of the stoichiometries may differ. In some cases for example mass conservation is not guaranteed by the reaction. This can be caused when similar metabolites were lumped together, the production/consumption of H2 O or electrons were ignored [50, 41]. The stoichiometric matrix stores the stoichiometric coefficients. See Figure 1.1 for a small real-world example and 1.2 for an even smaller made-up example.

Figure 1.1: Network of E. coli core energy metabolism. Figure taken from [9] • Exchange Reactions E ⊂ R: Exchange reactions model the inflow of nutrients and products in and out of the network. Hence, these reactions do not satisfy mass balance on purpose (else there couldn’t be any inflow of nutrients). The reactions that are not exchange reactions are called internal reactions I := R \ E.

One very important exchange reaction is the biomass reaction. It is possible to measure the types and amounts of amino acids the cell is build out of. If the cell grows, all those amino acids need to be produced in exactly that ratio, else necessary amino acids may be missing and essentially proteins cannot be build. The biomass reaction consumes the amino-acids in exactly the necessary ratio. Hence, the biomass reaction measures how fast the cell grows [37].

• Metabolite Concentrations c ∈ RM : It is very hard to measure (or estimate) the concentrations of all metabolites, but the area of metabolomics shows great advances and more and more metabolite concentration can be measured [30, 33, 55, 38]. This is valuable information from which you can estimate reaction directions and sometimes even reaction rates. • Flux J ∈ RR : The reaction rates are often called flux. The notion is analogous to flow in graphs. As with metabolite concentrations, the flux rates are often unknown and will often also depend on the environment. Some flux rates can be measured using atom-labeling experiments [134, 107]. • Kinetics f : RM → RR : To compute flux of reactions from metabolite concentrations properly, information on the kinetics of the reaction are needed. The kinetics are usually 2

1.1. METABOLIC NETWORKS

a A −1 B  5 C  0 S= D  0 E 0 F 0 

b −1 1 2 0 0 0

c 0 1 0 −1 0 0

d 0 −1 0 3 −1 0

e 0 0 −1 0 1 0

f 0 0 0 −1 0 1

g 0 0 0 0 −1 2

x1 1 0 0 0 0 0

x2 0 1 0 0 0 0

x3  0 0   0   0   0  −1

Figure 1.2: Example of a metabolic network with corresponding stoichiometric matrix S. The exchange reactions are denoted by x1 , x2 , x3 and are simply indicated by arrows in the Petri-net type of drawing. Metabolites are drawn by ellipses and reactions by rectangles. Numbers on the arrows are the stoichiometric coefficients of the metabolite in the reaction connected by the arrow. If there is no number on an arrow, the stoichiometric coefficient is 1. Reaction g is split into two reactions to indicate that it is reversible. functions that basically take metabolite concentrations as input, and compute the reaction rate. Good kinetic information is the rarest of all. Usually there are a few default types of kinetic functions that are adapted using constants that specify enzyme concentration, temperature, pH etc [86]. Because kinetic data is so hard to obtain, methods were developed that also work with reduced kinetic information [123]. • Regulatory Information: Not every reaction that can happen, happens. For example, it may be the case that regulatory control of the cell inhibits that the needed enzyme is produced. Available transciptory data of gene expression is used to integrate these behaviors in the model [29, 130]. If an enzyme that catalyzes a reaction is not produced, this means that the reaction cannot carry any flux. In this case, this kind of information can simply be integrated into the model by deleting the corresponding reaction [99]. If kinetic information is available, a so called kinetic model can be build. Since the reactions modify metabolite concentrations, the first derivative of the concentrations is determined by flux: ∂ c = SJ = Sf (c) ∂t

(1.1)

Hence we obtain an ODE describing the evolution and growth of the cell. As already mentioned, good kinetic parameters are very hard to obtain. Hence constraint-based methods were developed that do not need kinetic parameters. 3

CHAPTER 1. INTRODUCTION

1.2

Steady-State Assumption and Thermodynamic Constraints

Constraint-based methods require the additional steady state assumption, which is motivated as follows [99, 125]: If the environment does not change, the dynamic system (1.1) will reach a fixed point. In cases where the environment changes only slowly (compared to flux speeds), the error will not be very big. But the effect of this additional assumption is enormous, since Equation 1.1 simplifies to 0 = SJ,

(1.2)

which is also called the steady state assumption or flux conservation. It can be considered as the metabolic equivalent to the flow conservation constraint in graphical flow problems. You simply need to replace the stoichiometric matrix by the incidence matrix of a digraph. See Figure 1.3 for an example.

Figure 1.3: An example of a steady state flux in a metabolic network. The reactions that carry flux are indicated in red. The amount of flux through the reactions can be read off from the red numbers next to the reactions. The main difference of constraint-based methods is, that we are not interested in the evolution of one flux vector but on the set of all feasible fluxes under the given constraints (i.e. at steadystate) [99, 112]. In this context additional information on the biological feasible fluxes becomes vital. Some reactions for example can only proceed in one direction; hence, sign constraints on J are added. Sometimes it is also possible to measure how fast one reaction can proceed; this leads to upper (and lower) bounds on the feasible fluxes. One additional constraint that will be in the focus of this thesis is thermodynamic feasibility, which was first discussed by Oster in 1971 [95]. To understand thermodynamically feasible fluxes, we will have to include concentrations in a certain sense back into our formulation. This is no restriction on applicability, since it does not mean that we need to measure concentrations. As with flux, we can treat metabolite concentrations as unknown variables. A flux is called thermodynamically feasible if it does not violate the second law of thermodynamics. The second law of thermodynamics states that a reaction carries flux if and only if it reduces Gibbs free energy [10, 4, 101]. It is the same law that prohibits electrical current to go around a cycle if no energy source is attached [103]. We will see that this cycle-correspondence also holds for metabolic networks; for example the steady state flux shown in Figure 1.4 is not thermodynamically feasible because it contains an internal cycle. The reduction of Gibbs free energy is formulated using potential differences. Every metabolite 4

1.2. STEADY-STATE ASSUMPTION AND THERMODYNAMIC CONSTRAINTS

Figure 1.4: The steady state flux shown in this metabolic network (again in red) is not thermodynamically feasible, because it contains an internal cycle (dark red) i ∈ M in the network has a biochemical potential µi . This biochemical potential can usually be computed from its concentration ci by µi = µ0i − RT ln(ci ),

(1.3)

where R is the gas constant, T is temperature, and µ0i is the equilibrium potential, which is different for each type of metabolite [101, 95]. The second law of thermodynamics can then be formulated using the chemical potentials as for every internal reaction r ∈ I

∆µr Jr < 0 or ∆µr = 0 = Jr

(1.4)

where ∆µr := µSr is the potential difference induced by reaction r ∈ I (Sr denotes the r.th column of S). This formulation simply states, that a reaction carries positive (resp. negative) flux, if and only if it has negative (resp. positive) potential difference. Note, that lumped metabolites (because they are not interesting by-products etc.) may significantly perturb potential differences. Thus, excellent stoichiometric information is necessary for these constraints to be applicable. But this is already everything that is needed to add these additional constraints. This is the strength of thermodynamic constraints. Without much more additional information (except good quality stoichiometries), we can obtain more realistic results. In particular, we do not need to know information on reaction kinetics. Since these thermodynamic constraints do not require much more additional information, they can easily be added to genome-scale metabolic networks [109]. But thermodynamic constraints are not only used to obtain more consistent flux results [109, 36, 41, 100], but also to derive information on the feasible concentration space and putative regulatory sites [11, 76, 138, 27]. Systems biologists use several different constraint-based methods to analyze metabolic networks. These methods also have different areas of application. For example sampling methods are used to understand the set of all feasible fluxes and relationships between reactions using statistical methods. In elementary flux mode analysis (EFM, or sometimes EM) the flux space is analyzed by its extreme pathways. Extreme pathways contain a minimal number of reactions and thus are easier understandable by humans. In flux balance analysis (FBA) [94] on the other hand, one is interested in a pathway with optimal production of some target metabolite like for example biomass. This is a very common method, since biologists assume that by evolution the species evolves to the most efficient and thus competitive form. This assumption is also supported by experimental results, like on the growth rate of E. coli [35]. Optimization based approaches are also used if we want to determine the maximal flux that is possible through a given reaction. This information can then be used to analyze knock-out candidates. 5

CHAPTER 1. INTRODUCTION

Although thermodynamic constraints are easy to integrate from a biological perspective, we will see that it is hard to integrate thermodynamic constraints into these analysis methods from a computational point of view.

1.3

Mathematical Context

Analyzing thermodynamically feasible fluxes is not only interesting from a biological perspective. First of all metabolic networks are nothing else than directed, weighted hypergraphs. Hence, we can use Petri-nets [104, 133] and oriented matroids [17, 9] to describe these networks. In particular, we will look at oriented matroids more closely. The key notion of oriented matroids are circuits and cycles. Every steady state flux through a metabolic network can be seen as a cycle in an oriented matroid. In particular, elementary modes (one of the favorite tools of system biologists) are nothing else than circuits of an oriented matroid. One indicator of thermodynamic infeasibility will be the existence of internal circuits - again, we can use oriented matroid theory to analyze this. The set of steady state fluxes is described by linear equalities. In a popular analysis method, flux balance analysis (FBA), the only additional constraint are bounds on the fluxes; under these constraints, using a linear objective function the flux through a set of reactions is maximized. It is easy to see that these constraints can be formulated as an LP and thus solved using linear programming. On the other hand, we can also formulate every LP as an FBA problem. If we are given an LP in standard form max cx : Ax = b, x ≥ 0, we can add exchange fluxes that have to carry flux corresponding to the right hand side vector b (here, I denotes the identity matrix): max cx : Ax − Iy = 0, x ≥ 0, y = b
Clearly, this is an FBA problem of a metabolic network with stoichiometric matrix A −I . Observe that the metabolic network obtained this way may contain internal cycles that generate, respectively consume flux. This is a property that proper metabolic networks should not have. However, if metabolites were omitted in the reconstruction of a metabolic network, this may also occur in real world networks, too. This can be fixed by adding additional boundary metabolites (i.e. metabolites with an exchange reaction for the metabolite) to every reaction that participates in an internal cycle that consumes or produces flux. Thus, FBA is equivalent to solving LPs and by analyzing fluxes in metabolic networks we may also get new insights on LP-theory, because we are looking at LPs from a different perspective. Computing thermodynamically feasible flows is NP-hard (Theorem 5). Thus, in general, we cannot solve the problem by simply solving a couple of LPs. The most commonly used thermodynamically constrained analysis methods involve nonlinear programming and mixed integer programming. Although nonlinear programming formulations are flexible, non-convex domains make analyzing these methods difficult. In this thesis, we will mostly focus on mixed integer linear programming and constraint programming. Even on large scale networks, current MILP solvers are able to solve thermodynamically constrained problems in a reasonable amount of time, although straight-forward analysis methods indicate otherwise. Parametrized complexity theory is needed to capture these practical results also in a theoretical context. These theoretical results can then be used to build specialized algorithms that exploit common properties of metabolic networks. 6

1.4. STRUCTURE OF THIS THESIS

1.4

Structure of this Thesis

First of all, we want to understand the properties of thermodynamically feasible fluxes. Hence, Chapter 2 will introduce fundamental definitions and notation. We will prove equivalent definitions of thermodynamic feasibility and discuss methods for testing thermodynamic feasibility. It is highly recommended to read Chapter 2 before continuing with the other chapters. System biologists, however, do not only want to know whether a certain pathway is feasible or not, but they want to understand how the set of all feasible pathways looks like. This is also important for computational methods. In Chapter 3, we will analyze topological properties and convexity of the flux space. The results of the convexity analysis will lead to an NP-hardness result. To motivate more sophisticated analysis methods, we will also discuss sampling as a very direct method for analyzing the flux space and the problems arising from it. In Chapter 4, we will then talk about one of the System Biologists favorite method: Elementary Modes. In particular, we will see a direct link between thermodynamic constraints and elementary modes. The main part will be Chapter 5, where the effect of thermodynamic constraints of optimization methods will be analyzed. Several solution strategies from nonlinear programming over mixed integer programming and constraint programming will be discussed. The Chapter will end with a comparison of a constraint programming approach to other current implementations. Chapters 3, 4, and 5 have only few interconnections and it should be possible to read them independently of each other. Finally, in Chapter 6, we will summarize the results and discuss open problems as well as additional demands from system biologists.

7

CHAPTER 1. INTRODUCTION

8

Chapter 2

Thermodynamically Feasible Pathways In the introduction, we gave a brief intuition of thermodynamically feasible fluxes. Since the thermodynamic feasibility of a flux does not only depend on the flux itself but also on metabolite concentrations, the term “thermodynamically feasible flux” per se is not well defined. We will now develop a clean mathematical definition of thermodynamically feasible fluxes. This definition will then be used to investigate thermodynamically feasible fluxes in more detail. We will start off with the definition of a metabolic network in the mathematical sense. Definition 1 (Metabolic network) A three-tuple N = (M, R, S) , S ∈ RM×R is called a metabolic network. Its metabolites are denoted by M, the reactions by R, and the stoichiometric matrix by S. 2 The definition is similar to the one of a graph with the exception that we encode the incidence information via the stoichiometric matrix and not in the reactions themselves, see Figure 1.2. To work properly with the stoichiometric matrix, we will write S∗r to denote the column corresponding to reaction r and Sm∗ to denote the row corresponding to metabolite m respectively. Smr then denotes the entry corresponding to reaction r and metabolite m. If it is clear from the context, we will omit the ’*’ and simply write Sr to denote the r.th column or the r.th row of S respectively. We will also use sets of indices to denote submatrices; SI for example will denote the stoichiometries of all internal reactions I. We will also use similar notation on vectors, like the flux vector J; there, Jr denotes the flux through reaction r. In the Appendix Notation you will find a summary of all notation used in this thesis. Through out the thesis we will work with metabolic networks that satisfy the following two assumptions: Assumption 1 (S does not contain zero columns) In the whole thesis, we will assume that the stoichiometric matrix S of any metabolic network does not contain any zero columns (i.e. columns that contain only zeros). This is not a restriction for real-world applications, since zero columns correspond to reactions that do not involve any metabolites. Assumption 2 (Finite Models) We will also assume that R, M are finite. It may be strange to mention this explicitly, but there are works by Hatzimanikatis et al. [60] that implicitly define metabolites and reactions. In such a setting this assumption may not be true. 9

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

The first property we derived was the steady state assumption, which lead to the notion of steady state fluxes [99, 125] . An example can be seen in Figure 1.3. Definition 2 (Steady state flux) Given a metabolic network N = (M, R, S), we call a flux vector J ∈ RR a steady state flux in N if SJ = 0 In the following, we will partition the set of reactions into internal reactions I and exchange reactions E. The internal reactions will be subject to thermodynamic constraints, but the exchange reactions will be not. We can now define thermodynamically feasible fluxes with respect to given metabolite concentrations by adding the condition that a reaction carries flux if and only if it has negative potential difference [10, 4, 101]: Definition 3 (Strongly thermo. feasible flux w.r.t. metabolite concentrations) ˙ S) and metabolite potentials µ ∈ RM , we call a Given a metabolic network N = (M, R = I ∪E, R flux vector J ∈ R strongly thermodynamically feasible in N with respect to µ if • J is a steady state flux in N . • µS∗r Jr < 0 or µS∗r = 0 = Jr for all r ∈ I.

2

Figure 2.2: In this example, there exists no strongly thermodynamically feasible flux with these potentials, since the metabolite with potential 5 would have to be a sink, which is not allowed. But, there exists a weakly thermodynamically feasible flux, which carries only flux on the continuous arcs and none on the dashed ones (strongly thermodynamics would force flux on the red, dashed arcs).

Figure 2.1: According to the potentials on the nodes and the induced potential differences on the reactions, every reaction must carry flux to be strongly thermodynamically feasible. This is possible, since for every metabolite that is produced, there also exists a reaction that consumes it.

See Figure 2.1 for an example of metabolite potentials that allow a strongly thermodynamically feasible flux. For practical reasons, it is useful to relax the second condition slightly such that reactions do not have to carry flux if they have non-zero potential difference, but can carry flux (see Figure 2.2). For example if a reaction is catalyzed by an enzyme and the enzyme is inhibited, the reaction will not be able to carry flux although there may be non-zero potential difference. Another, mathematically motivated, reason can be observed in Section 3. Hence, we will call the thermodynamically feasible fluxes w.r.t. Definition 3 strongly thermodynamically feasible in contrast to weakly thermodynamically feasible fluxes that only have to satisfy the relaxed condition. For the definition of weakly thermodynamically feasible fluxes, we will follow the definition proposed by Beard et al. [9]. Since this will actually be the definition we will work with the most of the time in this thesis, we will often simply omit the term “weakly”: 10

2.1. WITHOUT INFORMATION ON METABOLITE POTENTIALS

Definition 4 ((Weakly) thermo. feasible flux w.r.t. metabolite concentrations) ˙ S) and metabolite potentials µ ∈ RM , we call a Given a metabolic network N = (M, R = I ∪E, R flux vector J ∈ R thermodynamically feasible with respect to µ if • J is a steady state flux in N • µS∗r Jr < 0 or Jr = 0 for all r ∈ I.

2

In many applications not all metabolite potentials are known. The potentials may actually also change when regulatory mechanisms modify metabolite concentrations. Hence, we will first consider the case when nothing on the metabolite potentials is known and define and analyze properties of thermodynamically feasible fluxes in that context. Since it is often the case that some information is available on metabolite potentials, we will also give a corresponding definition for that scenario. In addition, we will generalize the results from the first part to the more general case.

2.1

Without Information on Metabolite Potentials

If we have no knowledge on metabolite potentials, it seems as we could not derive any information on flux directions. But we can; although we do not know anything on the potentials, we have some information on the potential differences, which can be used. Definition 5 (Thermodynamically feasible flux) Given a metabolic network N = (M, R = ˙ S) we call a flux vector J ∈ RR thermodynamically feasible if there exist metabolite poI ∪E, tentials µ ∈ RM such that J is thermodynamically feasible w.r.t. µ. 2 Sometimes, when biologists apply a method called flux balance analysis (FBA) on internal reactions, they observe a very unrealistic result. FBA tries to find a maximal steady state flux through the metabolic network. See Section 5 for more information on FBA. These unrealistic results have the form that they contain a large (nearly unbounded) internal circulation. Such an internal circulation is unrealistic in the same sense as that electrical current does not flow in a circuit if you do not attach a battery. It turns out that these unrealistic fluxes are not thermodynamically feasible by Definition 5. This correspondence was probably known to systems biologists for several years but Beard et al. [9] were the first who investigated the problem from a mathematical perspective using oriented matroids. Oriented matroids were originally developed by Bland [18] as a generalization of linear programming. Independently to Bland, Las Vergnas developed oriented matroids as a generalization of directed graphs [79]. This generalization is similar to the one by matroids, which can be seen as a generalization of graphs. In particular, metabolic networks can be described by oriented matroids, which was first done by Oliveira et al. in 2001 [89]. We will now enjoy a small excursion into oriented matroid theory since the notion of circuit and cocircuit are native to oriented matroids or matroids in general. Although Beard et al. [9] used oriented matroids to prove the correspondence between internal circuits and thermodynamically infeasible fluxes, we will also give an alternative proof that simply uses LP duality theory. 11

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

2.1.1

Oriented Matroids

Oriented matroids can be seen as matroids with orientation information. As with matroids, oriented matroids describe the metabolic networks by the reactions R, which are called elements, and cycles, which are signed sets of reactions. In particular, the notion of metabolite is lost in oriented matroid theory. The book by Ziegler has a very nice introduction to oriented matroids in chapter 6 [141]. In contrast to matroids, oriented matroids also store information on the sign of the reactions that form the circuit. For example if (−2, 3, 4, −5, 1) is a cycle (because it is in the nullspace of S), the corresponding matroidal vector is (−, +, +, −, +). Hence, the vectors of oriented matroids cannot be simply represented by subsets of R, but we need to use signed subsets. Notation Let E be a set. Based on the notation used in [17], we call C = (C + , C − ) with ˙ − ⊆ E a signed subset of E. Alternatively, we also use the incidence vector notation for C + ∪C C, i.e. C ∈ {−, 0, +}E , where Ci = − ⇔ i ∈ C − and Ci = + ⇔ i ∈ C + for all i ∈ E. These two forms of notation will be used interchangeably. Since oriented matroids are a generalization of ˙ − will denote the matroids, we sometimes are not interested in the signs and hence, C := C + ∪C support of C. We will now give an abstract definition of oriented matroids and state some of the main theorems. See the book by Bj¨ orner et al. [17] for a comprehensive survey on oriented matroids. Definition 6 (Oriented Matroid) A tuple M = (E, C) with elements E and circuits C ⊆ {−, 0, +}E is called an oriented matroid if the following circuit axioms are satisfied: C0 ∅ ∈ 6 C, where ∅, respectively 0, are shorthand notations for (∅, ∅) in the context of signed sets. C1 C = −C, where −C := {−C : C ∈ C} C2 for all X, Y ∈ C if X ⊆ Y , then X = Y or X = −Y . C3 for all X, Y ∈ C with X 6= −Y and e ∈ X + ∩ Y − , there is a Z ∈ C such that Z + ⊆ (X + ∪ Y + ) \ {e}

Z − ⊆ (X − ∪ Y − ) \ {e}. Note that if we remove the orientation property, these axioms are the circuit axioms of ordinary matroids, i.e. M := (E, C) is a matroid.

It can be shown that also a stronger form of axiom C3 can be derived.

Theorem 1 (Strong elimination, Theorem 3.2.5 in [17]) Let C be a collection of signed subsets of a set E satisfying (C0), (C1), (C2). Then (C3) is equivalent to C3’ for all X, Y ∈ C with X 6= −Y and e ∈ X + ∩ Y − and f ∈ (X + \ Y − ) ∪ (Y − \ X + ), there is a Z ∈ C such that Z + ⊆ (X + ∪ Y + ) \ {e},

Z − ⊆ (X − ∪ Y − ) \ {e}, and f ∈ Z.

12

2.1. WITHOUT INFORMATION ON METABOLITE POTENTIALS

Now we can use this definition to interpret metabolic networks as oriented matroids. Since we are only working with signs in oriented matroid theory, we have to translate the real flux vectors into sign vectors: Definition 7 (Sign function) For x ∈ Rn , define sign(x) := s ∈ {−, 0, +}n with:   − if xi < 0 si := 0 if xi = 0 for all i ∈ 1, . . . , n   + if xi > 0 We will show, given a metabolic network without irreversible reactions, that the steady state flows with minimal support form an oriented matroid. We need the condition that the steady state fluxes must have minimal support to satisfy property C2. To work with minimality of circuits in a mathematically formal way, we define for signed subsets X = (X + , X − ), Y = (Y + , Y − ) of some base set E: X ⊆ Y :⇔ X + ⊆ Y + ∧ X − ⊆ Y − . Obviously ⊆ is a partial ordering on the signed subsets of E. Let V be a collection of signed subsets of E. We can define the set of minimal signed subsets: M in(V) := {v ∈ V : @w ∈ V \ {v} s.t. w ⊆ v} In the interpretation of metabolic networks as oriented matroids, the collection V be the set of flux directions of steady state fluxes and C = M in(V \ ∅) will be the circuits: Proposition 1 (Metabolic networks are oriented matroids) Let N = (M, R, S) be a metabolic network, then (R, C) defines an oriented matroid, where C = M in (V) and V = {sign(J) : SJ = 0, J 6= 0}. Proof This result is not a special result related to metabolic networks. To be precise, generating an oriented matroid from a real valued matrix (like the stoichiometric matrix) is one of the most standard ways to obtain an oriented matroid. Because this method is so important, oriented matroids that can be generated this way are called realizable [141, 17]. We will give the proof anyway, because it helps to understand the oriented matroid axioms. In particular the proof is rather constructive, hence it will be helpful when implementing algorithms for metabolic networks that are stated in terms of oriented matroids. We have to check the circuit axioms: C0 Satisfied, since C ⊆ {sign(J) : SJ = 0, J 6= 0}. C1 Assume this property is violated, then there must exist a X ∈ C such that −X 6∈ C. Since X ∈ C, there exists a J 6= 0 with SJ = 0 and sign(J) = X. Since S(−J) = 0 and sign(−J) = −X, it follows that −X ∈ V. Hence, there must exist a Y ∈ C with Y ⊂ X. By the same argument it follows that −Y ∈ V, which is a contradiction to X ∈ C, since −Y ⊂ X. C2 Assume there exists X, Y ∈ C with X ⊆ Y and X 6= Y 6= −X. We know that there exist J X , J Y with SJ X = 0 = SJ Y and X = sign(J X ), Y = sign(J Y ). Now consider Z(λ) = J Y + λJ X . Obviously Z(0) = J Y and thus sign(Z(0)) = Y . For λ large, respectively small enough we also know that sign(Z(λ)) ⊇ X respectively 13

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

sign(Z(λ)) ⊇ −X. Since −Y 6= X 6= Y and X ⊆ Y it follows by continuity of Z that there exists a λ such that sign(Z(λ)) ⊆ Y and Z(λ)e = 0 for some e ∈ X. This is a contradiction to the minimality of Y . C3 Since X, Y are circuits, there exist J X , J Y such that X = sign(J X ), Y = sign(J Y ). Since e ∈ X + ∩ Y − , it follows that JeX > 0, JeY < 0. Hence, there exists a λ > 0 such that JeX + λJeY = 0. Obviously, Z := J X + λJ Y , is also a steady state flux with Z + ⊆ (X + ∪ Y + ) \ {e}, Z − ⊆ (X − ∪ Y − ) \ {e}. By definition of C, it follows that there exists a Z 0 ∈ C with Z 0 ⊆ Z, since X 6= −Y and thus Z 6= ∅. Hence, (R, C) is an oriented matroid.



The following two examples of oriented matroids will be important for us: ˙ S), let Definition 8 (Flux Mode Matroid) Given a metabolic network N = (M, R = I ∪E, (R, C) denote the oriented matroid obtained from (M, R, S) by Proposition 1. We will call (R, C) the flux mode matroid. 2 Note, that the circuits of the flux mode matroid are the minimal sign vectors of steady-state fluxes. Hence, we have a direct correspondence to elementary modes, which are discussed in Section 4 in more detail. ˙ S), Definition 9 (Internal Cycle Matroid) Given a metabolic network N = (M, R = I ∪E, let (I, C) denote the oriented matroid obtained from (M, I, SI ) by Proposition 1. We will call (I, C) the internal cycle matroid. 2 The internal cycle matroid consists of all the circuits that we do not want to have in thermodynamically feasible solutions. As we will see, it is intrinsically linked to the space of potential differences by oriented matroid duality. See Figure 2.3 for an example on circuits of the two matroids. c

c d e

a

f

b

f

b

C 1 = ({a, b, c, e, f }, ∅)

C 2 = ({a, c, d, e, f }, ∅)

c

c d e

a

d e

a

f

d e

a

b

f

b

C 3 = ({a, b, e, f }, {d})

C 4 = ({b, c}, {d})

Figure 2.3: C 1 , C 2 , C 3 , C 4 , −C 1 , −C 2 , −C 3 , −C 4 are the circuits of the flux mode matroid. C 4 , −C 4 are the only circuits of the internal cycle matroid. We have now defined the oriented matroids we want to work with by their circuits. Oriented Matroid theory on the other hand offers many more possibilities to define an oriented matroid. One such alternate notion are vectors. Intuitively, vectors are simply circuits where we relaxed the minimality condition. 14

2.1. WITHOUT INFORMATION ON METABOLITE POTENTIALS

Definition 10 (Definition 3.7.1 in [17]) Let (E, C) be an oriented matroid. We say that, Y ∈ {−, 0, +}E is a vector of the oriented matroid, if there exist circuits X 1 , . . . , X n with Y = X 1 ◦ X 2 ◦ . . . , ◦X n , where ( (X ◦ Y )e =

Xe Ye

if Xe 6= 0 else

The following proposition follows directly from this definition: Proposition 2 Let O be an oriented matroid with base set E, circuits C and vectors V. Let V ∈ V be a vector and e ∈ V . Then there exists a circuit X ∈ C with X ⊆ V .

2

It can be shown that we can also decompose null-space vectors into circuits (Lemma 6.7 in [141]), hence the set V that we used to define the circuits of metabolic networks are precisely the vectors of the oriented matroid. Proposition 3 Let O = (E, C) be an oriented matroid realized by a matrix S. Then the set of vectors of O is V = {sign(J) : SJ = 0, J 6= 0}.

2

In graph theory, the union of circuits is called a cycle; hence, we will often call the vectors of an oriented matroid also the cycles of the oriented matroid. As already mentioned, oriented matroids are interesting for metabolic networks, because it gives us a relationship between the internal cycle matroid and the space of potential differences. This relationship is characterized by oriented matroid duality. We defined realizable matroids by the null-space U of a matrix. The dual matroid will be the matroids that are realized by the orthogonal subspace U ⊥ , and this is also the motivation for the following abstract definition (see Proposition 3.7.12 in [17]): Definition 11 (Dual oriented matroid) Let O be an oriented matroid on elements E and circuits C). The set of vectors V ∗ of the dual oriented matroid O∗ is V := {Y ∈ C ∗ : X ⊥ Y ∀X ∈ C} where X ⊥ Y if • X ∩ Y = ∅ or • there exist e, f ∈ X ∩ Y such that Xe Ye = −Xf Yf . O∗ = (E, C ∗ ) is called the dual oriented matroid of O. The elements of C ∗ are called the cocircuits of O and V ∗ is called the set of covectors (cocycles). 2 Proposition 4 It holds that O∗∗ = O. Proof See Proposition 3.4.1. in the book by Bj¨orner et al. [17]. 15



CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

For realizable oriented matroids, the orthogonality relation ⊥ has the following intuitive meaning: If X ⊥ Y holds, then there exist vectors J X , J Y such that X = sign(J X ), Y = sign(J Y ) with (J X )T J Y = 0. In particular, it can be shown (see Proposition 6.8 in [141]) that for a given matrix S ∈ Rm×n C ∗ = M in ({sign(µS) : µ ∈ Rm } \ {∅})

is precisely the set of cocircuits of the oriented matroid realized by S in the way of Proposition 1. Using linear algebra, we can realize the dual matroid by the null space matrix of S, hence the dual matroid is also realizable. We conclude:
Corollary 1 (I, C ∗ ) with C ∗ = M in {sign(µSI ) : µ ∈ RM } \ {∅} is the dual matroid of the internal cycle matroid. 2

2.1.2

Internal Cycles

We have now nearly all the tools at the ready to prove the main theorem of this subsection. This theorem will be fundamental for all further work in this thesis, since it relates thermodynamic feasibility to the existence of internal cycles. We only need one more little but very useful and important lemma that tells us that if a flux J is thermodynamically feasible, we always can get a thermodynamically feasible flux with all potential differences unequal to zero. Lemma 1 (Potential Perturbation Lemma) Given a metabolic network N = (M, R = ˙ S), then a flux vector J ∈ RR is thermodynamically feasible if and only if there exists I ∪E, potentials µ ∈ RM such that µSi 6= 0 for all i ∈ I and J is thermodynamically feasible w.r.t. µ. Proof ⇐: It follows directly from the definition.

⇒: Assume J is thermodynamically feasible. Hence, there exists µ0 ∈ RM such that J is thermodynamically feasible w.r.t. µ0 . Since S does not contain any zero columns (Assumption 1), there exists for every i ∈ I a j(i) ∈ M such that ej(i) Si 6= 0 (where ej denotes the j.th unit vector). Since I, M are finite (Assumption 2), S is a continuous operator and hence, we can find δ ∈ lin{ej(i) : i ∈ I} such that δSi 6= 0 for all i ∈ I and |δSi | < |µ0 Si |. Hence J is thermodynamically feasible w.r.t. µ = µ0 + δ and µSi 6= 0 for all i ∈ I.  As a consequence of Corollary 1 we can now derive the main theorem: Theorem 2 (Equivalence of Thermodynamic Infeasibility) Given a metabolic network ˙ S) a flux vector J ∈ RR is thermodynamically feasible if and only if it holds N = (M, R = I ∪E, I for all v ∈ R with sign(v) ⊆ sign(JI ) and SI v = 0 that v = 0. Proof (Using oriented matroid theory) An idea of the proof was given by Beard et al. [9].

⇒: Assume J is thermodynamically feasible and there exists a circuit X ∈ C with X ⊆ sign(JI ). Because J is thermodynamically feasible, there exists a µ ∈ RM such that sign(µSI ) ⊇ sign(JI ).

Select e ∈ X. Since sign(µSI ) ⊇ X, there exists by Proposition 2 a cocircuit Y ∈ C ∗ with e ∈ Y and Y ⊆ sign(µSI ). Since Y ⊥ X by duality and Yf = Xf for all f ∈ Y ∩ X, it follows that X ∩ Y = ∅, which is a contradiction. Hence there cannot exist a v 6= 0 with SI v = 0 and sign(v) ⊆ sign(JI ), because else there would exist a circuit X ⊆ sign(v). 16

2.1. WITHOUT INFORMATION ON METABOLITE POTENTIALS

⇐: Assume there exists no circuit X ∈ C with X ⊆ sign(JI ). Define E = I\supp(JI ). By Assumption 2, E is finite, thus we can write E = {e1 , e2 , . . . , en }, where n = |E|. Define Z0 := sign(JI ). We will now define Zi inductively such that there exists no circuit X ∈ C with X ⊆ Zi . Obviously this property is satisfied by Z0 . Define + − Yi := (Zi−1 ∪ {ei }, Zi−1 ),

+ − Yi0 := (Zi−1 , Zi−1 ∪ {ei }).

Claim 2.1.1 For at least one of Yi , Yi0 it holds that there exists no circuit X ∈ C with X ⊆ Yi (or X ⊆ Yi0 respectively). Proof Assume there exist circuits X, X 0 ∈ C with X ⊆ Yi , X 0 ⊆ Yi0 . By construction ei ∈ X, X 0 . Since S does not contain any zero columns by Assumption 1, it follows that ei is not a loop; hence, X 6= −X 0 and we can apply axiom C3. Since ei is the only entry where X, X 0 have different signs, it follows by axiom C3 that there exists a circuit Z with e 6∈ Z and Z ⊆ X ∪ X 0 ⊆ Zi−1 . This is a contradiction to the condition that Zi−1 does not contain any circuits.  Define Zi to be Yi if there exists no circuit X ∈ C with X ⊆ Yi , else set Zi := Yi0 . It follows by construction that there exists no circuit X ∈ C with X ⊆ Zn . Since Zn also has no zero-entries, it follows that X ⊥ Zn for all X ∈ C. Hence, Zn is a covector by definition of duality and thus, there exists a µ ∈ RM with sign(µSI ) = Z ⊇ sign(JI ). We conclude that J is thermodynamically feasible.  Proof (Using LP-duality theory) We have the following chain of equivalences: ∃v ∈ RI : sign(v) ⊆ sign(JI ), SI v = 0, v 6= 0 ( ≥ 0 if Ji ≥ 0 I ⇔∃v ∈ R : 0 < sign(JI )v, SI v = 0, vi ∀i ∈ I ≤ 0 if Ji ≤ 0 ( ) ( ≥1 Ji ≥ 0 M ⇔∅ = µ ∈ R : µSi ∀i ∈ I ≤ −1 Ji ≤ 0

(2.1) (by LP duality theorem)

Since we do not have any additional constraints on µ, we can multiply µ by −1 and obtain that by Lemma 1 and scaling the last line is equivalent to thermodynamic infeasibility.  Internal cycles are sometimes also called Type III pathways due to the pathway classification of Schilling et al. [113]. They should not be confused with type II pathways, which are sometimes called futile cycles and are considered biologically important [74, 102]. The second proof of Theorem 2 already gives a method for checking whether a given steady state flux J is thermodynamically feasible. We simply have to solve the LP (2.1), i.e. maximize sign(JI )v and check if the optimal solution value is greater than 0. 17

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

This idea can also be used to compute a thermodynamically feasible flux from a given infeasible flux J, see [98, 85]: Input: A steady state flux J repeat I + := {i ∈ I : Ji > 0} I − := {i ∈ I : Ji < 0} L := arg max{1LI + − 1LI − : S(I + ∪I − ) L(I + ∪I − ) = 0, JI − ≤ LI − ≤ 0, JI + ≥ LI + ≥ 0} J := J − L until 1LI + − 1LI − = 0 return J Algorithm 1: Removing thermodynamically infeasible cycles Note that this algorithm works if there are no additional constraints on the flux J, like negative upper bounds or positive lower bounds. Such constraints may get violated during the update steps of the algorithm, but the sign of the flux values on the arcs will stay the same or the flux gets reduced to 0: Proposition 5 (Correctness) Given a steady state flux J as input, when Algorithm 1 terminates, it returns a thermodynamically feasible flux J ∗ with the property that: sign(J ∗ ) ⊆ sign(J). Proof Let J i denote the value of J at the end of the i.th iteration of the algorithm. Let Li denote the value of L computed in the i.th iteration. Hence we have by definition of the algorithm J0 = J J i = J i−1 − Li

J n = J ∗ if Ln = 0. Li always exists because the LP that computes L is always feasible since 0 is a feasible solution. Sine Ln = 0, it follows by Thm. 2 that J n is thermodynamically feasible. Since Li is feasible, it satisfies JIi−1 ≤ LiI − ≤ 0 and JIi−1 ≥ LiI + ≥ 0 and thus, sign(J i ) ⊆ sign(J i−1 ) for all i = 1, . . . , n. − + The proposition follows by induction.  Proposition 6 (Termination and running time) Algorithm 1 runs in polynomial time. Proof In each iteration one cycle is reduced by its maximal amount. Thus, the flux of one reaction (which wasn’t zero before) becomes zero. Cycles are only searched on those reactions that have nonzero flux. Thus, reactions that didn’t carry any flux will also not carry any flux after the iteration. Hence the number of reactions carrying flux strictly decreases with each iteration. Hence, the algorithm terminates after at most |I| iterations. Since we can solve LPs in polynomial time, the algorithm runs in polynomial time.  Since the LP that needs to be solved in each iteration only changes slightly, LP warm starting may give huge performance improvements. Nigam and Liang [85] propose a method for identifying violated cycles where they do not need to solve an LP in each iteration. Instead, they compute minimal bases of the null space matrix (these are circuits) and check, whether there exists such a basis vector that has the same signs as the flux vector. They state, that this is a necessary condition for thermodynamic infeasibility (although the lemma they cite only gives that it is a sufficient condition). Indeed, it is not a necessary condition, as can be seen in Figure 2.4. 18

2.2. WITH INFORMATION ON METABOLITE POTENTIALS

d c

a b

d e f

g

c

a

c

d e

b

e f

h

g

c

g

b h

Figure 2.4: The network is drawn in black, the circuits of the basis in blue. As indicated by the blue arrows, every circuit contains at least one negative entry, although there exists a circulation that is strictly positive on every arc. Nigam and Liang [84, 85] also try to find strongly thermodynamically feasible fluxes, but due to topological properties of the feasible flux space of strongly thermodynamically feasible fluxes (see chapter 3), they are doomed to fail. Hence, their algorithm returns a weakly thermodynamically feasibly flux, if it is not able to find a strongly thermodynamically feasible flux. Finally. we also observe that Algorithm 1 only updates fluxes of reactions that participate in internal cycles. This is a very nice property, since in practice often only the fluxes through exchange reactions are of interest. Observation 1 (Fluxinvariance of exchange reactions) Algorithm 1 keeps the flux through reactions that do not participate in internal cycles invariant. In particular, the flux through exchange reactions is kept invariant, i.e. JE∗ = JE for input flux J and output flux J ∗ .

2.2

With Information on Metabolite Potentials

With increasing advances in metabolomics, more and more systems biologists want to include metabolite concentration information into metabolic pathway analysis [11, 138, 100, 76, 77, 36, 41, 27]. As already explained in the introduction, we usually can compute the chemical potential µi of a metabolite i ∈ M from its concentration value ci by the following formula: µi = µ0i − RT ln(ci )

(1.3)

In practice R, T can be considered as known constants, but this is usually not the case with µ0i and the concentration ci . Although, µ0i values have only been measured for very few metabolites, they can be estimated for nearly all metabolites using the group contribution method developed by Mavrovouniotis [82], which was further improved by Jankowski et al. [67]. These estimates are usually surprisingly precise, however, to produce reliable results one has to work with error tolerances. Also, most methods that measure concentration values do not operate precisely. Hence, we will assume that for every metabolite we are given lower and upper bounds on the potential. Definition 12 (Thermodynamically feasible flux w.r.t. potential bounds) Given a me˙ S) and lower and upper bounds on metabolite potentials tabolic network N = (M, R = I ∪E, `, u ∈ RM , we call a flux vector J ∈ RR thermodynamically feasible with respect to `, u if there exists µ ∈ RM such that • ` ≤ µ ≤ u and 19

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

• J is thermodynamically feasible in N w.r.t. µ.

2

Some of the following results hold on an even more generalized form of the above definition, where the metabolite potentials can be restricted to an arbitrary set. Definition 13 (Thermodynamically feasible flux w.r.t. potential space Q) Given a me˙ S) and a potential space Q ⊆ RM , we call a flux vector tabolic network N = (M, R = I ∪E, R J ∈ R thermodynamically feasible with respect to Q if there exists µ ∈ Q such that J is thermodynamically feasible in N w.r.t. µ. 2 Inspired by the results on the case where nothing was known on metabolite potentials, we now want to derive similar results. Indeed, there exists a similar theorem to Theorem 2, which was already discovered by Mavrovouniotis in 1996 [83]. Although Mavrovouniotis was a bit sloppy regarding strict and weak inequalities and the statement of the theorem as well as his proof are not completely correct, he had the correct idea anyway. Since the proof of Mavrovouniotis does not work out properly (distinction between strict and weak inequalities is essential when applying Farkas Lemma), we will give a new proof. But this proof is similar in the sense, that we also use LP-duality theory. Elemental is the notion of minimally infeasible sets (MIS) called bottlenecks by Mavrovouniotis, which in this case can be represented as v ∈ {−, 0, +}I like the circuits of the internal cycle matroid. Indeed, if all lower bounds are −∞ and all upper bounds are ∞, the minimally infeasible sets will be precisely the circuits of the internal cycle matroid. Similar to the circuits, we also use the signed set notation v = (v + , v − ) with v + = {i : vi = +} and v − = {i : vi = −}. In Figure 2.5, you can see a network with an infeasible set that is not a circuit of the internal circuit matroid. e a

c

b 5

[0, 10]

d 5

Figure 2.5: Because the potential difference between the left and right metabolite is 5 − 5 = 0, no flux is possible through reactions b, c at the same time. Hence, ({b, c}, ∅) is an infeasible set and the only feasible pathways are ({a, b, e}, ∅) and ({c, d}, {e}). To capture the notion of infeasible sets more precisely, we introduce the notion of a subnetwork . Each vector v ∈ {−, 0, +}I represents a subnetwork of the original network in the following way: ˙ S) be a metabolic network. Definition 14 (Subnetwork) Let N = (M, R = I ∪E,

Given a vector v ∈ {−, 0, +}I we call the network N v = (M, I v := I ∩ v, S v ) the subnetwork of N w.r.t. v. The stoichiometric matrix S v of the subnetwork contains the same stoichiometries v of the forward reactions as S, i.e. S∗i = S∗i for all i ∈ v + . For the backward reactions the signs v of the stoichiometric coefficients are flipped, i.e. S∗i = −S∗i for all i ∈ v − . 2 If we are given a flux J and are interested in properties on the reactions that carry flux, we can simplify things by looking at the subnetwork of flux carrying reactions. This subnetwork is exactly the network N v , where v = sign(J). ˙ S) be a metabolic network. Let J be a flux in N . Proposition 7 Let N = (M, R = I ∪E, 20

2.2. WITH INFORMATION ON METABOLITE POTENTIALS

Let v = sign(J). Then J is a steady state flux, if and only if J v := (vJ)v is a steady state flux in N v . Additionally J v > 0. 2 Thus, it suffices to check whether it is thermodynamically feasible if all reactions in N proceed in forward direction. If it is not feasible, we have found an infeasible direction configuration that also cannot occur in the original network. This infeasible direction configuration will by captured by the notion of infeasible set (IS): ˙ S) be a metabolic network with Definition 15 (infeasible set) Let N = (M, R = I ∪E, M metabolite potential space Q ⊆ R . A vector v ∈ {−, 0, +}I is called an infeasible set, if v there exists no µ ∈ Q such that µS∗i < 0 for all i ∈ I v of the subnetwork N v . 2 These infeasible sets behave as one naturally would expect infeasible sets to behave. I.e. if you make the set larger, it still stays infeasible: Proposition 8 Let v ∈ {−, 0, +}I be an infeasible set in a metabolic network N . Let v ⊆ w ∈ {−, 0, +}I . Then w is also an infeasible set. Proof By inclusion follows that I v ⊆ I w and S v = SIwv . Hence, if there exists a µ ∈ Q such that µSiw < 0 for all i ∈ I w , the same µ would also disprove infeasibility of v.  Assume now, we are given a flux vector J. If J is thermodynamically feasible, all reactions where J has positive flux must have negative potential differences. Translated to the theory of infeasible sets, sign(J) must not be an infeasible set. As a consequence of Proposition 8 we now only have to check for each (inclusion) minimal set v whether v ⊆ sign(J). ˙ S) be a metabolic network with metabolite potential space Theorem 3 Let N = (M, R = I ∪E, Q ⊆ RM . Then there exists a collection of minimally infeasible sets C ⊆ {−, 0, +}I such that a flux vector J is thermodynamically feasible w.r.t to Q if and only if v 6⊆ sign(JI ) for every v ∈ C. Proof According to Definition 15, J is thermodynamically feasible if and only if sign(−JI ) is not an infeasible set. Since I is finite (Assumption 2) there exists by Proposition 8 a finite number of minimal infeasible sets C such that sign(JI ) is an infeasible set if and only if there exists a v ∈ C with v ⊆ sign(JI ).  In the case, where the potential bounds are −∞ and +∞, we know that the minimal infeasible sets are the circuits of the internal cycle matroid. Since oriented matroids are well studied, we know how to find those circuits and how to check, whether sign(J) of a flux J contains such a circuit. The following Theorem 4, which was originally stated by Mavrovouniotis [83] in a similar fashion in 1996, gives answer to the generalized question. The theorem as stated here, corresponds to statements (I) and (V) of Theorem 3 in the article by Mavrovouniotis [83]: Theorem 4 Let N = (M, R, S) be a metabolic network with lower and upper bounds on metabolite potentials `, u ∈ RM . Then the following are equivalent: 1. There exists a µ ∈ RM with ` ≤ µ ≤ u such that µSi < 0 ∀i ∈ R. P 2. For all v ∈ RR with i∈R vi = 1 it holds that there exists a µ ∈ RM with ` ≤ µ ≤ u such that µSv < 0. 21

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

Proof The first statement holds if and only if the following minimization problem has a solution less than 0: min{α : α − µSi ≥ 0 ∀i ∈ R, ` ≤ µi ≤ u ∀i ∈ M}

(2.2)

Since primal and dual LP of (2.2) are feasible it follows from the duality theorem that the solution value is equal to

= =

max{x` + yu : −Sv + x + y = 0, 1v = 1, v ≥ 0, x ≥ 0, y ≤ 0}

(2.3)

max{x` − yu : −Sv + x − y = 0, 1v = 1, v ≥ 0, x ≥ 0, y ≥ 0}

max{(x − y)` − y(u − `) : −Sv + x − y = 0, 1v = 1, v ≥ 0, x ≥ 0, y ≥ 0}

We will now do a component wise analysis of the solution vector. If `i = ui , then xi `i − yi ui = Si∗ v`i . If `i 6= ui we can observe for fixed v that (xi − yi )`i − yi (ui − `i ) = Si∗ v`i − yi (ui − `i ) will only be maximized if yi is minimized (since ui − `i > 0). Hence, we have ( 0 if Si∗ v ≥ 0 yi = −Si∗ v if Si∗ v < 0 ( Si∗ v if Si∗ v ≥ 0 ⇒ xi = 0 if Si∗ v < 0 Since Sv is a fixed vector for fixed v, the minimum of min{µSv : ` ≤ µ ≤ u} is obtained by the component wise minimums, i.e. µi = `i if Si∗ v > 0 and µi = ui if Si∗ v < 0. It follows that (for both cases `i = ui and `i < ui ) ⇒ xi `i − yi ui

= Si∗ v · (µmin )i , where µmin = arg min µSv : ` ≤ µ ≤ u µ

⇒ x` − yu = µmin Sv = min{µSv : ` ≤ µ ≤ u} We conclude that the following are equal: min{α : α − µS∗i ≥ 0 ∀i ∈ R, ` ≤ µi ≤ u ∀i ∈ M}

= max{min{µSv : ` ≤ µ ≤ u}, 1v = 1, v ≥ 0} Thus, the lemma holds.



In the theory of infeasible sets, we can restate Theorem 3 as follows: ˙ S) be a metabolic network with lower and upper bounds Corollary 2 Let N = (M, R = I ∪E, M on metabolite potentials `, u ∈ R .

A steady state flux J ∈ RR is thermodynamically feasible in N w.r.t. the bounds `, u if and only if for all v ∈ RI with sign(v) ⊆ sign(JI ) and v 6= 0 exists a µ ∈ RM with ` ≤ µ ≤ u such that µSI v < 0. 2 So the main result of Corollary 2 is, that if J is not thermodynamically feasible, we can find a witness that proves infeasibility. In the cases of circuits these are the internal circulations v, which satisfy SI v = 0. Now it is easy to see, that all internal circulations are also infeasible sets. 22

2.2. WITH INFORMATION ON METABOLITE POTENTIALS

Another special case is used by Nolan, Lee, Xu, Boghihian et al. [87, 137, 19]. They are only given potentials of boundary metabolites. In this case, all minimal infeasible sets are elementary modes where the exchange reactions have been removed. For general metabolite concentration bounds we observe: Observation 2 The proof of Theorem 4 also gives us an LP formulation that either proves thermodynamic feasibility of a flux J or returns a witness that proves infeasibility. We only need to solve the LP (2.3) on the subnetwork N sign(JI ) . If the solution value is less than zero we obtain feasibility. If the solution value is larger or equal to zero, we can use v, which is part of the solution vector, as the witness for infeasibility.

23

CHAPTER 2. THERMODYNAMICALLY FEASIBLE PATHWAYS

24

Chapter 3

Flux Space In Chapter 2 we discussed methods on determining thermodynamic feasibility of steady state fluxes. Here we will take a more general perspective and try to understand the space of all thermodynamically feasible steady state fluxes.

3.1

Structural Properties

To start of, we will discuss the space of steady state fluxes without thermodynamic constraints. There are several degrees of information detail that require different kinds of analyses. The most basic information that is assumed to be known is whether a reaction is reversible or not, i.e. whether it is possible to send flux backwards through the reaction. Using this information the so called flux cone can be defined. ˙ S) be a metabolic network. Let D denote Definition 16 (flux cone) Let N = (M, R = I ∪E, the set of irreversible reactions of N . Then F := {J ∈ RR : SJ = 0, Ji ≥ 0 ∀i ∈ D} is the flux cone of N with irreversible reactions D.

(3.1) 2

It follows straight from the definition that the flux cone F of a metabolic network D is a closed polyhedral cone. Often, also lower and upper bounds `, u on the flux values are given. Obviously, we can model reaction directions also using flux bounds. Hence, it follows for the space P of all feasible reactions with respect to the flux bounds that P := {J ∈ RR : SJ = 0, ` ≤ J ≤ u} is a closed polyhedron . If we do not allow `i = −∞ or ui = ∞ for some reactions i, then P will also be bounded and hence, be a polytope . In practice, usually, only very few bounds are known (except directions), so the only reactions that are constrained with lower and upper bounds are exchange reactions. This means that if N contains internal cycles; flux through these cycles is usually not bounded. This shows one major 25

CHAPTER 3. FLUX SPACE

reason, why thermodynamic constraints should be used in the analysis of metabolic networks. Unbounded, or very large flux will falsify results obtained by analyzing the steady state flux space significantly [100, 109]. As we know from Chapter 2, thermodynamic constraints prohibit flux through internal cycles. This important property in practice comes with drawbacks on the mathematical properties. While the space of steady state fluxes P was closed and convex, these properties are lost if we add strong thermodynamic constraints. Proposition 9 There exist metabolic networks, where the space Tstrong := {J ∈ P : ∃µ s.t. µSi Ji < 0 or µSi = 0 = Ji ∀i ∈ I} of strongly thermodynamically feasible fluxes is neither closed nor convex.    1 0 −1 1 , Proof Consider the network N = M = {a, b}, R = {1, 2, 3, 4}, S = 0 −1 1 −1 where I = {3, 4} and E = {1, 2}. The network is depicted in Figure 3.1 and a projection of the strongly thermodynamically feasible fluxes can be seen in Figure 3.2. J4 , ∆µ4

J1 = 1 = J2 J3 , ∆µ3 3 1

a

b

2

∆µ

4 Figure 3.1: Already this network has a nonFigure 3.2: The set of possible values for the convex feasible flux space potential differences ∆µ = µSI is shown by the red line. The projection of the thermodynamically feasible fluxes on to the (J3 , J4 )plane is marked in light green. The open interval of strongly thermodynamically feasible fluxes with J1 = 1 = J2 is marked in dark green. As can easily seen by the pictures, the thermodynamically feasible flux space is not convex, since J 1 = (3, 3, 2, −1)T is strongly thermodynamically feasible w.r.t. µ1 = (1, 0) and J 2 = (−3, −3, −1, 2)T is strongly thermodynamically feasible w.r.t. µ2 = (0, 1) but 21 (J 1 + J 2 ) = 1 T 2 (0, 0, 1, 1) contains a inner circulation and hence is not strongly thermodynamically feasible. This is illustrated by the light green area in Figure 3.2. To see that the strongly thermodynamically feasible flux space is not closed, consider the sequence
T {J n }n∈N with J n = 1, 1, 1 − n1 , − n1 . Every J n , n ∈ N is strongly thermodynamically feasible w.r.t. µ1 . But J ∗ = limn→∞ J n = (1, 1, 1, 0)T is not strongly thermodynamically feasible, since it holds for all µ ∈ RM that sign(∆µ3 ) = −sign(∆µ4 ). This is illustrated by the dark green line in Figure 3.2.  26

3.2. THERMODYNAMICALLY FEASIBLE FLUX THROUGH A GIVEN (INTERNAL) REACTION The relaxation to weakly thermodynamically feasible fluxes solves at least one problem: The space of weakly thermodynamically feasible fluxes is closed. This is a very important property for the optimization methods later on (c.f. Chapter 5), because the existence of optimal solutions is not clear if the feasible domain of the optimization problem is not closed. We will now give a proof where the set of possible metabolite potentials µ ∈ RM can be arbitrarily constrained. Lemma 2 Let Q ⊆ RM . The space T := {J ∈ P : ∃µ ∈ Q s.t. µS∗i Ji < 0 or Ji = 0 ∀i ∈ I} of (weakly) thermodynamically feasible fluxes is closed.

Proof Let (J n )n∈N be a sequence in T converging to J ∗ . Since J n → J ∗ , there exists an n ∈ N such that ∀i : Ji∗ > 0

Jin > 0

∀i : Ji∗ < 0.

Jin < 0

Hence sign(JI∗ ) ⊆ sign(JIn ) and thus J ∗ is also thermodynamically feasible by Proposition 8.



But in general the set of weakly thermodynamically feasible fluxes is not the closure of the strongly thermodynamically feasible fluxes. An easy example can be seen on Figure 3.2 if P = {J : J1 = 1, J4 ≥ 0}. Then it holds that Tstrong = ∅ but T = {(1, 1, 1, 0)T }.

The space of weakly thermodynamically feasible fluxes is also not convex, as can easily be seen from Figure 3.2, too. This is the property that turns thermodynamic constraints computationally hard and there have been attempts to simplify these constraints. For example, Yang et al. [140, 139] derived irreversibility properties of reactions from thermodynamic constraints and thus computed a convex set that contains all thermodynamically feasible fluxes.

3.2

Thermodynamically Feasible Flux through a Given (Internal) Reaction

Although there are many cases in practice where only the flux through exchange reactions is of interest, there are also many applications where we also want to have information or constraints on fluxes through internal reactions. One such situation is a method called flux variability analysis, where we want to know how much the flux through any given reaction varies under all possible optimal fluxes. We will discuss flux variability analysis in greater detail in Chapter 5. Other biological applications are concerned with the rate of ATP production and hence the flux through the ATPase reaction; see [10, 9]. Hence, we will now deal with the problem of finding a thermodynamically feasible flux that has positive flux on a given (internal) reaction. We will only consider internal reactions, since we can answer the problem easily for exchange reactions by Observation 1. We will also assume that all reactions are irreversible, since this is the case for many reactions in genome-scale models nowadays. If the network has reversible reactions, these reactions can be substituted by forward and backward reactions. Problem 1 (Thermoflux) Given: • Metabolic network N = M, R = I ∪ E, S ∈ RM ×R 27




CHAPTER 3. FLUX SPACE

• objective (internal) reaction r ∈ R Question: Does there exist a thermodynamically feasible flux J ≥ 0 with Jr > 0 ? If our network is a directed graph, Thermoflux is closely related to the two-vertex-disjointpath-problem, Disjoint-Path for short (see Figure 3.3 for examples). In Disjoint-Path, we want to find two vertex disjoint paths from s1 to t1 and from s2 to t2 respectively. DisjointPath was shown to be NP-hard by Fortune, Hopcroft, Wyllie in 1980 [46]. I will use this to show that even this small additional restriction (positive flux on an internal reaction) makes the thermodynamic flux problem NP-hard. Problem 2 (Disjoint-Path) Given: • Digraph G = (V, A)

• s1 , t1 , s2 , t2 ∈ V

Question: Does there exist vertex disjoint paths from s1 to t1 and s2 to t2 ?

Figure 3.3: Corresponding Instances of Thermoflux and Disjoint-Path Theorem 5 Thermoflux is NP-hard. Proof We will reduce Disjoint-Path to Thermoflux. Let G = (V, A) be the given digraph and s1 , t1 , s2 , t2 ∈ V the vertices we want to find disjoint paths for. Let I be the incidence matrix of G. We now construct a metabolic network N by adding an input reaction a, an output reaction o and a target reaction r. The reaction a will deliver flow to s1 . Reaction r will be the internal reaction we are interested in in Thermoflux and which will transport flow from t1 to s2 . Reaction o will then serve as the output flow reaction from t2 . Formally, we have N = {V, A ∪ {a, o, r}, S}, with SA = I, Sa = es1 , So = −et2 , Sr = es2 − et1 , where ei denotes the i.th unit-vector. The internal reactions of N are I = {A, r} and the external reactions are E = {a, o}. Claim 3.2.1 There exists disjoint paths from s1 to t1 and from s2 to t2 in G if and only if there exists a thermodynamically feasible flow through r in N . Proof ⇒: Let P1 , P2 ⊂ A be the two disjoint paths. Set Ji = 1 for all i ∈ P1 ∪ P2 ∪ {a, r, o} and Ji = 0 else. Since P1 , P2 are paths, we have IJA = −es1 + et1 − es2 + et2 . It follows that SJ = 0, hence J is a steady-state flux. J is also thermodynamically feasible, since P1 , P2 are vertex disjoint; hence, every vertex on the path P1 ∪ P2 ∪ {a, r, o} has degree 2. Thus, no internal cycles are possible. 28

3.3. SAMPLING

⇐: Let J be a thermodynamically feasible flux with Jr > 0. Since a, o are the only exchange reactions, both must carry flux as well. Define P := {e ∈ A : Je > 0}. Assume there exists a path in (V, P ) from s1 to t2 that does not contain r. Let jP := min Je : e ∈ P
the flux through this path. Since J is thermodynamically feasible, J˜ := J − jP eP ∪{a,o} is also a thermodynamically feasible flux with Jr > 0. Thus, we can assume w.l.o.g. that every path from s1 to t2 contains r. Since Jr > 0, there exists such a path. Let us denote this path by P . Since J is thermodynamically feasible, P is simple, i.e. it does not contain any cycles. It follows that P \ {r} is the union of two vertex disjoint paths, one going from s1 to t1 and the other going from s2 to t2 . 

We have now shown the reduction of Disjoint-Path to Thermoflux. Hence, by the NPhardness result of Fortune et al. [46] it follows that Thermoflux is NP-hard.  As a side remark, note that the construction used in the reduction can also be done in reverse. Hence, we can also reduce the thermodynamic flux problem on directed graphs to the disjoint path problem. Although the disjoint path problem is NP-hard in general, it becomes tractable under practical conditions like small directed tree width [68]. Hence, the extensions of concepts like directed tree width to metabolic networks may also lead to tractability result for thermodynamic flux and algorithms that are efficient in practice. Another interesting remark is concerning the problem definition of Problem 1. If we do not want to require that J ≥ 0, we cannot adapt the proof accordingly. A reduction of the undirected vertex disjoint path problem will not prove NP-hardness, since the undirected vertex disjoint path problem is solvable in polynomial time [119].

3.3

Sampling

Probably the most direct way to understand the space of feasible fluxes is to simply generate feasible points at random and check if this set of points contains interesting properties. Sampling methods were for example successfully used to analyze the impact of diabetes, ischemia and diet on the metabolic network of human mitochondria [127]. Also, for understanding the interaction between the human alveolar macrophage and M. tuberculosis, sampling methods proved useful [20]. More examples and a comprehensive overview can be found in the article by Schellenberger and Palsson [110]. The randomly sampled points are used to estimate correlations and variances of flux variables. A high correlation between two flux-variables means that flux through these two reactions is very likely to be coupled [97]. This is important information needed to design experiments that give as much new information as possible. For example, if flux variance analysis tells you that two reactions are known to be coupled, it is not likely (unless your model is incorrect) that you will obtain new insights if you disable one reaction and then measure the flux through the other reaction. Sarıyar et al. also use this method to reduce the high dimensional flux space to a lower dimension that still reflects most of the information to recognize modular patterns in the network [106]. For these analysis methods to work properly and reliable, it is necessary that the sampled points are distributed as uniformly as possible in the feasible flux space. Genome-scale networks, however, consist of thousands of reactions, thus the feasible space has high dimension and the curse 29

CHAPTER 3. FLUX SPACE

of dimensionality kicks in. This means that classical decomposition methods fail, because the feasible space has to be decomposed into an exponentially growing amount of parts. Also rejection methods that sample on larger domains and simply reject infeasible points are not practical in most cases, since the relative amount of feasible points that can be retained goes to zero quickly. Because of this, a modified Markov Chain Monte Carlo method, called artificial centering hit-and-run (ACHR) [69], is used for sampling the flux space [110]. The artificial centering hit-and-run method bases on the Markov Chain Monte Carlo method called “hit-and-run” which works as follows: we start at an interior initial point x0 . We then choose a random direction and walk a random distance along this direction. The new point is then used as the starting point for the next iteration. It can be shown that the sample distribution converges in total variation to the target distribution (which in this case is the uniform distribution) [69]. The convergence rate, however, depends on the distribution with which the walking direction is chosen. To get a close to optimal walking direction distribution, the artificial centering step is added which uses the previously computed iterates and thus the Markov chain property is lost and thus also the theoretical convergence result. While the general hit-and-run algorithm works for arbitrary sets (also non-convex ones), the artificial centering hit-and-run method is motivated by point-symmetric polytopes. Practical experiments, however, showed that the method works also well for arbitrary polytopes [69]. But, as we have already seen, in the case of thermodynamically feasible flows, the set of feasible points is not convex anymore. To deal with this problem, Cogne et al. [27] simply sampled the unconstrained flux space and then filtered out only the thermodynamically feasible flux distributions. They observed that only 4% of all sampled flux distributions were thermodynamically feasible. Schellenberger et al. [109] followed a similar method but postprocessed all sampled points to turn them into the closest thermodynamically feasible points. This was done by solving a mixed integer linear program (MILP) or mixed integer linear quadratic program (MIQP) for each sampled point. This method then leads to many points that lie on the boundary of the thermodynamically feasible domain (in fact, their sampling method did not produce any cycle-free points initially, see Figure 4b in [109]). It is questionable whether the obtained distribution is even close to the uniform distribution of the thermodynamically feasible space. At first sight it even looks as if the method may only be sampling points in a space of measure zero, but, since the thermodynamically feasible space may be of lower dimension and thus the MILP method or MIQP method can find points in the relative interior, this does not need to be the case. A better method is suggested by Price et al. [100]. Instead of sampling the space of thermodynamically unconstrained fluxes, they first perform thermodynamically constrained flux variability analysis (see chapter 5) to determine the flux bounds on each of the reactions. They observe that for the H. pylori network only 0.1% of all samples violate thermodynamic constraints. These constraints can simply be rejected and thus the thermodynamically feasible flux space can be uniformly sampled.

30

Chapter 4

Generating Vectors, Extreme Rays and Elementary Modes Elementary modes and generating vectors are one of the most commonly used tools to analyze metabolic networks, infer potential regulatory sites and find knock-out targets [115, 117, 61, 32, 92]. Basically, generating vectors, extreme rays and elementary modes are reformulations of the the flux cone that can give deeper insights to the functions of metabolic networks [96]. They all have their foundation in Minkowski-Weyl’s theorem for polyhedral cones. See [58, 22] for an introduction into convex polytope theory. Cones can be described in two different ways: By an outer description (as Definition 16 of the flux cone) or an inner description (as generating vectors, etc.). A comprehensive study of elementary modes can be found in the PhD-Thesis of Terzer [124]. We will now briefly mention the main properties that are needed to understand the effect of thermodynamic constraints on inner descriptions. Definition 17 (Outer description of a cone C) A cone C ⊆ Rn is called finitely constrained or polyhedral if there exists a m ∈ N and matrix A ∈ Rm×n such that C = {x : Ax ≤ 0} Definition 18 (Inner description of a cone C) A cone C ⊆ Rn is called finitely generated if there exists a k ∈ N and matrix B ∈ Rn×k such that C = {Bλ : λ ≥ 0}. The column vectors in B are called the generators of C. Observe that every point of a cone can be written as a conical combination of its generators. The theorem of Minkowski-Weyl implies that, since the flux cone is finitely constrained, there also exists a finite number of generators of the flux cone: Theorem 6 (Minkowski-Weyl’s Theorem for Polyhedral Cones) A cone C ⊂ Rn is finitely constrained if and only if C is finitely generated. 2 It would be nice if the minimal set of generating vectors of a cone C would be uniquely determined. But this is only the case if C is pointed . A cone is pointed if it does not contain a linear subspace. 31

CHAPTER 4. GENERATING VECTORS, EXTREME RAYS AND ELEMENTARY MODES

Every cone can be decomposed as the Minkowski sum of a pointed cone (the recession cone) and a linear subspace (lineal) [58]. It can easily be seen that there exist flux cones that are not pointed, for example this is the case if all reactions are reversible. In practice biologists want to associate biological functions to the generating vectors of the flux cone. Hence, it is not very convenient if the latter cannot be uniquely determined. Hence, Larhlimi and Bockmayr proposed to analyze the flux cone in terms of minimal metabolic behaviors, which are basically the generating vectors of the recession cone [78]. An alternative workaround are so called extreme pathways. In extreme pathways every internal reaction is made irreversible. This can be done by adding reversed reactions for reversible reactions [96]. ˙ S) be a metabolic network with Definition 19 (Extreme Pathway) Let N = (M, R = I ∪E, irreversible reactions D ⊆ R. The extreme pathway cone is defined by o n F (EP) := (J + , J − ) ∈ RR × RI\D : SJ + − SI\D J − = 0, J − ≥ 0, Ji+ ≥ 0 ∀i ∈ I ∪ D . Let F (EP) = L+C be a decomposition of the extreme pathway cone into the lineal L = lin(F (EP) ) and a pointed cone C. Let R denote the set of extreme rays of C and a set of generators of L. The projections of the elements of R onto the flux cone, i.e.  R∗ := prEP (J + , J − ) : (J + , J − ) ∈ R , with + + prEP (J + , J − ) = J, with JE∪D = JE∪D , JI\D = JI\D − J −,

are called extreme pathways of N .

2

It holds by construction of the extreme pathway cone that   F = prEP F (EP) Hence, each set of extreme pathways (which are not unique, if lin(F (EP ) ) 6= {0}) generate the flux cone. Note that in general the extreme pathway cone is not guaranteed to be pointed, since only internal reactions are made irreversible. Exchange reactions may still be reversible and hence if reversible exchange reactions are linearly dependent, it holds that lin(F (EP ) ) 6= {0} and thus F (EP ) is not pointed. But in practice it can be assumed that reversible exchange reactions are not linearly dependent. If also the exchange reaction reactions are split into irreversible reactions, one ends up with elementary modes (EM, also called elementary flux modes) [96]. Definition 20 (Elementary Mode) Let N = (M, R, S) be a metabolic network with irreversible reactions D. Call U := R \ D the set of reversible reactions. The elementary mode cone is defined by +  U − F (EM) := (J + , J − ) ∈ RR + × R+ SJ − SU J = 0 . Let R denote the set of extreme rays of F (EM) . The projections of the elements of R onto the flux cone, i.e.  + R∗ := prEP (J + , J − ) : (J + , J − ) ∈ R , prEP : (J + , J − ) 7→ J, with JE∪D = JD , JU = JU+ − J − are called elementary modes of N .

2

32

4.1. ENUMERATION OF ELEMENTARY MODES

U Since F (EM) ⊆ RR + × R+ , the elementary mode cone is pointed, hence the set of extreme rays of (EM) F generates the elementary mode cone. By the same argument as for the extreme pathways it follows that the elementary modes are generators of the flux cone. In particular, they are also extreme pathways. Note that, because of the projection step, not all EMs may be needed to generate the flux cone.

Hence, the set of elementary modes is uniquely defined for every metabolic network. Several enumeration algorithms exists for elementary modes which allows easy post processing. This makes them one of the most favorite tools of systems biologists.

4.1

Enumeration of Elementary Modes

An extensive study of enumeration algorithms can be found in the PhD-Thesis by Terzer [124] and in [48]. Enumerating elementary modes is not an easy task, since first of all, the number of elementary modes usually grows exponentially with the number of reactions [48]. For example, consider the network shown in Figure 4.1. Every reaction is irreversible; hence, there are no internal cycles and thus the number of elementary modes equals the number of paths through the network. For every node xi , i > 2 there are two possibilities for reaching xi : either from xi−1 or xi−2 . Hence, we get for the total number of paths fi to node xi that fi = fi−1 + fi−2 . Since f1 = f2 = 1, it follows that fn is the n.th Fibonacci number. the Hence,   number of pathways in this example growth exponentially in n with fn = Θ x2

x1

xn−2

x4

x3

√ 1+ 5 2

n

.

xn

xn−1

x5

Figure 4.1: The number of elementary modes growth exponentially with n . But, the enumeration problem is not only hard, because there may be exponentially many elementary modes. Assume an enumeration algorithm already enumerated a set S of elementary modes. A result by Acu˜ na et al. [2] states that it is NP-hard to decide whether there exists an elementary mode x 6∈ S. This makes the enumeration problem #P -complete.

In particular, Acu˜ na et al. [3] also showed that enumerating all elementary modes x with xi > 0 for a given reaction i ∈ R cannot be done in polynomial total time. This follows by reduction from the following problem: Decide whether a given weighted directed graph has a negative weighted cycle containing a specified arc. This problem was shown to be NP-hard by Khachiyan et al. [70].

4.2

Elementary Modes and Oriented Matroids

Elementary modes and the circuits of the flux mode matroid are closely related [91]. This relationship will be very useful for further results and will motivate why we will often write 33

CHAPTER 4. GENERATING VECTORS, EXTREME RAYS AND ELEMENTARY MODES

elementary modes in the form of signed sets; like we write circuits. Although the relationship between elementary modes and the circuits of the flux mode matroid is mentioned by Oliveira et al., I did not find a mathematically rigorous proof. Hence, we will now show that the set of elementary modes equals the set of circuits of the flux mode matroid that satisfy the sign constraints (required by irreversible reactions). Elementary modes are flux rays, but circuits are sign-vectors. So we will show in a first step, that we do not loose any information by just looking at the signs. This is mainly due to the minimality of elementary modes. Lemma 3 Let x, y be elementary modes of a metabolic network N = (M, R, S). If sign(x) = sign(y), it follows that x = y. Proof Since x, y are elementary modes, there exist extreme rays x ¯, y¯ of the elementary mode cone that project to x, y respectively. If we can show that supp(¯ x) = supp(¯ y ) implies x ¯ = y¯, we are done, because then x = y follows immediately. Thus, we can assume w.l.o.g. that all reactions N are irreversible. Then the elementary mode cone equals the flux cone and we have F (EM ) = {J ∈ RR : SJ = 0, J ≥ 0}. We can also conclude that supp(x) = supp(y). Define the hyperplane H := {J : 1J = 1}. Since 0 6∈ H and F (EM ) ∩ H 6= ∅, it follows that F (EM ) ∩ H is a polytope and that there exists a bijective mapping q of extreme rays of F (EM ) to vertices in F (EM ) ∩ H: for all extreme rays r of F (EM )

q(r) = r ∩ H

q −1 (v) = {λv : λ ∈ R+ }

for all vertices of F (EM ) ∩ H

Also, the support is invariant under the transformation q. It is a basic fact from LP-theory (see for example [114]) that each vertex is uniquely identified by a basis. Let Bx , By be the bases of q(x), q(y) respectively. It follows that supp(x) ⊆ Bx , supp(y) ⊆ By . If supp(x) = supp(y), it follows that Bx is also a basis of q(y) and thus q(x) = q(y). By bijectivity of q we conclude x = y.  It follows that elementary modes are rays with minimal sign-support. Corollary 3 Let x be an elementary mode of a metabolic network N = (M, R, S). It holds that there exists no flux vector y 6= 0 of N with sign(y) ⊂ sign(x). Proof Follows directly from the basis-representation of Lemma 3.



We can even prove a stronger result that states that all rays of the flux cone with minimal support are also elementary modes. For this, the following lemma comes in very handy: Lemma 4 Let J ∈ F be a flux vector. Then J can be written as a conical combination of generators of elementary modes ei , i = 1, . . . , k with sign(ei ) ⊆ sign(x) for all i = 1, . . . , k. Proof Since the elementary modes span the flux cone, there exists a vector x = (x+ , x− ) ∈ RR × RU in the elementary mode cone that projects to J. 34

4.2. ELEMENTARY MODES AND ORIENTED MATROIDS

Claim 4.2.1 There exists a vector y = (y + , y − ) with prEM (y) = J and satisfies yi+ = 0 ∨ yi− = 0 for all i ∈ U. − + − Proof Assume there exists an i ∈ U s.t. x+ i , xi > 0. Define ε = min{xi , xi }. By definition of the projection onto the flux cone, also y with ( ( x+ if j 6= i x− if j 6= i − + j j ∀j ∈ R, y = ∀j ∈ U yj = j + − xj − ε if j = i xj − ε if j = i

projects onto J. Since y does not violate any sign constraints, it follows that y ∈ F . Additionally y satisfies yi+ = 0 ∨ yi− = 0 and the property is still satisfied for all j ∈ U for which it was already satisfied by x. This process can be repeated for all i ∈ U for which the property is violated.  − Hence, let us assume w.l.o.g. that x satisfies x+ i = 0 ∨ xi = 0 for all i ∈ U. Since x is in the elementary mode cone, it can be written as a conical combination of generators of extreme rays ri , i = 1, . . . , k. Since F (EM ) is a subset of the positive orthant, it follows that supp(x) ⊇ supp(ri ) for all i = 1, . . . , k. Thus sign(J) ⊇ sign(ei ) for all i = 1, . . . , k, where ei is the projection of ri to the flux cone. 

Corollary 4 The set of nonzero flux vectors with minimal sign-support are the generators of elementary modes. Proof It follows directly from Corollary 3 that the set of elementary modes is a subset of the set of flux vectors with minimal support. Let x be a flux vector with minimal support. By Lemma 4 it follows that x can written as a conical combination of elementary modes ej , j = 1, . . . , k with sign(ej ) ⊆ sign(x). Since x 6= 0 and x has already minimal support, it follows that x must be the conical combination of generators of exactly one elementary mode. Hence, x generates an elementary mode.  Since it is also the case here that it suffices to look at the signs, it is interesting to compare the relationship of elementary modes to oriented matroids. Define E := {sign(v) : v is an elementary mode} as the set of sign vectors of elementary modes. First of all we have to observe that (R, E) is not an oriented matroid, since axiom (C2) can be violated by the sign constraints. On the other hand, consider the flux mode matroid (R, C) (Definition 8). This matroid is studied extensively by Oliveira et al. in [89, 90, 91, 6]. If the metabolic network N does not contain any irreversible reactions, it follows that E is the set of circuits of the flux mode matroid. In general, we only have E ⊆ C. Theorem 7 Let N = (M, R, S) be a metabolic network with irreversible reactions D. Let E be the set of elementary modes e. Let C be the set of circuits of the flux mode matroid. It holds that there exists a bijective mapping ∼ between E and {C ∈ C : C − ∩ D = ∅} with e ∼ sign(e).

35

CHAPTER 4. GENERATING VECTORS, EXTREME RAYS AND ELEMENTARY MODES

Proof By definition of the flux mode matroid, every circuit C ∈ C has minimal sign support and there exists a vector v such that Sv = 0 and sign(v) = C. If additionally C − ∩ D = 0, v also satisfies the flux directions and hence, is a feasible flux vector. Since C 6= 0 by circuit axiom (C0), it follows that v 6= 0; thus, by Corollary 4, it follows that v generates an elementary mode. Let e ∈ E, it follows from the definition of the flux mode matroid that sign(e) ∈ C. Since e satisfies reaction directionalities, it follows that C − ∩ D = ∅. 

Because of this result, we will often simply use the bijective relation ∼ to identify elementary modes with circuits of the flux mode matroid; i.e. we will not distinguish in notation between the ray of the elementary mode and the corresponding circuit. The relationship between oriented matroids and elementary modes has also another nice consequence: We can use the same algorithm that we used for enumerating elementary modes to enumerate circuits of realizable oriented matroids. In particular, we can use elementary mode enumeration to enumerate the circuits of the internal cycle matroid (see Definition 9). In these cases however, intelligent preprocessing of the network will speed up the computation of the internal circuits significantly [136]. The greatest performance enhancement is obtained by removing all the reactions that cannot participate in internal circuits. These reactions can easily be identified using flux variability analysis, which will be discussed in Section 5.

4.3

Thermodynamically Infeasible Elementary Modes

In the case, where we do not have any constraints on metabolite concentrations, the characterization of thermodynamically infeasible elementary modes is easy. As we have seen, every elementary mode is a circuit of the flux mode matroid (M, C). Let CI denote the circuits of the internal cycle matroid. By the definition of these matroids (see Definitions 8, 9), it follows that CI ⊆ C. Let E again denote the set of elementary modes. We can conclude that the thermodynamically feasible elementary modes are E \ CI . In particular every thermodynamically feasible flux vector is contained in the cone spanned by elementary modes in E \ CI .

So, how does this look like if the range of metabolite concentrations is constrained? From Proposition 8 we know that if a flux v is thermodynamically infeasible w.r.t. to a metabolite potential space Q and sign(v) ⊆ sign(w), it follows that flux w is also thermodynamically infeasible. This property is directly captured by the notion of feasibility classifier introduced by Terzer: Definition 21 (Feasibility classifier, see Definition 6.2. in [124]) A function f (v) : RR → {0, 1} is called feasibility classifier for flux vectors v ∈ RR if for any flux vector w ∈ RR with sign(v) ⊆ sign(w), infeasibility (f = 0) of v implies infeasibility of w, i.e. sign(v) ⊆ sign(w), f (v) = 0 ⇒ f (w) = 0. 2

Hence, there exists a feasibility classifier for thermodynamic feasibility and we can apply the following Theorem 8 to thermodynamically feasible fluxes: Theorem 8 (see Theorem 6.1 in [124]) Let f be a feasibility classifier. Let F be the set of all feasible elementary modes w.r.t. f . Then every feasible flux (w.r.t. f ) is contained in the cone spanned by the elementary modes in F . 36

4.3. THERMODYNAMICALLY INFEASIBLE ELEMENTARY MODES

Proof The previous work on properties of elementary modes plays out nicely now; it allows us to give a short an elegant proof, compared to the proof given by Terzer [124]. Let J be a feasible flux with respect to the feasibility classifier f . By Lemma 4, there exist elementary modes ej , j = 1, . . . , k with J ∈ conv{ej : j = 1, . . . , k} and sign(ej ) ⊆ sign(J) for all j = 1, . . . , k. By the definition of feasibility classifier, it follows for every j = 1, . . . , k that ej is feasible.  We conclude that the convex hull of all thermodynamically feasible fluxes is spanned by the thermodynamically feasible elementary modes. ˙ S) be a metabolic network with irreversible reactions D Corollary 5 Let N = (M, R = I ∪E, and metabolite potential space Q. Then it holds for the space T of thermodynamically feasible fluxes w.r.t. Q conv(T ) = conv{e : e ∈ F }. where F is the set of thermodynamically feasible elementary flux modes w.r.t. Q.

37

2

CHAPTER 4. GENERATING VECTORS, EXTREME RAYS AND ELEMENTARY MODES

38

Chapter 5

Optimization Based Approaches Survival of the fittest, one of the fundamental principles in evolutionary theory is also the motivation behind a large branch of metabolic analysis methods. Given the set of possible reactions that may occur in a metabolic network, we assume that the species has evolved its regulatory control in such way that the metabolism optimizes its natural target: growth. For a cell to grow, certain bio-molecules like amino-acids, lipids, ATP need to be produced in a specific ratio. This ratio can be determined by analyzing the biomass decomposition of the cell. This biomass decomposition can be analyzed on different levels of detail and is also still focus of ongoing research [37]. From the mathematical perspective we are satisfied to know that such a biomass decomposition exists, can be formulated as an additional artificial reaction and that biologists are interested in maximizing biomass production [132, 99, 125, 34]. This analysis method is called flux balance analysis (FBA). But Biomass is not the only possible objective function. In the area of bioengineering, we do not want to maximize biomass, but by-products, like, for example, ethanol. Although FBA is very precise in predicting growth and yield rates [35], one optimal solution is very unreliable in predicting the flux value of by-products [71]. For such applications flux variability analysis (FVA) is a useful tool, which can be used to predict the range of possible by-product production rates under maximal biomass production. Hence, FVA serves as a compromise to elementary mode analysis. Practically the number of elementary modes (EM) grows exponentially with the number of reactions in a network [73], so the method becomes infeasible for genome scale networks. Usually, we are also not interested in the raw data of all elementary modes but only in specific properties. For example, we are interested in the maximal or minimal flux through a given reaction of the network. This question is answered by FVA [81, 16]. For maximization to make sense, we must not operate on the flux cone (see Definition 16) as we did for elementary modes. Every time we could have flux through the objective reaction, the flux would be unbounded. Hence we need to add bounds on the fluxes (see figure 5.1). There are two slightly different approaches possible: Growth maximization and yield maximization. If we have data on maximal flux rates (e.g. nutrient uptake rates), we can directly add these as upper bounds on the flux. This is usually referred to as growth maximization. If we do not have data on flux rates, we can still ask the following question: ”How much biomass can the cell produce per mole nutrient?” Mathematically, we also restrict the uptake rate of a nutrient, but we silently ignore flux bounds on all internal metabolites that may actually prove to be bottlenecks for maximal growth. Hence, the most efficient pathway (the one that maximizes 39

CHAPTER 5. OPTIMIZATION BASED APPROACHES

biomass r max biomass production

nutrient unutrient Figure 5.1: Since the nutrient uptake rate is bounded, the flux through reaction r is bounded and in particular also the production of biomass. yield) may use bottleneck reactions whereas a cell that wants to grow as fast as possible will not use pathways through bottleneck reactions and thus will not be yield-optimal [116, 126]. After all, we obtain for a metabolic network N = (M, R, S) with flux bounds `J , uJ a formulation of the form max cJ s.t. SJ = 0

(5.1)

J

J ≥`

J ≤ uJ

This is a nice LP that can be solved quickly even for genome scale metabolic networks.

5.1

Thermodynamically Infeasible Optimal Solutions

As we already observed in previous chapters not all of the fluxes in the feasible domain of (5.1) are thermodynamically feasible. In some cases the difference is not very huge, but in others, the results are simply completely unrealistic and unusable in practice. Lets start with the nice case. We want to maximize flux over exchange reactions (typically biomass production) and we do not have flux forcing internal reactions. A flux forcing reaction is a reaction that by the flux bounds is forced to carry flux. Definition 22 (flux forcing reaction) Given lower and upper flux bounds `J , uJ a reaction r is called flux forcing w.r.t. `J , uJ , if `Jr > 0 or uJr < 0. 2 Assume we are now only interested in the flux value on the objective reactions and do not have any information on metabolite potentials. In this case the following theorem tells us that we can safely ignore that the computed flux may be thermodynamically infeasible, because there exists another optimal flux that is feasible. ˙ S) be a metabolic network with flux bounds `J , uJ . Assume Theorem 9 Let N = (M, R = I ∪E, further, that no internal reaction i ∈ I is flux forcing. Let c ∈ RR be a cost function with cI = 0. Let J be an optimizer of (5.1). 40

5.2. GENERAL DIFFICULTY

Then exists an optimal steady state flow J ∗ through the metabolic network N with JE∗ = J E and `J ≤ J ∗ ≤ uJ that is thermodynamically feasible. Proof By Proposition 6, we can use Algorithm 1 to compute a flux J ∗ from J. By Propositions 5 it follows that J ∗ is thermodynamically feasible. By Observation 1 it follows that JE∗ = J E . Hence, J ∗ is an optimal flow, since it has the same optimal value as J. 

Price et al. [100] analyzed the effect of thermodynamic constraints on the genome-scale network of Helicobacter pylori. As expected by Theorem 9, they observed that thermodynamic constraints do not have much effect on reactions that are not involved in interior cycles. On the other hand they observed drastic changes for reactions contained in interior cycles. So, what can happen if one of the assumptions for the nice case is violated (e.g. if we optimize on a reaction that is contained in an interior cycle)? The result is most impressive, if we want to optimize over an internal reaction r ∈ I that participates in an internal circuit C, as we might want to do in FVA. An additional realistic assumption is that we do not have any flux bounds on the reactions in C. In this case, we can send unbounded amount of flux through the reactions in C. Thus, the LP (5.1) is unbounded and we obtain an entirely useless result. In the case of internal flux forcing reactions, the LP (5.1) may have a feasible solution (for example a circulation), but there may not exist any thermodynamically feasible flux that satisfies the flux bounds `J , uJ . If we are given information on metabolite concentrations, the minimal infeasible sets of reaction directions can be very wild. Hence, it will be very difficult to prove theorems like Theorem 9.

5.2

General Difficulty

Incorporating thermodynamic constraints into elementary flux mode analysis can easily be done using filtering techniques [124]. But this method has the serious disadvantage that it usually requires to create an exponential amount of data. For genome-scale networks this is a very unpractical property. In optimization based approaches, we can constrain our feasible domain to thermodynamically feasible fluxes right away. Without going into details on how to formulate these constraints, we can already state something on the difficulty of the problem. In Theorem 5 we already showed that it is NP hard to decide whether there exist a thermodynamically feasible flux through a given internal reaction. This is a very negative result regarding flux optimization methods, since it means that we will not even be able to find any polynomial time approximation algorithm for computing the maximum flow through internal reactions: Corollary 6 There exists no polynomial time approximation algorithm (of any quality) that solves max cJ : SJ = 0, `J ≤ J ≤ uJ , J is thermodynamically feasible in N ˙ S) with flux bounds `J , uJ ∈ RR and objective for a given metabolic network N = (M, R = I ∪E, R function c ∈ R . 2 Luckily, practical metabolic networks have nice structural properties that make the problem solvable (in a reasonable amount of time) even for genome-scale networks. We will now discuss 41

CHAPTER 5. OPTIMIZATION BASED APPROACHES

all currently known methods for formulating strong and weak thermodynamic constraints. In ˙ S) will be the metabolic network for which we will build the following, N = (M, R = I ∪E, mathematical programming formulations. `J , uJ ∈ RR will denote lower and upper bounds on the fluxes and `µ , uµ ∈ RM will denote lower, respectively upper bounds on the potentials. We will only analyze linear optimization of flux, hence c : RR → R will be a linear function and denote the cost function.

5.3

Nonlinear Programming

Nonlinear formulations are closest to the physical theory the models were build out of. Hence it is natural to solve these models directly. There also exists a rich theory on solving nonlinear programs and several implementations like IPOPT [15], Knitro [25], Lancelot [28] and many others [131, 44, 51, 120]. The most direct formulation follows directly from the definition of strongly thermodynamically feasible flux: max cJ s.t. SJ = 0 µS∗i Ji < 0 or µS∗i = 0 = Ji J

J

` ≤J ≤u

∀i ∈ I

(5.2)

`µ ≤ µ ≤ uµ

This formulation, however, is problematic in several ways: • The projection of the feasible domain onto the J variables is not closed (see Proposition 9). Hence, the existence of an optimal solution cannot be guaranteed even if the space of feasible flux values is bounded. • Equation 5.2 is quadratic, hence a linear programming solver is insufficient. By Proposition 9, it also follows that the program is not convex. Since Equation 5.2 is the only non-linear constraint, this means that the non-convexity is captured implicitly be the quadratic term. Global optimization of nonlinear programs, however, is very difficult [56, 13]. • Equation 5.2 contains an “or”. Such constraints can only be modeled by general constraint solvers or by mixed integer programs using boolean auxiliary variables. Observe that we would significantly restrict the solution space if we would force Ji 6= 0 for all i ∈ I. If the network for example contains isolated cycles or dead ends, the program would become infeasible immediately. • The constraint in (5.2) is strict. Although, this causes that the feasible domain is not open, solvers that use barrier method may well be able to cope with this property. On the other hand, the barrier method only finds local optima. Most of the criticism (except the part about closedness) also applies to the weak thermodynamic constraint, i.e. where Equation 5.2 is relaxed to µS∗i Ji < 0 or Ji = 0 ∀i ∈ I 42

5.3. NONLINEAR PROGRAMMING

The most promising approach in nonlinear formulations is followed by Heuett and Qian [63] and Fleming et al. [42]. It is motivated by the fact that physically every reaction actually carries forward and backward flux simultaneously. Hence, the flux vector is split into a forward and a backward part J = J + − J − . From kinetic rate laws the following equation can be derived [12]:  + Ji ∆µi = −RT ln Ji− where ∆µ = µS are the potential differences. The nonlinear program can then be formulated as follows (where k = RT is a constant): min cJ s.t. SJ = 0  ∆µi = −k ln

Ji+ Ji−

 ∀i ∈ I

∆µ = µS J = J+ − J−

`J ≤ J ≤ uJ `µ ≤ µ ≤ uµ

Although this formulation looks very different to (5.2), the two formulations are equivalent. ˙ S) be a metabolic network. Let k > 0 be a constant. A Theorem 10 Let N = (M, R = I ∪E, flux vector J ∈ RR is thermodynamically feasible w.r.t. metabolite potentials µ ∈ RM if and only if there exists J + , J − ≥ 0 such that J = J+ − J−  + Ji ∀i ∈ I µSi = −k ln Ji−

(5.3) (5.4)

Proof ⇐: Let J + , J − be given such that (5.3) and (5.4) are satisfied. Observe that it suffices to check strong thermodynamic feasibility for each reaction separately. Hence, for every reaction i ∈ I there are 3 cases to analyze J > 0 ⇔ Ji+ > Ji− : It follows that

Ji+ > 1 ⇒ ln Ji−



Ji+ Ji−

 > 0 ⇒ µS∗i < 0

and thus the sign constraint is satisfied, since J = J + − J − > 0. J < 0 ⇔ Ji+ < Ji− : It follows that  + Ji+ Ji < 1 ⇒ ln < 0 ⇒ µS∗i > 0 − Ji Ji− and thus the sign constraint is satisfied, since J = J + − J − < 0. J = 0 ⇔ Ji+ = Ji− : It follows that  + Ji Ji+ = 1 ⇒ ln = 0 ⇒ µS∗i = 0 Ji− Ji− and thus the sign constraint is satisfied, since J = J + − J − = 0. 43

CHAPTER 5. OPTIMIZATION BASED APPROACHES

⇒: Since we have to satisfy Equation 5.4, we have to choose J + , J − in a clever way. But again, we can do the analysis separately for every i ∈ I:

Ji = 0: Choose Ji+ = Ji− . Any positive value does the job. Ji < 0: Set Ji+ = α, Ji− =  α − Ji . Clearly Ji+ , Ji− ≥ 0 and Ji = Ji+ −  Ji− for any positive α. − − J J For α ↓ 0 we get k ln Ji+ → ∞. For α → ∞ we get k ln Ji+ → 0. By continuity i  − i J of the logarithm there exists an α such that ∆µi = k ln Ji+ . i

− + − Ji > 0: Set Ji+ = α + Ji , J α. Clearly Ji+ , Ji− ≥ 0 and Ji = J i = i −J i for any positive α. Ji− Ji− For α ↓ 0 we get k ln J + → −∞. For α → ∞ we get k ln J + → 0. By continuity i  −i J of the logarithm there exists an α such that ∆µi = k ln Ji+ .  i

This way Heuett and Qian got rid of the ugly “or-formulation” and were able to apply sequential quadratic programming [63] to find local optima. Of course, they do not magically get rid of the open-domain problem completely. Consider the network (shown here again in Figure 5.2), which we already used as an example in Proposition 9. Assume that reaction 3 shall be maximized, while flux through reaction 1 is set to 1. As we already saw in Proposition 9, the strongly thermodynamically feasible domain is relatively open and hence, the maximum does not exist. Let (J j )j∈N be a sequence of strongly thermodynamically feasible fluxes converging to the supremum (the weakly thermodynamically feasible flux (1, 1, 1, 0)T ). Since also reaction j − 4 must satisfy the thermodynamic constraint 5.4, it follows that (J j )+ 4 , (J )4 → 0 for j → ∞. Hence, the solution will converge to a weakly thermodynamically feasible flux without numerical explosions of J + , J − . 3 1

a

b

2

4 Figure 5.2: There exists no strongly thermodynamically feasible maximal flux through reaction 3 We showed that to every optimal steady state flux there exists a thermodynamically feasible flux that is also optimal if the cost function does only involve exchange reactions and we do not have any flux forcing internal reactions (Theorem 9). Fleming et al. [43] show a very similar result for strongly thermodynamically feasible fluxes using the formulation (5.4). ˙ S) with lower and upper bounds Theorem 11 (Theorem 1 in [43]) Let N = (M, R = I ∪E, on the exchange reactions `J , uJ ∈ RE . Let c ∈ RE be a cost function on the exchange reactions. Let J¯ ∈ RR be an optimal steady state flux in N , i.e. J¯ ∈ arg min{cJ : SJ = 0, `J ≤ JE ≤ uJ }. Then there exists a strongly thermodynamically feasible flux J ∈ RR with JE = J¯E .

Proof In the proof they formulate a strictly convex minimization problem on the internal flux variables JI+ , JI− :

+ 

 min JI , log(JI+ ) + d + 1I + JI− , log(JI− ) + d + 1I + − JI ,JI >0

s.t. SI (JI+ − JI− ) = b, 44

5.4. MIXED INTEGER LINEAR PROGRAMMING

where b = −SE J¯E , and d ∈ RI can be chosen arbitrarily.

Using KKT-theory they show that the optimal solution is strongly thermodynamically feasible. From the Lagrangian multipliers at the optimal solution the corresponding metabolite potentials can be computed. 

The theorem has the nice side result that we can efficiently compute thermodynamically feasible fluxes, since we only need to solve a convex minimization problem. The parameter d can be used to compute different thermodynamically feasible solutions; indeed, Fleming et al. also prove that every strongly thermodynamically feasible flux (with given exchange fluxes) can be obtained as a result of the minimization problem by appropriate choice of d (Theorem 2 in [43]).

5.4

Mixed Integer Linear Programming

An entirely different approach to the nonlinear formulation is followed by [9, 109, 27, 64, 62]. The weak thermodynamic constraint can be formulated as the following two implications on the flux J and the potential differences ∆µ = µSI : Ji > 0 ⇒ ∆µi < 0∀i ∈ I

Ji < 0 ⇒ ∆µi > 0∀i ∈ I

(5.5) (5.6)

Note, that if Ji = 0 then ∆µi can be anything. Hence, if and only if a steady state flux J satisfies (5.5) and (5.6) then J is weakly thermodynamically feasible.

5.4.1

Mixed Integer Linear Programming (MILP) formulations

To get an MILP-formulation, these implications are modeled using artificial boolean variables. There are several possibilities on doing this: Using 1 boolean variable per reaction For each reaction i ∈ I, we introduce an artificial variable ai . If Ji > 0, ai = 1. If Ji < 0, ai = 0. Using the lower and upper bounds `J , uJ on the flux, we can formulate this relationship as follows: `Ji (1 − ai ) ≤ Ji ≤ uJi ai As we can see, if ai = 1, this constraint takes the form 0 ≤ Ji ≤ uJi . If ai = 0, the constraint takes the form `Ji ≤ Ji ≤ 0. Hence, we have to set ai = 1 if Ji > 0 and set ai = 0 if Ji < 0. If Ji = 0 both ai = 1 and ai = 0 are feasible. Since we only can model weak inequalities in mixed integer linear programs, we cannot exactly reproduce the strict inequalities on ∆µ. Instead, we introduce an artificial ε > 0 that forces ∆µ away from zero. Since we also do not have bounds on the potential differences, we introduce a large artificial bound M that should not be restricting in practice: −M ai + ε ≤ ∆µi ≤ M (1 − ai ) − ε, 45

CHAPTER 5. OPTIMIZATION BASED APPROACHES

We can see that if ai = 1, it follows that −M + ε ≤ ∆µi ≤ −ε, hence ∆µ < 0. If ai = 0, we get ε ≤ ∆µi ≤ M − ε forcing ∆µi to be positive. We conclude that this formulation satisfies

Ji > 0 ⇒ ai = 1 ⇒ ∆µi ≤ −ε ⇒ ∆µi < 0∀i ∈ I Ji < 0 ⇒ ai = 0 ⇒ ∆µi ≥ ε ⇒ ∆µi > 0∀i ∈ I.

Hence, every feasible solution of this MILP formulation is weakly thermodynamically feasible. On the other hand it is not that clear whether the converse is also true, since we introduced the artificial variables ε, M . In particular, thermodynamically feasible solutions allow ∆µi = 0 if Ji = 0 which is not allowed in the MILP formulation. Thus, given any pair (J, µ) of a flux J that is thermodynamically feasible w.r.t. µ, we cannot always choose ε, M sufficiently small respectively large enough to make (J, µ) also feasible in the MILP formulation. Since we are only interested in objective functions on the flux in this chapter, only the existence of a µ ∈ RM is sufficient for thermodynamic feasibility (see Definition 4). Hence, given a thermodynamically feasible flux J, we only need to find a µ ∈ RM such that there exist ε, M such that (J, µ) is a feasible solution in the MILP. By the perturbation lemma (Lemma 1) this is indeed possible. Theorem 12 (Equivalence of MILP formulation on one variable) There exist ε, M > 0 such that the following to maximization problems are equivalent, i.e. αThermo = αMILP , where αThermo is the optimal solution of the thermodynamically constrained problem αThermo := max cJ : SJ = 0, `J ≤ J ≤ uJ , J is thermodynamically feasible in N and αMILP is the optimal solution of the MILP formulation αMILP := max cJ s.t. SJ = 0 ∆µ = µS `J ≤ J ≤ uJ

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + ε ≤ ∆µi ≤ M (1 − ai ) − ε ∀i ∈ I a ∈ {0, 1}I

Proof We already observed that if (J, µ) is a feasible solution of the MILP, J is also thermodynamically feasible. Hence, we have αThermo ≤ αMILP for any choice of ε, M .

For the reverse direction, remember that only the signs of J determine whether J is thermodynamically feasible, i.e. whether there exists a µ such that J is thermodynamically feasible w.r.t. µ. Since there is only a finite number of sign configurations of J, a finite collection Ω ⊂ RM of chemical potentials exists with the following property: For every thermodynamically feasible J there exists a µ ∈ Ω such that J is thermodynamically feasible w.r.t. µ.

By the perturbation lemma (Lemma 1), we can find for each µ ∈ Ω an equivalent µ0 such that ∆µ0i 6= 0 for all i ∈ I. Let Ω0 denote the finite collection containing one equivalent µ0 for each µ ∈ Ω. Since Ω0 is finite, we can define ε :=

min

µ∈Ω0 ,i∈I

|µSi |,

M := 46

max |µSi |.

µ∈Ω0 ,i∈I

5.4. MIXED INTEGER LINEAR PROGRAMMING

Now we can find for every thermodynamically feasible flux J a µ ∈ Ω0 such that J is thermodynamically feasible w.r.t. µ and by definition of ε, M it follows that (J, µ) is a feasible solution of the MILP formulation. Hence αThermo ≥ αMILP .  Using 2 boolean variables per reaction Alternatively, we can also formulate a MILP model, where we use two boolean variables ai , bi per reaction i ∈ I. This has the advantage that we do not need to use the perturbation lemma to show equivalence. Each of the two variables now only indicates one sign of Ji . If ai = 1, it means that Ji ≥ 0 is possible. If bi = 1, it means Ji ≤ 0 is possible. Since Ji cannot be positive and negative at the same time, only one of ai , bi can be set to one, which is enforced by the constraint ai + bi ≤ 1

`Ji bi ≤ Ji ≤ uJi ai

Hence, the main difference to the previous formulation is that ai = 0 = bi is also a valid configuration which does not allow any flux through reaction i. This, however, also allows ∆µi = 0 by formulating the sign constraint on ∆µi as follows: −M ai + εbi ≤ ∆µi ≤ M bi − εai The argumentation is the same as above. If ai = 1, it follows that ∆µi ≤ −ε; if bi = 1, it follows that ∆µi ≥ ε; and if ai = 0 = bi , it follows that ∆µi = 0.

Hence, we can conclude again that for every feasible point (J, µ) of this formulation, J is thermodynamically feasible w.r.t. µ. Since this formulation is weaker than the previous formulation using only one variable per reaction, it follows that this formulation is also equivalent to the original thermodynamic flux problem. Corollary 7 (Equivalence of MILP formulation on two variables) There exist ε, M > 0 such that the following two maximization problems are equivalent, i.e. αThermo = αMILP , where αThermo is the optimal solution of the thermodynamically constrained problem αThermo := max cJ : SJ = 0, `J ≤ J ≤ uJ , J is thermodynamically feasible in N and αMILP is the optimal solution of the MILP formulation αMILP := max cJ s.t. SJ = 0 ∆µ = µS `J ≤ J ≤ uJ

`Ji bi ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + εbi ≤ ∆µi ≤ M bi − εai ∀i ∈ I ai + bi ≤ 1 ∀i ∈ R a, b ∈ {0, 1}I

Note that we can also prove the corollary directly as we did prove Theorem 12 with the simplification, that we do not have to perturb µ such that µSi 6= 0 for all i ∈ I. 47

CHAPTER 5. OPTIMIZATION BASED APPROACHES

5.4.2

Quality of the LP Relaxation

Although solving MILPs is NP-hard in general, even large MILPs can be solved quickly by current solvers like Gurobi, CPLEX and SCIP, just to name a few. Solving MILPs has been an active research area for many years and several solving strategies have been developed. The book by Wolsey [135] gives a good overview on these methods. One of the main workhorses of mixed integer programming is the branch and bound algorithm . It is a very universal approach which cannot only be applied to MILP but also to more general constraint programming problems. For branch and bound to work, you need a method with which you can compute upper and lower bounds on the solution value of problems. In the case of maximization problems the upper bounds are computed by relaxations of the problem that are easier to solve. In MILP this bound is usually computed by the LP relaxation. The LP relaxation simply relaxes the integrality constraints, for example instead of a ∈ {0, 1} we solve the problem for a ∈ [0, 1]. This bound is often also called the dual bound , since every dual feasible solution of the LP-relaxation gives an upper bound. For the lower bounds already computed feasible solutions are used. The best solution that was computed is also called the incumbent solution. In MILP such solutions are either obtained if the optimal solution of the LP relaxation is integer, or by so called primal heuristics that simply guess some (primal) solutions. If for a problem the lower bound equals the upper bound, we know that we already solved the problem to optimality and we are done. If the lower bound is smaller than the upper bound, the problem is branched into several subproblems. One of these subproblems then must contain the optimal solution of the original problem. Branch and Bound is usually very effective, because not all of the subproblems need to be solved. If the upper bound of a subproblem is lower than the current value of the incumbent solution, no better solution can be obtained by the subproblem and it can be cut off. Input: Problem P incumbent = −∞ (no incumbent found yet) Q := {P } while Q 6= ∅ do choose p ∈ Q b := upper bound on p computed by a relaxation if b > incumbent then s := best solution of p computed by a primal heuristic (or −∞ if no solution found) if s > incumbent then incumbent := s end divide feasible domain of p into subproblems S Q := (Q \ {p}) ∪ S end end return incumbent Algorithm 2: General branch and bound routine for constraint solvers

48

5.4. MIXED INTEGER LINEAR PROGRAMMING

1 x, y, z ∈ {0, 1} incumbent = −∞ dualbound = 10

2

3 x = 0, y, z ∈ {0, 1} incumbent = 0 dualbound = 0

x = 1, y, z ∈ {0, 1} incumbent = 0 dualbound = 8

Problem 2 solved to optimality 4

6 x = 1, z = 0, y ∈ {0, 1} incumbent = 0 dualbound = 6

x = 1, z = 1, y ∈ {0, 1} incumbent = 6 dualbound = 4 Cut off, since dualbound ≤ incumbent Problem 1 solved to optimality

5 x = 1, y = 1, z = 0 incumbent = 6 dualbound = 6 Problem 4 solved to optimality Figure 5.3: Example run of a branch and bound algorithm on decision variables x, y, z. The problems are run according to their numbering. The incumbent at each node is the best solution computed by this and all previous nodes. The subproblem with x = 1, y = 0, z = 0 does not need to be solved, since Problem 4 is already solved to optimality after Problem 5.

We see, the quality of the LP relaxation is very important for the performance of MILP solvers. Let us now analyze the LP-relaxation of the MILP formulations. We will only investigate the formulation using one decision variable for each reaction, since the other formulation behaves analogously.

αLP := max cJ s.t. SJ = 0 ∆µ = µS `J ≤ J ≤ uJ

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + ε ≤ ∆µi ≤ M (1 − ai ) − ε ∀i ∈ I a ∈ [0, 1]I

49

CHAPTER 5. OPTIMIZATION BASED APPROACHES

By choosing ai =

1 2

for every i ∈ I we obtain αLP ≥ max cJ

s.t. SJ = 0 ∆µ = µS `J ≤ J ≤ uJ

`Ji uJ ≤ Ji ≤ i ∀i ∈ I 2 2 M M − + ε ≤ ∆µi ≤ − ε ∀i ∈ I 2 2 Since M is a large constant and `Ji , uJi are usually also very large (small in the case of `Ji ) because no data on maximal fluxes on internal reactions is known, the relaxation simply corresponds to FBA that is not thermodynamically constrained. Hence, if we optimize over an internal reaction r and if the network contains a circuit C that contains r, we can have a large circulation through r. This large circulation then leads to a very bad LP-relaxation value and hence the dual bound doesn’t say anything useful about the thermodynamically constrained problem. This result is actually not that surprising. In Corollary 6 we already observed that the problem is not polynomial time approximable. Hence, no relaxator that runs in polynomial time can give any information on the thermodynamically feasible flux for each possible network.

5.4.3

Cutting Planes for MILP formulations

A second workhorse of MILP solvers are cutting plane algorithms. Given a fractional (infeasible) solution, a cutting plane is generated that cuts of the solution but that does not cut off any of the feasible integer solutions. Algorithms that compute such cutting planes are also called separation algorithms. It follows by the ellipsoid method that if we can separate every infeasible solution in polynomial time that we also obtain a polynomial time algorithm for solving the optimization problem [57]. It follows that for NP-hard problems it is very hard (if not impossible) to find such a polynomial separation algorithm that can cut off any infeasible solution. Hence, one tries to find cutting plane algorithms that can cut off some of the infeasible solutions. A very successful strategy is to identify known substructures in a general problem and to derive cutting planes for these substructures. Flow Cover Cuts Such a family of cutting planes for mixed 0-1 programs like our thermodynamically constrained flow problems are flow cover cuts. Flow cover cuts are inspired by the following model: We consider a vertex of a graph with ingoing arcs N1 and outgoing arcs N2 . Each arc i ∈ N1 ∪ N2 has a maximal capacity ci for the flow which is only available if the arc is turned on, i.e. if the corresponding boolean variable ai = 1, else the capacity is 0, i.e. it holds for the flux x: xi ≤ ci ai ∀i ∈ N1 ∪ N2 50

5.4. MIXED INTEGER LINEAR PROGRAMMING

In addition we have the constraint that the inflow into the vertex must be less than the outflow, see Figure 5.4. Hence, in total the feasible set has the form  X := (x, a) ∈ Rn+ × {0, 1}n : x(N1 ) ≤ b + x(N2 ), xi ≤ ci ai ∀i ∈ N1 ∪ N2 , where b ∈ R is an unconditional flux. b N1

0 ≤ xi ≤ ci ai

0 ≤ xi ≤ ci ai

N2

Figure 5.4: Mixed 0-1 flow model of flow cover Proposition 10 (Flow cover inequality) Let C1 ⊆ N1 , C2 ⊆ N2 , L2 ⊆ N2 \ C2 and λ > 0. Then the flow cover inequality X x(C1 ) + max(cj − λ, 0)(1 − aj ) ≤ b + c(C2 ) + λa(L2 ) + x(N2 \ (C2 ∪ L2 )) j∈C1

is a valid inequality for X. Proof See for example the book by Wolsey [135].



So, where can flow cover cuts be applied in thermodynamically constrained flow problems? Of course, we can directly apply them to the flux variables; this is precisely the scenario they were designed for. But, we can also apply them to the potential differences. Let N be the null space matrix of the internal stoichiometries SI . Then we have ∆µ = µSI ∃µ ∈ RM

⇔ ∆µN = 0

Since the columns of N are by definition internal cycles, this leads us to the following interpretation. Let K be an internal cycle. We can find flow cover inequalities for X := {(x, a) ∈ RI × {0, 1}I : x ({i : Ki > 0}) ≤ x ({i : Ki < 0}) , xi ≤ |Ki |Mi ai }, which is simply a scaled instance (with xi = −∆µi |Ki |) of ˜ := {(∆µ, a) ∈ RI × {0, 1}I : ∆µK ≥ 0, ∆µi ≥ −Mi ai } X ⊇ {(∆µ, a) ∈ RI × {0, 1}I : ∆µK = 0, ∆µi ≥ −Mi ai }

⊇ {(∆µ, a) ∈ RI × {0, 1}I : ∆µK = 0, −M ai + ε ≤ ∆µi ≤ Mi (1 − ai ) − ε} =: T

(5.7)

˜ Since X ˜ is a Hence, we can translate the flow cover cuts for X to valid inequalities C for X. relaxation of an original circuit constraint (5.7) from the MILP formulation, the inequalities C are also valid for the original MILP formulation. 51

CHAPTER 5. OPTIMIZATION BASED APPROACHES

Observe that MILP solvers can only find these flow cover cuts easily if the the MILP contains the equation ∆µN = 0 If this is not the case, i.e. if we use the following (affinely transformed) formulation to constrain ∆µ by ∆µ = µS∗I we have no (in)equality that contains more than one variable of ∆µ. Hence, flow cover inequalities cannot be applied directly in this case. In Table 5.2 of Section 5.6 you can see that the way how ∆µ is constrained leads to significant changes in the running time. It is plausible that flow cover cuts may be responsible.

Combinatorial Benders’ Cuts and the Method of Beard et al. Note that the flow cover inequality also contains ∆µ next to the boolean variable a. But in the thermodynamically constrained flow problem as we defined it, we are actually not interested in ∆µ. We are only interested in feasible values of a. In such cases combinatorial Bender’s cuts turn out to be very useful [26]. As with the ordinary Bender’s reformulation the idea is to project out difficult variables of the MILP formulation; i.e instead of solving min cx : (x, y) ∈ P, x integer for some polytope P , we can also solve min cx : x ∈ projx (P ), x integer. In contrast to ordinary Benders’ cuts, we can use integrality of x to additionally strengthen the formulation by replacing projx (P ) by a polytope Q ⊆ projx (P ) that still contains all the integer points of projx (P ). In the optimal case, we can even find an integral Q, i.e. a polytope Q that only has integral vertices. We can apply combinatorial Benders’ cuts as follows to the thermodynamically constrained problem: For ease of notation define sign : {0, 1} → {−, +} with sign(0) = − and sign(1) = +. Hence, sign gives us the directions the 0/1-entries of a encode. If we apply sign to a vector, consider this as component wise application. Define the set  T := a ∈ {0, 1}I : sign(a) = −sign(µSI ) ∃µ ∈ RM Observe that sign(a) = −sign(µSI ) if and only if −M ai +ε ≤ µSI ≤ M (1−ai )−ε for sufficiently large M > 0 and small ε > 0. 52

5.4. MIXED INTEGER LINEAR PROGRAMMING

Using this, the thermodynamically constrained flux problem can be written as: max cJ s.t. SJ = 0 `J ≤ J ≤ uJ

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I a ∈ conv(T )

a ∈ {0, 1}I

To solve this problem efficiently, we need an outer description of conv(T ), or at least a separation algorithm that cuts of integer points that are not in T . The usual way to obtain this is by application of duality theory or Farkas Lemma to the inner description of T . In this case however, we already did all the necessary work in chapter 2. Remember that a simply encodes the signs of the flux vector, i.e. we have sign(JI ) ⊆ sign(a). Thus, we can describe the infeasible a using internal cycles. This is also the observation of Beard et al. [9]. The theorem that gave the relationship to the circuits was Theorem 2. We saw two different proofs of that theorem. The proof using oriented matroids is similar to the proof by Beard et al., while the direct proof, which uses duality theory, is the proof you would come up with if you analyze the problem from a combinatorial Benders’ cuts perspective. ˙ S) be a metabolic network and Corollary 8 Let N = (M, R = I ∪E,  T := a ∈ {0, 1}I : sign(a) = −sign(µSI ) ∃µ ∈ RM . Let a ∈ {0, 1}I . Then holds that a 6∈ T if and only if there exists a v ∈ RI with SI v = 0, v 6= 0 and sign(v) ⊆ sign(a). Proof If sign(v) ⊆ sign(a) it follows directly from Theorem 2 that there exist no µ ∈ RM such that sign(a) = sign(µSI ), hence a 6∈ T .

If a 6∈ T , it follows that no flux vector J with sign(J) = sign(a) is thermodynamically feasible, hence there exists by Theorem 2 a v ∈ RI with SI v = 0, v 6= 0 and sign(v) ⊆ sign(a).  As a side result we also get valid inequalities for T . If v is an internal cycle, then it follows that X

(1 − ai ) +

i:vi 0

ai ≤ #{i : vi 6= 0} − 1

(5.8)

is a valid inequality for T . By Corollary 8 we can find for every integral a 6∈ T such an internal cycle that witnesses infeasibility of a and hence, we find a separating inequality. We can rewrite (5.8) and obtain using (C = sign(v)) X

(1 − ai ) +

i:vi 0

X

ai ≤ #{i : vi 6= 0} − 1 ai ≤ |C + | − 1

i:Ci =+ −

⇔ a(C ) − a(C ) ≤ |C + | − 1 53

(5.9)

CHAPTER 5. OPTIMIZATION BASED APPROACHES

Observe that if C is not a circuit, there exists a circuit C˜ ⊂ C. The inequality for C˜ is stronger than the inequality for C, because a(C˜ + ) − a(C˜ − ) ≤ |C˜ + | − 1 ⇒ a(C˜ + ) + |C + \ C˜ + | − a(C˜ − ) ≤ |C + | − 1

⇒ a(C + ) − a(C˜ − ) ≤ |C + | − 1

⇒ a(C + ) − a(C − ) ≤ |C + | − 1

Thus, we conclude that it is sufficient to check the circuit inequalities (5.9):  T = {0, 1}I ∩ a ∈ [0, 1]I : a(C + ) − a(C − ) ≤ |C| − 1 for every circuit C ∈ C . So, we have a nice formulation for T using inequalities. But how good is this formulation? It would be ideal, if conv(T ) would be completely described by the bounds on a and the circuit inequalities, i.e.  conv(T ) = a ∈ [0, 1]I : a(Ci+ ) − a(C − ) ≤ |C + | − 1 for every C ∈ C Unfortunately, this does not hold in general. For example, consider the metabolic network (depicted in Figure 5.5) N6 = ({a, b, c}, {1, 2, 3, 4, 5, 6}, S) with 1 2 3  a 0 0 1 S = b 0 1 0 c 1 0 0

3

4 −1 −1 0

5 −1 0 −1

6  0 −1 . −1

5 a

1 c

4

6 b

2 Figure 5.5: The polytope conv(T ) corresponding to this metabolic network is not completely described by the circuit inequalities (5.9). The boxes represent the reactions, the circles represent the metabolites. The circuits of the network can easily be read of the figure. Obviously, there are no loops and no parallel reactions. For each choice of using exactly one of the reactions 4, 5, 6, we obtain one circuit (and its opposite). If we use two of reactions 4, 5, 6, we are also only left with one possible circuit (and its opposite) that we can build. If we use all three of reactions 4, 5, 6 we 54

5.4. MIXED INTEGER LINEAR PROGRAMMING

have three possibilities by using one of 1, 2 or 3. These are all possibilities because there are no circuits containing 5 or 6 reactions by circuit axiom (C2) of oriented matroids. Hence we have the following circuits (up to opposites): (+, +, 0, 0, 0, +)

(+, 0, +, 0, +, 0)

(0, +, +, +, 0, 0)

(+, 0, −, −, 0, +)

(+, −, 0, −, +, 0)

(0, +, −, 0, −, +)

(+, 0, 0, −, +, +)

(0, +, 0, +, −, +)

(0, 0, +, +, +, −)

Now consider the polytope described by the circuit inequalities:  P := a ∈ [0, 1]I : a(C + ) − a(C − ) ≤ |C + | − 1 for every C ∈ C Claim 5.4.1 The point a = ( 21 , 12 , 12 , 1, 1, 1) is a vertex of P . Proof We need to check that a is a basic feasible solution. a is feasible: Obviously a ∈ [0, 1], so we only have to check the circuit inequalities. Since a is symmetric, it suffices to check one of each type of circuit inequality: 1 ≤ a1 + a2 + a6 ≤ 2 −1 ≤ a1 − a3 − a4 + a6 ≤ 1 0 ≤ a1 − a4 + a5 + a6 ≤ 2

1 1 + +1=2 2 2 1 1 satisfied since a1 − a3 − a4 + a6 = − − 1 + 1 = 0 2 2 3 1 satisfied since a1 − a4 + a5 + a6 = − 1 + 1 + 1 = 2 2

satisfied since a1 + a2 + a6 =

a is basic: If a is a vertex, there must exist a basis for a, i.e. we must find 6 inequalities that are satisfied by a with equality and have linearly independent left hand sides. Such inequalities are: a4 ≤ 1

a5 ≤ 1 a6 ≤ 1

a1 + a2 + a6 ≤ 2

a1 + a3 + a5 ≤ 2

a2 + a3 + a4 ≤ 2

These inequalities are linearly independent, because using Laplace expansion we can easily compute   0 0 0 1 0 0  0 0 0 0 1 0      1 1 0  0 0 0 0 0 1   =2   1 0 1 det = −2 ⇒ det    1 1 0 0 0 1  0 1 1  1 0 1 0 1 0  0 1 1 1 0 0 55

CHAPTER 5. OPTIMIZATION BASED APPROACHES

Thus, a is a basic feasible solution and hence, a is a vertex of P .



We conclude that P is not an integer polytope and hence, cannot be equal to the convex hull of the feasible integer points T . Of course, we can try to derive additional inequalities to strengthen the formulation. The polytope is closely related to the set covering polytope; hence, we can use cuts that were developed for the set covering polytope (see for example, [7, 108]) also for T . Although we do not have a complete description of conv(T ) (because there is no complete description of the set covering polytope), the reformulation still is a strengthening in contrast to the original formulation. In the original formulation the direction indicator variables (a) only need to be decreased slightly (depending on the artificial ε, M ) to allow a potential difference in the opposite direction (than originally encoded by a). Before we can apply combinatorial Bender’s cuts in practice, there is one more problem: The network may contain exponentially many internal circuits. Thus, enumerating all internal circuits is not a method that scales well. Thus, we also need to have a method with which we can quickly check feasibility and separate infeasible solutions. As already observed in the direct proof of Theorem 2, given an integer vector a we can easily check if a ∈ T , by solving the following LP: ( ≥ 0 if ai = 1 ∀i ∈ I max sign(JI )v : SI v = 0, vi I v∈[−1,1] ≤ 0 if ai = 0 If the optimal value is less or equal to zero, we have found no internal cycle and a is feasible. If the optimal value is greater 0 the solution vector v is an internal cycle and hence we also found an inequality that we can separate. Since v need not be a circuit, we may want to improve the cut by searching for circuits contained in sign(v). This can be done by a local search routine that heuristically blocks some reactions and resolves the LP or by modifying the objective function [39, 52]. Combinatorial Benders’ Cuts with Information on Metabolite Potentials In Chapter 2 we also analyzed thermodynamic feasibility of steady state fluxes if the metabolite potentials are constrained to intervals. In Theorem 3 we proved that there also exist minimal infeasible sets of reaction directions such that a ∈ {0, 1}I feasible ⇔ C 6⊆ sign(a) for each minimal infeasible set C ∈ {−, 0, +}I . The cuts (5.9) we derived in the circuit case can also be applied to infeasible sets C ∈ {−, 0, +}I by a(C + ) − a(C − ) ≤ |C + | − 1 and it follows directly for x ∈ {0, 1}I that  x feasible ⇔ x ∈ a ∈ [0, 1]I : a(C + ) − a(C − ) ≤ |C + | − 1 for all minimal infeasible sets C . The only question that we need to solve to use this in practice is again a separation routine. In the case of circuits this was very easy, because the infeasible sets for circuits were easy to find. But also in the more general setting, it suffices to solve an LP and we also already did all the work in Chapter 2. The LP is part of the proof of Theorem 4, which established the result that 56

5.4. MIXED INTEGER LINEAR PROGRAMMING

we can show that a set C ∈ {−, 0, +}I is an infeasible set by finding a witness and solving an LP . In the case of an all positive direction vector we observed (Observation 2) that we only need to solve the following LP to determine feasibility: max{x`µ + yuµ : −Sv + x + y = 0, v 1 = 1, v ≥ 0, x ≥ 0, y ≤ 0} If the optimal solution is less than zero, the network is feasible. Else, we find a solution vector v that proves infeasibility. It turned out that supp(v) is an infeasible set. Thus, we only have to translate this result also for direction vectors with mixed signs. We can do this by simply reversing all the reactions i ∈ I with ai = 0. Reversing reaction i means that we flip the signs of all the entries of row i in S. This again is equivalent to changing the sign of vi . Hence, given a reaction vector a ∈ {0, 1} we can test its feasibility by solving ( ) ( ≥ 0 if a = 1 i max x`µ + yuµ : −Sv + x + y = 0, v 1 = 1, x ≥ 0, y ≤ 0, vi ∀i ∈ I . ≤ 0 if ai = 0 Since the only change we made was on v which is not part of the objective, it follows that if the optimal solution is less than zero, then the network is feasible. Since sign(v) ⊆ sign(a) it follows that sign(v) is an infeasible set if the objective is greater or equal to zero. Thus we are also able to solve the separation problem in the case where we want to use metabolic potential information: Theorem 13 (Separation) The separation problem for integer valued flux directions a ∈ {0, 1}I ˙ S) with given bounds `µ , uµ on the metabolite poin a metabolic network N = (M, R = I ∪E, tentials can be solved in polynomial time. The same holds true if no bounds on the metabolite potentials are given. Proof The result follows directly from the observations above.



Analogously to the circuits case, we can strengthen the cut using a local search routine that minimizes the infeasible set. Combinatorial Bender’s Cuts for Flux Variables We can also introduce combinatorial Benders’ cuts on the flux variables. Fischetti et al. call such cuts optimality cuts [39]. Since we optimize on the flux variables, this is not as easy as introducing cuts on the potentials. But, we can consider the following feasibility problem, which checks if a solution with value greater than some constant jopt exists: ∃a, J, µ s.t. cJ ≥ jopt SJ = 0

`Ji (1 − ai ) ≤ Ji ≤ uJi ai

−1000ai + (1 − ai ) ≤ µSi ≤ −ai + 1000(1 − ai ) a ∈ {0, 1}I

By adaptively increasing jopt we can use this feasibility problem to compute the optimal flux. 57

CHAPTER 5. OPTIMIZATION BASED APPROACHES

Since µ, J are only linked by the flux directions a, we can project out µ and J separately and introduce combinatorial Benders’ cuts for each. See the article by Codato et Fischetti [26] on how to compute combinatorial Benders’ cuts for ordinary systems described by weak inequalities. Let T denote the polytope described by the combinatorial Benders’ cuts for µ. Let Q denote the polytope described by the combinatorial Benders’ cuts for J respectively. Thus, we can rewrite the feasibility problem as: ∃a ∈ {0, 1}I : a ∈ T, a ∈ Q Hence, we only have to solve a series of integer programs (IPs) to solve the problem. On the other hand, experiments in practice showed, that the minimal infeasible sets that were generated to cut of infeasible solutions of Q consisted of many reactions, hence these cuts were not very strong and many cuts had to be generated. We will now argue, why this method is likely not to perform well. For the following, consider the case where jopt is very small, so that feasibility only depends on the sign of the bounds `J , uJ and cr = 1 for one reaction r ∈ R and cp = 0 for all p 6= r. Assume further that `J ≤ 0, uJ ≥ 0. By reversing reactions, we can assume w.l.o.g. that uJi > 0 for all i ∈ R. Let D denote the set of irreversible reactions, i.e. D = {i ∈ R : `Ji = 0}. We will now establish a connection between minimal cut sets , introduced by Klamt et al. [73], and the combinatorial benders cuts on the flux variables.

Definition 23 (Minimal Cut Set (MCS)) A cut set is a set of reactions X such that if all the reactions in X are blocked (removed from the network), then there is no positive flux possible through the target reaction r. A minimal cut set is a cut set that is inclusion-wise minimal. (The restriction to positive flow through r is sometimes omitted, but we need it and Ballerstein et al. [8] need it, too.) 2 Ballerstein et al. [8] showed, that each minimal cut set for a target reaction r can be identified to a irreducible inconsistent subsystem (IIS) of the following linear system: SJ = 0

α

J =0

β

Ji ≥ 0 ∀i ∈ D

δ

cJ ≥ jopt



Bellerstein et al. proved that if X is a minimal cut set for r, then there exists a irreducible infeasible subset that contains exactly the rows  and βX (βX denotes the rows setting flux on reaction X to zero) and some rows of α, δ. We will now slightly modify this system to SJ = 0

α

J ≥0

β

Ji ≥ 0 ∀i ∈ D

δ

J ≤0

γ

cJ ≥ jopt

(5.10)



This means, we are now only blocking one direction of the reactions, this leads to the novel notion of oriented cut sets. 58

5.4. MIXED INTEGER LINEAR PROGRAMMING

Definition 24 (Oriented Minimal Cut Set (OMCS)) An oriented cut set is a signed set of reactions (X + , X − ) such that there exists no flux vector J with JX + ≤ 0, JX − ≥ 0 and Jr > 0.

An oriented minimal cut set is a oriented cut set that is inclusion-wise minimal (according to the subset relation we defined for signed sets in Section 2.1.1). 2

For the following we will use the result that elementary modes can be identified with circuits of the flux mode matroid , c.f. Theorem 7, and use circuit-notation as well as oriented matroid theory. This leads to the following, equivalent alternative definition of oriented cut set: Definition 25 An oriented cut set is a signed set of reactions (X + , X − ) such that there exists no elementary mode E with E + ∩ X + = ∅ = E − ∩ X − and r ∈ E + . 2 Lemma 5 Each OMCS (X + , X − ) can be identified with an IIS of (5.10) that contains exactly the rows , βX − , γX + and some rows of α, δ. Proof The proof works analogously to the proof of Lemma 1 in [8]. Since (X + , X − ) is an OMCS, it follows that the following system is inconsistent: SJ = 0 Ji ≥ 0 ∀i ∈ X − Ji ≤ 0 ∀i ∈ X + Ji ≥ 0 ∀i ∈ D

cJ ≥ jopt

Since this system may not be irreducible, it remains to show that there exists an irreducible subsystem of it that also contains the rows , βX − , γX + . Assume, this would not be the case. Obviously, every system that does not contain  is consistent (by J = 0). So, there exist Y − ⊆ X − , Y + ⊆ X + such that one of the inclusions is strict and there exists an infeasible subsystem that only uses rows , βY − , γY + in addition to some of the rows of α, δ. It follows that (Y + , Y − ) is an OMCS, a contradiction to the minimality of (X + , X − ).  The system (5.10) can also be used to describe the set of feasible integer points in Q (i.e. the signs of steady state fluxes through r) a ∈ Q ⇔ {βi:ai =0 , γi:ai =1 , , α, δ} is not an infeasible subsystem of (5.10). Combinatorial Benders’ cuts cut off all the configurations of a such that {βi:ai =0 , γi:ai =1 , , α, δ} becomes infeasible. Thus, we have a one-to-one correspondence of combinatorial benders cuts for Q and oriented cut sets. Hence, also the minimality translates and we can bijectively associate every irredundant combinatorial benders cut to an OMCS. The following result shows that the number of OMCSs is greater than the number of MCSs. This is bad, because already the number of MCS explodes corresponding to the network size in practical examples [73]. Since T ∩ Q needs not be empty, even if it does not contain integral points, there may be cases where the method ends up in the very unfortunate situation that a large portion of cuts needs to be generated to exclude all integral points. Lemma 6 Let N = (M, R, S) be a metabolic network with irreversible reactions D ⊆ R. Let ˙ − r ∈ R. Let X be a minimal cut set of N w.r.t. r. Then there exists a partition X = X + ∪X + − such that (X , X ) is an OMCS of N w.r.t. r. 59

CHAPTER 5. OPTIMIZATION BASED APPROACHES

˙ − there exists an elementary mode E with Proof Assume for every partition X = X + ∪X + + − − + E ∩ X = ∅ = E ∩ X and r ∈ E .

By Assumption 2, we can order the elements of X and write w.l.o.g. X = {r1 , r2 , . . . , rn }.

˙ − there exists an elemenClaim 5.4.2 For every i ∈ {0, . . . , n} and every partition X = X + ∪X tary mode E such that {r1 , . . . , ri } ∩ E = ∅,

E ∩ X + = ∅ = E− ∩ X −, r ∈ E+ +

Proof by induction on i. For i = 0 the claim equals the assumption, hence it is true. Assume now, the claim is true for i. We will now prove it for i + 1. Let X + , X − be an arbitrary but fixed partition of X. Define Y + := X + 4{ri+1 }, Y − := X − 4{ri+1 }, where 4 denotes the symmetric difference. By induction hypothesis, there exist elementary modes A, B s.t. {r1 , . . . , ri } ∩ A = ∅

{r1 , . . . , ri } ∩ B = ∅

A ∩ X + = ∅ = A− ∩ X − , r ∈ A+ +

B+ ∩ Y + = ∅ = B− ∩ Y −, r ∈ B+

If A ∩ {ri+1 } = ∅ or B ∩ {ri+1 } = ∅, we have found the elementary mode we were looking for.

Hence, we only need to analyze the case where ri+1 ∈ A and ri+1 ∈ B. By construction it follows that ri+1 ∈ A+ ∩ B − or ri+1 ∈ A− ∩ B + . With r ∈ A+ ∩ B + , we can apply circuit axiom (C3’) (see Theorem 1) and deduce that there also exists a circuit C of the flux mode matroid with C + ⊆ (A+ ∪ B + ) \ {ri+1 }

C − ⊆ (A− ∪ B − ) \ {ri+1 } r ∈ C. It follows that {r1 , . . . , ri } ∩ C ⊆ {r1 , . . . , ri } ∩ (A ∪ B) = ∅, ri+1 6∈ C. Since C only has signs that were used by A, B, it is also an elementary mode and C + ∩ X + = ∅ = C − ∩ X − . Since r ∈ C, it follows that r ∈ C + .  With i = n, we can conclude that there exists an elementary mode that does not contain any of the blocked reactions X but carries flux through r. This is a contradiction to X being a minimal cut set.  We now showed that the combinatorial cuts on the flow-polytope Q have a one-to-one correspondence to OMCSs; and the number of OMCSs explodes in practice. If we want to check if a given flow through the target reaction is possible, the cuts will of course not equal the OMCSs anymore. But, since the concept is very similar, it is very likely, that also in this case the number of cuts will be huge. Hence, I consider it unlikely that this approach will lead to an efficient practical method. 60

5.4. MIXED INTEGER LINEAR PROGRAMMING

Practical Applicability of Combinatorial Benders’ Cuts We successfully strengthened the MILP formulation by using Combinatorial Benders’ cuts. So this theoretical result should also be observable in practical experiments. I ran thermodynamically constrained FBA experiments in the E. coli iAF1260 model and maximized flux through the reaction with id NDPK1. I used NDPK1 because it is one of the reactions that is contained in internal cycles. The naive model I used is the model proposed by Schellenberger et al. [109] which is also used in the COBRA-Toolbox [112]. The main difference is that it encodes the space of feasible potential difference using the null space matrix N of SI : max JNDPK1 s.t. SJ = 0 ∆µN = 0 `Ji (1 − ai ) ≤ Ji ≤ uJi ai

−1000ai + (1 − ai ) ≤ ∆µi ≤ −ai + 1000(1 − ai ) a ∈ {0, 1}I

In addition I implemented a constraint handler that generated needed combinatorial Benders cuts. The performance of the combinatorial Benders’ cuts was then tested • by using only the constraint handler to generate cuts to enforce thermodynamic feasibility. • by pre-computing all combinatorial benders cuts (possible using elementary mode enumeration, see Section 4.2) and adding all these cuts to the initial unconstrained problem. I obtained the running times shown in Table 5.1. Formulation, Solver Naive model 1 in Gurobi Naive model 1 in Scip Combinatorial Benders’ cuts in Scip Combinatorial Benders’ cuts in Scip pre-generated Combinatorial Benders’ cuts in Gurobi pre-generated

Running time 5 sec. 556 sec. does not finish 1.53 sec 0.24 sec

Table 5.1: Running times for different formulations and solvers for maximization of NDPK1 in the E. coli iAF1260 model. The result is very interesting. First of all, Gurobi already solves the naive formulation very well; Scip on the other hand, has great difficulties in finding the optimal solution. Experiments with a smart primal heuristic (see Section 5.5.1) showed that a good heuristic is able to reduce the running time to 15 seconds in Scip. Note also that Gurobi is able to use all 4 CPU cores, while Scip is single-threaded. But, this is not the interesting part. Surprisingly, the constraint handler for the combinatorial Benders’ cuts does not improve the running time at all, where as the formulations with pregenerated cuts are extremely quick. I see only two explanations: Either my implementation of the constraint handler was not very good, or heuristics and branching strategies had to work without knowing the other cuts and hence, made bad branching decisions. For example, NDPK1 is contained in only 4 out of 38 internal circuits. We will see later that it will suffice to branch on 61

CHAPTER 5. OPTIMIZATION BASED APPROACHES

the 4 circuits that contain NDPK1 to obtain an optimal solution. If we, however, start branching on the other circuits, we end up doing much useless work. Both, Scip and Gurobi, only create 5-10 branching nodes in the formulation with pre-generated combinatorial Benders’ cuts. Apparently, the correct way of branching is important. To understand why combinatorial Benders’ cuts alone do not lead to a good formulation, consider the networks shown in Figure 5.6.

Figure 5.6: Benders Cuts do not restrict the flow through the objective arc much. If the objective arc is only contained in one circuit consisting of n reactions the flow through the objective arc is only reduced by n1 . If the objective arc is contained in two or more circuits, flux equal to the upper bound value is possible. We introduced the Benders’ cuts to disable internal circuits. For a positive circuit (C, ∅) of size n, the cut looked as follows (assuming w.l.o.g. that all reactions are in forward direction): a(C) ≤ |C| − 1 This means that in the LP relaxation, ai = |C|−1 |C| for all i ∈ C is feasible. From this, it follows that the cut only reduces the maximal circulation (possible in FBA without thermodynamic 1 constraints) by a factor of at most |C| . It gets even worse if the reaction is contained in 2 circuits. Then the same maximal flux through the reaction as in the unconstrained case can be achieved (see Figure 5.6). Hence, we need to branch on the internal circuits.

5.5

Constraint Programming

We observed that in the end, the IP solvers had to branch in a clever way on the internal circuits. But if we have to branch anyway, why do we need to compute the cuts in the first place? With constraint programming, we can enforce thermodynamic feasibility on the fly by implementing a branching routine. This has the main advantage that we do not need to precompute all internal cycles. I will now propose a new method for solving thermodynamically constrained flux problems. In constraint programming, we can add arbitrary constraints to our problem. The solver only requires that we also specify how the solver ought to enforce a constraint if it is violated. This means that we have to implement a constraint handler . In a cutting plane algorithm we told the solver to generate a cutting plane, but in constraint programming we can also tell the solver how to branch the current node to resolve the infeasibility. We do not even need to branch on integer variables, but we can also branch on continuous variables. 62

5.5. CONSTRAINT PROGRAMMING

Remember that in our thermodynamically constrained problem we introduced the boolean variable a only because we needed to somehow capture the sign of the flux variable. In constraint programming, we can directly use the flux variables to determine infeasibility and enforce feasibility. I implemented such a constraint handler that works basically as shown in Algorithm 3. Input: A flux J Find violated minimal infeasible set C = (C + , C − ) ⊂ sign(J) (see Theorem 13). if we found no infeasible set C then return feasible else for i ∈ C + do if `Ji ≤ 0 then Create child node with uJi = 0. end for i ∈ C − do if uJi ≥ 0 then Create child node with `Ji = 0. end Add created child nodes to branch and bound tree and continue solving nodes. end Algorithm 3: Constraint handler that enforces thermodynamic feasibility First we search for an infeasible set, which can be done in polynomial time by Theorem 13. This infeasible set can then be minimized iteratively by blocking reactions until a minimal infeasible set C is found (note that this will not necessarily give an infeasible set with minimum number of reactions). Since C is minimal, every reaction in C carries flux. We know that to obtain a feasible solution at least one reaction in C must not carry flux in the same direction as it does in the current solution. This can be achieved by simply setting the corresponding bound to zero and thus blocking flux in the current direction. Since we do not know which of the reactions we must block to obtain an optimal feasible solution, we create child nodes where in each child node one reaction of C is blocked.

5.5.1

A Heuristic

Remember that we can transform every steady state flux using Algorithm 1 into a thermodynamically feasible flux. Of course, we need not get an optimal solution, but (if we do not have any flux forcing reactions) we will at least get a feasible solution for our problem. Hence, Algorithm 1 can be used as a primal heuristic. First of all, we can run Algorithm 1 before we even start with the branch and bound. Even if the network has flux forcing reactions, we may be lucky and those reactions do not participate in internal cycles. Usually, only exchange reactions (like nutrient uptake reactions) are flux forcing. Thus, we very often already have a feasible solution that can be used by the IP solver during presolving. In many practical cases, the objective function does not involve flux variables of reactions in internal cycles. In these cases, Algorithm 1 finds an optimal solution and this optimal solution has even the same objective value as the optimal steady state flow (see Theorem 9). Hence, in this case the heuristic finds a solution with the same value as the LP-relaxation; thus, we (and 63

CHAPTER 5. OPTIMIZATION BASED APPROACHES

also the IP solver) have a proof that the solution is optimal. If the objective function does contain flux variables of reactions contained in internal cycles, we may be lucky and Algorithm 1 also gives an optimal solution. Usually, however, we do not know whether this solution is already optimal or not. Hence, we will still have to branch on the current node to prove optimality of the solution that we already found. When we run the heuristic also on the nodes of the branch & bound tree, it is this reason that makes it not very useful to run the heuristic on every node. We may even get optimal solutions in some nodes, but we will not know it, until these nodes are branched out. Hence, in practical applications, it is most useful to run the heuristic only, when the objective flux variables are not contained in any internal cycles (of the subproblem).

5.5.2

Parametrized Complexity

We observed that in general the thermodynamic flux problem is hard (NP-hard), but in certain cases (Observation 1) the problem is solvable in polynomial time. The condition of Observation 1 was that we do not have flux forcing reactions and no objective flux variables that are contained in internal circuits. In practice it is usually the case that there exist some internal circuits, but there are only few of them. Parametrized Complexity is the theoretical approach to understand these kinds of problems (see for example [31, 45]). A typical example for parametrized complexity analysis is SAT. SAT can be solved by brute force in O(2k m) time, where k is the number of variables and m is the number of clauses. This means, that if we have few variables, but many clauses, we can solve SAT efficiently. The decision problem to the thermodynamic flux optimization problem is the following: Problem 3 (Thermodynamic Flux Optimization) Given: ˙ S) with flux bounds `, u ∈ RR • A Metabolic Network N = (M, R = I ∪E, • A linear objective function c : RR → R. • A target flux jopt ∈ R

Question: Does there exists a thermodynamically feasible flux J with ` ≤ J ≤ u such that cJ > jopt ? Parametrized tractability highly depends on the parameter. If, for example, k counts the number of flux variables of internal reactions in the objective function, the problem is not tractable in k. We have seen that already for k = 1, we can find instances where it is NP-hard to decide if there exists a flux through the reaction (Theorem 5). On the other hand, let k count the number of circuits containing important reactions in the network that are defined as follows: Definition 26 (Important Reaction) A reaction r is called important, if it is flux forcing (by Definition 22) or if cr 6= 0. 2 Using Theorem 9 and the observations about the heuristic, we get 64

5.5. CONSTRAINT PROGRAMMING

Theorem 14 If the number of circuits containing important reactions is bounded, the thermodynamically constrained flux optimization problem (Problem 3) can be computed in polynomial time. Proof If none of the important reactions is contained in a circuit, we can compute the optimal flux by Theorem 9 in polynomial time, i.e. in time O ((|R| + |M|)p ) for some constant p.

By brute force we can compute all the fluxes where one reaction of each circuit that contains an important reaction is blocked. Assume the number of circuits containing important reactions is bounded by k. Thus, there are at most |R|k possibilities to block one reaction of each circuit. We conclude that Problem 3 can be solved in time O |R|k (|R| + |M|)p




Parametrized tractability results are, of course, only useful in a practical context if the parameter is also small in practice. This is indeed the case here. I used the metabolic networks that are available from the BiGG database [111]. For all networks, except H. sapiens recon 1, I was able to compute the set of internal circuits. The human reconstruction has 958 reactions that are involved in internal circuits. This amount was probably to large for the elementary mode algorithm by Terzer [125]. Also the improved algorithm of Wright et al. is not able to enumerate the internal circuits of the human reconstruction [136]. In all metabolic networks that I tested (except the H. sapiens network for which I wasn’t able to compute the circuits), there are only few, small internal circuits, as can be seen in Figure 5.7. This is already a very nice property of the network. If we were unlucky and the property is not satisfied by the network, we might have to test 2|R| possibilities of reaction directions. With this property we already get, by Theorem 14, a much better bound on the complexity. On the other hand 38 circuits (as for example in the E. coli iAF1260 model) are still quite a lot. Underestimating the size of the circuits, by assuming that all circuits have size 2, we still get 238 ≈ 2.7 · 1011 possibilities. Hence it is still not practical to branch through all possibilities of blocking reactions for blocking all internal circuits. But, by Theorem 14, we only need to consider the circuits containing important reactions. For the case that only one reaction is important (a typical scenario in FVA), we can see in Figure 5.8 that even less circuits are involved. Many reactions are only contained in 1 or 2 circuits. The most difficult case (in the number of circuits that contain a reaction) occurs in M. tuberculosis, where we have reactions that are contained in up to 11 circuits. By also looking at the sizes of the circuits, we can bound the number of cases that need to be analyzed to at most 62208. The number of possibilities is relatively small, because the circuits that are involved, consist of at most 5 reactions, as can be seen from Figure 5.7. In the case of S. cerevisiae, we have many large circuits (containing up to 13 reactions). For S. cerevisiae only an upper bound of approx. 6.58 · 108 can be shown by multiplying the cardinalities of the circuits that contain a common reaction (note, that we must not branch on the objective reaction, since it must be able to carry flux; this reduces the branching candidates per circuit by 1). We conclude that for most metabolic network that are encountered in practice, we can already apriori predict that the constraint programming method with Algorithm 3 as constraint handler for resolving thermodynamic infeasibility and the discussed heuristic (Algorithm 1) will compute a solution in a practically useful amount of time. 65

CHAPTER 5. OPTIMIZATION BASED APPROACHES

35 E. coli iJR904 E. coli iAF1260 E. coli iJO1366 H. pylori iIT341 M. barkeri iAF692 M. tuberculosis iNJ661 S. aureus iSB619 S. cerevisiae iND750

30

number of circuits

25

20

15

10

5

0 1

2

3

4

5

6 7 8 size of circuit

9

10

11

12

13

Figure 5.7: Sizes of the internal circuits in various metabolic network models

5.6

Performance Comparisons

As discussed, there exist several different ways on how to formulate thermodynamic constraints. Since the feasible domain of strong thermodynamic constraints is not closed, we will only compare running times of various implementations for weak thermodynamic constraints. Since the nonlinear programming formulations like the ones in [63, 42] do only give local optima, we will only compare running times to MILP formulations like the ones in [109, 27, 64] with the CP formulation introduced in the previous sections. The CP formulation was implemented with SCIP [1] in C++. For loading the metabolic network models, libSBML [21] is used. Both, libSBML and SCIP, are open source, and can be freely used for noncommercial purposes. For the CP formulation I implemented the constraint handler and heuristic mentioned in the previous sections. This code is available at http://page.mi. fu-berlin.de/arnem/fast-tfva.html and there also exists an interface for Matlab such that Matlab users can easily use the implementation.

5.6.1

Optimization on Single Reactions

First we will analyze the performance of the three different MILP formulations and my constraint programming formulation for optimization on a single objective reaction. As we have already seen from the theory part, different kinds of reactions pose different degrees of difficulty. Hence, we will not analyze just one reaction, but reactions with different difficulty. The difficulty will be measured in the number of internal circuits the reaction is contained in and in the theoretical upper bound on the number of cases that need to be analyzed by the brute force method described in Section 5.5.2. 66

5.6. PERFORMANCE COMPARISONS

40 E. coli iJR904 E. coli iAF1260 E. coli iJO1366 H. pylori iIT341 M. barkeri iAF692 M. tuberculosis iNJ661 S. aureus iSB619 S. cerevisiae iND750

35

number of reactions

30 25 20 15 10 5 0 1

2

3 4 5 6 7 8 9 number of circuits containing the reaction

10

11

Figure 5.8: number of reactions over the number of internal circuits each reaction is contained in

MILP 1 This is the formulation used by Schellenberger et al. [109]. K is the nullspace matrix of SI and M = 1000.

max cJ s.t. SJ = 0 K∆µ = 0 `J ≤ J ≤ uJ

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + ε ≤ ∆µi ≤ M (1 − ai ) − ε ∀i ∈ I a ∈ {0, 1}I

MILP 2 Instead of using the null space matrix, we can also directly formulate the space of potential differences using the potentials µ. A formulation of this kind is used by Hoppe et al. 67

CHAPTER 5. OPTIMIZATION BASED APPROACHES

[64] and Cogne et al. [27]. max cJ s.t. SJ = 0 ∆µ = µSI

`J ≤ J ≤ u

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + ε ≤ ∆µi ≤ M (1 − ai ) − ε ∀i ∈ I `µ ≤ µ ≤ uµ a ∈ {0, 1}I

MILP 3

This formulation is obtained from MILP 3 by eliminating the ∆µ variables. max cJ s.t. SJ = 0 `J ≤ J ≤ eJ

`Ji (1 − ai ) ≤ Ji ≤ uJi ai ∀i ∈ I

−M ai + ε ≤ µSi ≤ M (1 − ai ) − ε∀i ∈ I `µ ≤ µ ≤ uµ

a ∈ {0, 1}I

CP Is the constraint programming approach that is proposed in this thesis. Reaction CA2t3pp PPKr NDPK1 ACKr ACCOAL ABUTt2pp CYNTAH Biomass

number of circuits 6 5 4 3 2 1 0 0

theoretical bound 729 1680 100 140 4 2 1 1

CP 2.37s 0.86s 1.03s 0.78s 0.8s 0.9s 0.71s 1.01s

MILP 1 10.2s 0.78s 5.06s 1.19s 1.99s 4.17s 0.50s 1.93s

MILP 2 42.5s 5.72s 7.40s 7.13s 5.25s 6.33s 5.09s 10.8s

MILP 3 29.5s 4.86s 10.9s 7.87s 5.92s 5.19s 7.70s 18.0s

Table 5.2: Running times for maximizing flux through selected reactions of the E. coli iAF1260 model with different theoretical complexity bounds. The MILP formulations were run using gurobi on a Intel Core i5-2400S (2.5 GHz, 4 cores, 6 MiB L3 cache). The CP approach was implemented using SCIP 2.0.1 and run on a Intel Core i3-2310M (2.1 GHz, 2+2 cores, 3MiB L3 cache) In Figure 5.9 we can see that indeed the running time correlates with the theoretical upper bound. However, the deviations are very huge. For example, there are reactions with a low theoretical complexity bound that need more running time than some of the reactions with high theoretical complexity bound. Since these deviations can also be observed for the constraint programming approach, it is plausible that the theoretical bound is simply a bad estimation in 68

5.6. PERFORMANCE COMPARISONS

2

10

Milp 1 Milp 2 Milp 3 CP 1

running time

10

0

10

−1

10

0

10

1

10

2

10 theoretical upper bound

3

10

4

10

Figure 5.9: For each reaction, corresponding to its running time and theoretical bound, a dot is drawn (color according to the method). The lines are linear regression lines through these points. Because this plot also shows reactions that need no branching, the theoretical bound was increased by one for each reaction.

69

CHAPTER 5. OPTIMIZATION BASED APPROACHES

several cases. For example, if a reaction is used by several circuits, we would not have to check as many possibilities as the bound would estimate. We also see that the theoretical work still gives a nice reward: The constraint programming algorithm that explicitly uses the property of few internal circuits clearly outperforms the other formulations. In particular, in the most difficult case, the constraint programming approach is more than 4 times faster than the best one of the other approaches (see Table 5.2). But there exists also one case, where the constraint programming approach is a bit slower than the MILP formulation using the null space matrix (MILP 1). MILP 1 also outperforms the other MILP approaches - very often also with a huge running time difference. This is remarkable, since the three MILP formulations are only affine transformations of each other.

5.6.2

Flux Variability Analysis on Networks of Different Sizes

As we were able to see from Table 5.2 on MILP formulation running times, MILP 1 performed best. Hence, we will use the implementation by Schellenberger et al [109] as a reference. The implementation by Schellenberger et al. also has the advantage, that it is part of the COBRAToolbox which is published under an open source license [112]. The COBRA-Toolbox already has flux variability analysis implemented. All experiments are performed using metabolic models from the BiGG-Database [111]. At the time of writing the BiGG database contained models of • E. coli (iJR904, iAF1260), • H. pylori (iIT341), • H. sapiens (Recon. 1), • M. barkeri (iAF962), • M. tuberculosis (iNJ661), • S. aureus (iSB619), • S. cerevisiae (iND750), and a few very small networks (like a textbook version of E. coli ). We will compare the running times of flux variability analysis on the various metabolic networks from the BiGG-database. Although Schellenberger et al. claim to have run FVA successfully on the H. sapiens model [109], the current implementation of the COBRA Toolbox (version 2.0.3) throws an error on the H. sapiens network. The H. sapiens network is very problematic, since it contains many reactions that are involved in internal circuits. Also the implementation proposed in this thesis does not terminate within 2 hours of computation time. Flux variability analysis is usually applied to analyze the space of optimal solutions or solutions that are close to optimal. We will apply flux variability analysis to • the flux space of optimal fluxes (100%), • the flux space of at least 90% optimal fluxes, 70

5.6. PERFORMANCE COMPARISONS

• the flux space of at least 50% optimal fluxes, • the flux space of at least 10% optimal fluxes, and

seconds / (no. iterations * reactions in int. circuits)

1

E. coli iJO1366

E. coli iAF1260

S. cerevisiae iND750

M. tuberculosis iNJ661 E. coli iJR904

M. barkeri iAF692 S. aureus iSB619

H. pylori iIT341

• the complete flux space (0%).

CP COBRA

0.1

0.01

0.001

0.0001 1000

1500 no. reactions

2000

2500

Figure 5.10: The running time per iteration (optimization of flux through a reaction) also depends on the number of reactions in internal circuits. Hence, this plot shows the average running time per iteration divided by the number of reactions. All FVA times used in this plot are were taken for 0% optimal fluxes. The constraint programming approach is between 50 and 100 times faster than the MILP approach by Schellenberger et al. (see Table 5.3 and Figure 5.10). This is not explicable by the results that could be seen in Figure 5.9, alone. Indeed, one other property plays an important role in the achieved speedup. In many cases, it is sufficient to solve a single LP (as in FBA without thermodynamic constraints). This is for example the case if the computed LP solution does not contain any important circuits already (see Theorem 9). Using warm starting for the LP solver (we do not compute the LP solution from scratch, but use the feasible solution of a previous iteration), we gain an enormous speedup for many reactions. This was already observed in the context of FVA without thermodynamic constraints by Gudmundsson and Thiele [59]. 71

CHAPTER 5. OPTIMIZATION BASED APPROACHES

100000

E. coli iJR904 (CP) E. coli iAF1260 (CP) E. coli iJO1366 (CP) H. pylori iIT341 (CP) M. barkeri iAF692 (CP) M. tuberculosis iNJ661 (CP) S. aureus iSB619 (CP) S. cerevisiae iND750 (CP) E. coli iJR904 (COBRA) E. coli iAF1260 (COBRA) E. coli iJO1366 (COBRA) H. pylori iIT341 (COBRA) M. barkeri iAF692 (COBRA) M. tuberculosis iNJ661 (COBRA) S. aureus iSB619 (COBRA) S. cerevisiae iND750 (COBRA)

10000

Seconds

1000

100

10

1 0

20

40 60 Optimality

80

100

Figure 5.11: Running Times of FVA on different metabolic networks. For every network a different degree of optimality was required. The 0% optimal fluxes was FVA over the whole flux space. The 100% optimal fluxes are the fluxes, where optimal flux through the biomass reaction is achieved. Another interesting observation is the increased running time of at least 10% optimal fluxes (Figure 5.11). In the constraint programming approach this is only slightly visible, but in the COBRA implementation, the increase is significant. In the case of E. coli iJR904, the running time nearly doubles. Apparently, the condition that there needs to be some flow through the network makes the problem harder. The theoretical approach described in this thesis does not explain this behavior. The running time increase may be an artifact of the IP-solver “magic” that gurobi performs, but there also exists the possibility that this effect can be captured by a parameter and included in theoretic complexity analyses of thermodynamic constrained flux problems.

72

73

CP COBRA CP COBRA CP COBRA CP COBRA CP COBRA CP COBRA CP COBRA CP COBRA CP COBRA

method

43

3

66

31

?

8

46

38

19

internal circuits

1266 (64)

743 (7)

1025 (53)

690 (30)

3742 (958)

554 (22)

2583 (76)

2382 (68)

1075 (40)

reactions (in int. circuits)

12s 204s 38s 844s 5s 179s 74s 1662s

19s 963s 109s 11162s 138s 25869s 7s 122s > 2h

0%

24s 327s 46s 1444s 5s 221s 84s 3763s

21s 1852s 112s 17865s 140s 48100s 7s 197s > 2h

10% 21s 1706s 112s 18875s 146s 42476s 7s 151s > 2h Error 23s 260s 48s 1396s 4s 234s 86s 4219s

50%

21s 309s 47s 1495s 4s 197s 83s 4615s

20s 1732s 110s 17367s 137s Error 7s 162s > 2h

90%

12s 332s 28s 1316s 3s 198s 51s Error

18s Error 77s 15049s 114s Error 6s 176s > 2h

100%

Table 5.3: Running Times of FVA on different metabolic networks. For every network a different degree of optimality was required. 0% optimality means that we performed FVA over the whole flux space. The 100% optimal fluxes are the fluxes, where optimal flux through the biomass reaction is achieved.

S. cerevisiae iND750

S. aureus iSB619

M. tuberculosis iNJ661

M. barkeri iAF692

H. sapiens Recon. 1

H. pylori iIT341

E. coli iJO1366

E. coli iAF1260

E. coli iJR904

Model

5.6. PERFORMANCE COMPARISONS

CHAPTER 5. OPTIMIZATION BASED APPROACHES

74

Chapter 6

Future Work In this thesis we mostly discussed the case where the objective function only contains flux variables. For future research we can go further down this track and try to solve also difficult instances (like the human metabolic network reconstructions) efficiently. On the other hand biologists are increasingly interested in also deriving information about metabolite concentrations. In a certain aspect this is also possible using optimization methods containing thermodynamic constraints.

6.1

Optimization over Fluxes

All networks we analyzed only contained few internal circuits and thus were easy to analyze - except the human metabolic network reconstruction (H sapiens recon 1 ). This means that in practice we may indeed encounter networks that are hard to analyze using only the tools discussed in this thesis (although I was not able to validate it, Schellenberger et al. claim that their implementation was able to do FVA on H. sapiens recon 1 in roughly 32 hours [109]). To improve the running times, the following approaches could be considered: • Improving heuristics: By looking more closely on the reactions that are problematic in the H. sapiens network, we can see that in many cases, the solver has great difficulties finding the optimal feasible solution. This means that the solver does not spend much time in proving optimality (which is coNP-hard and hence cannot be done by smart guessing [unless coNP = NP]). An improved heuristic (that ideally immediately guesses the optimal solution) would lead to a very short running time for these reactions. We could implement an additional diving heuristic (see chapter 12.5 in [135]): The heuristic proposed in this thesis only subtracts flux through internal circuits. This has the effect that the flux of one reaction through the circuit is set to zero and thus, there remains no flux through the circuit. Instead of subtracting the flux through the internal circuit, a heuristic can instead set the flux through one of the circuit reactions to zero (a good choice should be the reaction that would have been set to zero flux by the subtraction operation). Then the flux through the remaining network is recomputed and the procedure repeated until no internal circuit remains. • Improving the dualbound: If the network contains a large complex of circuits containing 75

CHAPTER 6. FUTURE WORK

a reaction and there is no flux possible through this reaction, good heuristics are of no use in this case. Unless NP = coNP (which is considered unlikely [66]), there even exists a sequence of networks such that no sequence of witnesses (of optimality) exist that can be checked in polynomial time w.r.t. the size of the corresponding network. But, of course we only want to solve practical problems, thus we may be lucky and may be able to find witnesses that prove optimality in polynomial time. In this thesis, we only worked with the LP relaxation. Since we are actually doing constraint programming, we may be able to develop a sharper relaxation. For example if the objective reaction is part of a circuit, we could try to contract the circuit in such a way that all fluxes (except the flux in the circuit) through the objective reaction are still feasible. In the case of directed graphs, this operation is illustrated in Figure 6.1.

Figure 6.1: A possibility for contracting a circuit in a directed graph containing an objective reaction (thick arc). Observe that flux entering by arc e and leaving by f did not contribute to the objective in the original model, but does so in the relaxation. Another possible relaxation method is to transform the metabolic network into a directed graph. This way, we loose coupling information, but for digraphs there exist measures like directed tree width (Johnson et al. [68]), D-width (Safari [105]), DAG-width (Berwanger et al. [14]), and many more [65, 49] that in some sort of way measure the acyclicity of digraphs. It can be shown that the disjoint path problem is fixed parameter tractable w.r.t. to most of these acyclicity measures (see for example [68]). We have seen in Section 3.2 that deciding whether there exists a thermodynamically feasible flux (in directed graphs) through a given reaction is equivalent to deciding the disjoint path problem. Hence, it is plausible that it is also possible to extend the tractability results to thermodynamically constrained flows. This could then be used to solve the thermodynamically constrained flux problem efficiently on the relaxed network to sharpen the flux bound on the original thermodynamically constrained flux problem in the metabolic network. • Derive stronger tractability results: The tractability result (Theorem 14) we proved in this thesis has the main practical disadvantage that the number of internal circuits of a network (the tractability parameter itself) may grow exponentially with the size of the network. Hence, it would by nice, if concepts like directed tree-width, D-width, etc. [68, 105, 14, 65, 49] that are only defined for directed graphs could be translated to metabolic networks. This may lead to other tractability results for the thermodynamically constrained flux problem. 76

6.2. OPTIMIZATION OVER POTENTIALS, INTEGRATED MODELS

All computational methods developed in this thesis only dealt with the case that nothing is known about metabolite concentrations. We have seen (see Theorem 3) that we can use similar branching techniques as we used for internal circuits in the more general setting. But, we will not be able to use Theorem 9 anymore to solve a certain class of problems efficiently. To obtain similarly efficient algorithms we need to generalize the result to systems with a more general class of infeasible reaction sets. In very complex cases, it may also be the case that the original, naive MILP formulation is more efficient than the constraint programming approach. In such a case, it may be worthwhile to analyze hybrid systems that use both constraint programming and part of the MILP formulation (e.g. for difficult subnetworks).

6.2

Optimization over Potentials, Integrated Models

Considering works by Cogne et al. [27], Hoppe et al [64], Heuhett et al. [63] and many others [11, 138, 76], we can see great interest and effort in also obtaining information about potential concentrations using thermodynamically constrained analysis. This is useful, because from the metabolite potentials information on the regulatory control of the cell can be obtained [76]. It should be mentioned at this point, that integration of regulatory and transciptory information into metabolic networks is done as well, too [29, 130]. In terms of optimization programs this means that we have to deal with objective functions that also contain potentials. This is not only problematic in a computational sense. Remember the original physical formulation of pairs (J, µ) of thermodynamically feasible fluxes J w.r.t. µ. This formulation contained strict inequalities. In the case of weakly thermodynamically feasible fluxes where we only required the existence of a corresponding vector of potentials µ, we were lucky. Since we were not interested in µ, we could project out µ. The resulting domain of feasible fluxes was closed.

J4, ∆µ4

J1

∆µ3

∆µ4

J3

J4

J2

J1 = 1 = J2

J3, ∆µ3

∆µ3 + ∆µ4 = 1

Figure 6.2: In this example the flux through reactions 1, 2 is constrained to 1 and the potential difference between the left most and right most metabolite is 1. Such constraints are not unrealistic [87]. The feasible flux space is denoted by the black point. Because J4 , J3 > 0, ∆µ3 , ∆µ4 > 0. It follows that the feasible potential space (red line) is open. Now, we want to optimize on µ. Hence we cannot project out µ and thus the old trick does not work. Indeed, for a fixed J the domain of feasible µ is not closed (see Figure 6.2). In practice this is of course no reason to not optimize on such sets. In particular, since we may also have additional 77

CHAPTER 6. FUTURE WORK

bounds on the metabolite potentials, it is possible that solutions of such practical optimization problems exist. Instead of the possible result cases of ordinary minimization programs, we have to expect answers of the types • minimal solution found • problem is infeasible (no feasible solution exists) • problem is unbounded • infimum of solution found, but no minimum exists (the new case) For programs containing strict linear inequalities there already exists work by Kuhn from 1956 [75]. This work was recently extended by Goberna and Rodr´ıguez [54, 53] for systems containing linear inequalities, but infinitely many. Non-convex optimization problems of a specific form, with one strict inequality were analyzed by Lemaire and Volle [5]. One of the difficulties in such problems lies in finding the closure of a the feasible domain. For example, it does not hold in general for arbitrary sets A, B that cl(A) ∩ cl(B) = cl(A ∩ B). Goberna and Rodr´ıguez introduce the concept of evenly convex hulls [54] to decide whether a given system (of possibly infinitely many inequalities) contains a feasible solution. Since thermodynamic constraints are non-convex, a theory needs to be build that can handle the non-convexity of thermodynamically feasible fluxes and potentials. With such a theory available, we can then build algorithms that also solve optimization problems in the non-closed domain of thermodynamically feasible fluxes and potentials. If this is not difficult enough, we can start incorporating kinetic data into the steady state models. Fleming and Thiele [40] showed that for mass-action kinetics there always exists a steady state flux. Usually, kinetic data is not available, but in some special cases such data is already available today. Liebermeister and Klipp [80] proposed a general approach called convenience rate law , which could be used as theoretical starting ground for integration of kinetic parameters. This leads to additional nonlinear constraints . Currently there already exists a wide variety of solvers for mixed integer nonlinear programs and also some solvers like alphaBB, BARON, Couenne, LindoGlobal and Scip are also able to deal with nonconvex domains and objective functions [24]. Because MINLPs are a very general class of problems, targeted exploiting of structural properties will very probably lead to significantly more efficient algorithms.

6.3

Conclusion

In this thesis we discussed different notions of thermodynamic feasibility: strongly thermodynamically feasible fluxes and weakly thermodynamically fluxes. For weakly thermodynamically feasible fluxes, we analyzed the space of feasible fluxes in great detail and have seen methods that can be used to check thermodynamic feasibility efficiently. We showed that it is hard to decide whether there exists a thermodynamically feasible flux through a given reaction (for elementary modes this was easy). For one of the most often used tools, FBA, we discussed how thermodynamic constraints can be integrated efficiently. The performance of the algorithm was shown theoretically by a parametrized tractability result and also experimentally by a series of flux variability analyses on metabolic networks of different sizes. Furthermore, this is not everything that can be done and indeed the demands of system biologists reach higher: Integration of information on metabolite potentials and kinetic data, as well as optimization on metabolite potentials can make the problem arbitrary difficult and challenging to analyze. 78

Appendix Notation Here I summarize all the notation I use in the thesis: • vi denotes the i.th element of a vector v. • vA denotes the vector containing only the elements of the index set A. • Sij denotes the element at row i and column j of S. • Si∗ denotes the i.th row of matrix S. • S∗i denotes the i.th column of matrix S.

• SA denotes the matrix containing only the elements of the index set A ⊂ N2 . If it is clear from the context, A may also be of one dimension and then SA only contains the columns respectively the rows indexed by A. • 1 denotes the all ones vector of appropriate size. P • v(A) := i∈A vi = 1vA . • X ∈ {−, 0, +}E is called a signed subset of E. • (X + , X − ) with X + ∩ X − = ∅ is another notation of a signed set X ∈ {−, 0, +}E , where X + = {i ∈ E : Xi = +} and X − = {i ∈ E : Xi = −}. ˙ − is the support of a signed set (X + , X − ). • X := X + ∪X

• sign : RE → {−, 0, +}E is the component wise sign operation, see Definition 7. • sign : {0, 1}E → {−, +} assigns component wise 0 7→ − and 1 7→ +. • (X + , X − ) ⊆ (Y + , Y − ) ⇔ X + ⊆ Y + ∧ X − ⊆ Y − . • supp(v) denotes the support of v.

79

APPENDIX

Usual names of variables are: • R for the set of reactions, • M for the set of metabolites, • I for the set of internal reactions, • E := R \ I for the set of exchange reactions, • S for the stoichiometric matrix, • C set of circuits (either internal or all), • D for the set of irreversible reactions, • J for the flux vector, • µ for the potentials, • ∆µ for the potential differences, • c for the cost function, • ` for lower bounds, • u for upper bounds, • `J for lower bounds on the fluxes, • uJ for upper bounds on the fluxes, • `µ for lower bounds on the potentials, • uµ for the upper bounds on the potentials

80

Index ACHR, 30 artificial centering hit-and-run, 30 ATP, 27

evenly convex hull, 78 exchange reaction, 10, 40, 63 extreme ray, 31

backward flux, 43 Benders’ cut, 52 big-M, 50 biomass, 39 bottleneck, 20 branch and bound, 48

FBA, 39, 61 feasibility classifier, 36 finitely constrained, 31 finitely generated, 31 flow cover cut, 50 flow cover inequality, 51 flux, 10 flux balance analysis, 39, 61 flux cone, 25 flux coupling, 29 flux forcing, 40, 63 flux mode matroid, 14, 35, 59 flux polyhedron, 39 flux variability analysis, 27, 30, 39, 70 flux variance analysis, 29 fluxinvariance, 19 forward flux, 43 FVA, 39

circuit, 12, 36, 59 circuit axiom, 12 circuit inequality, 54 closed, 26 cocircuit, 15 cocycle, 15 combinatorial Benders’ cut, 52 concentration, 19, 77 concentrations, 75 constraint handler, 62 constraint programming, 62, 66 convenience rate law, 78 convex, 26 coupling, 29 covector, 15 CP, 62, 66 cutting plane, 50 cycle, 12, 15

generating vector, 31 generator, 31 group contribution, 19 growth maximization, 39 heuristic, 63, 75 IIS, 58 important reaction, 64 incumbent, 48 infeasible set, 21, 56 inner description, 31 integer program, 58 integrated model, 77 internal cycle, 16, 63 internal cycle matroid, 14 internal reaction, 10 IP, 58

digraph, 28 directed graph, 28 directed tree width, 29, 76 Disjoint-Path, 28 dual bound, 48, 75 dual oriented matroid, 15 elementary flux mode, 31, 32 elementary mode, 14, 31, 32, 39, 59 elementary mode cone, 32 EM, 32 81

INDEX

irreducible inconsistent subsystem, 58 irreversible, 25, 32

realizable oriented matroid, 13 recession cone, 32 reversible, 25, 32

kinetic, 78 sampling, 29 separation, 50, 56 sequential quadratic programming, 44 set covering polytope, 56 sign vector, 12 signed set, 20 steady state, 10 stoichiometric matrix, 9 strict inequality, 26, 78 strongly thermodynamically feasible, 10, 26 subnetwork, 20

lineal, 32 local search, 56 LP relaxation, 48 LP-duality, 17, 20 LP-relaxation, 63, 76 Markov Chain Monte Carlo, 30 MCS, 58 metabolic network, 9 metabolite, 9 metabolite concentration, 19 metabolite potential, 19 MILP, 45, 66, 77 minimal cut set, 58 minimal infeasible set, 56 minimal metabolic behavior, 32 minimally infeasible set, 20 MIS, 20, 56 mixed integer linear program, 30, 77 mixed integer linear programming, 45, 66 mixed integer quadratic program, 30

thermodynamically feasible, 10, 11 Thermoflux, 27 tractability, 29, 76 two vertex disjoint path problem, 28 type III pathway, 17 vector (of an oriented matroid), 15 weakly thermodynamically feasible, 11, 27 yield maximization, 39

nonlinear constraints, 78 nonlinear programming, 42 OMCS, 59 optimality cut, 57 oriented cut set, 58 oriented matroid, 12, 33, 59 oriented matroid duality, 15 oriented minimal cut set, 59 orthogonality, 16 outer description, 31 parametrized complexity, 64 perturbation lemma, 16 pointed cone, 31 polyhedral cone, 25, 31 polyhedron, 25, 31 polytope, 25, 31 potential, 19, 77 potential perturbation lemma, 16 primal heuristic, 48, 63 reaction, 9 82

Bibliography [1] T. Achterberg. SCIP: Solving constraint integer programs. Mathematical Programming Computation, 1(1):1–41, 2009. http://mpc.zib.de/index.php/MPC/article/view/4. [66] [2] V. Acu˜ na, F. Chierichetti, V. Lacroix, A. Marchetti-Spaccamela, M.-F. Sagot, and L. Stougie. Modes and cuts in metabolic networks: complexity and algorithms. BioSystems, 95:51–60, 2009. [33] [3] V. Acu˜ na, A. Marchetti-Spaccamela, M.-F. Sagot, and L. Stougie. A note on the complexity of finding and enumerating elementary modes. BioSystems, 99:210–214, 2010. [33] [4] R. A. Alberty. Thermodynamics of Biochemcial Reactions. Massachusetts Institute of Technology, Cambridge, MA, 2003. [4, 10] [5] L. B. and V. M. A general duality scheme for nonconvex minimization problems with a strict inequality constraint. Journal of Global Optimization, 13:317–327, 1998. [78] [6] C. G. Bailey, D. W. Gull, and J. S. Oliveira. Hypergraphic oriented matroid relational dependency flow models of chemical reaction networks. arXiv:0902.0847v1 [math.CO], 2009. [35] [7] E. Balas and S. M. Ng. On the set covering polytope: I. all the facets with coefficients in {0, 1, 2}. Mathematical Programming, 43:57–69, 1989. [56] [8] K. Ballerstein, A. von Kamp, S. Klamt, and U.-U. Haus. Minimal cut sets in metabolic networks are elementary modes in a dual network. Bioinformatics, 2011. in press, doi:10.1093/bioinformatics/btr674. [58, 59] [9] D. A. Beard, E. Babson, E. Curtis, and H. Qian. Thermodynamic constraints for biochemical networks. Journal of Theoretical Biology, 228:327–333, 2004. [i, iii, 2, 6, 10, 11, 16, 27, 45, 53]

[10] D. A. Beard, S. dan Liang, and H. Qian. Energy balance for analysis of complex metabolic networks. Biophysical Journal, 83:79–86, 2002. [4, 10, 27] [11] D. A. Beard and H. Qian. Thermodynamic-based computational profiling of cellular regulatory control in hepatocyte metabolism. American Journal of Physiology - Endocrinology and Metabolism, 288:E633–E644, 2005. [1, 5, 19, 77] [12] D. A. Beard and H. Qian. Relationship between thermodynamic driving force and one-way fluxes in reversible processes. PLoS ONE, 2007. [43] 83

BIBLIOGRAPHY

[13] T. Berthold, S. Heinz, and S. Vigerske. Extending a CIP framework to solve MIQCPs. preprint. http://vs24.kobv.de/opus4-matheon/frontdoor/index/index/docId/606. [42]

[14] D. Berwanger, A. Dawar, P. Hunter, and S. Kreutzer. DAG-width and parity games. Lecture Notes in Computer Science, 3884/2006:524–536, 2006. [76] [15] L. Biegler and V. Zavaka. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization. Computers & Chemical Engineering, 33(3):575–582, 2009. [42] [16] Y. Bilu, T. Shlomi, N. Barkai, and E. Ruppin. Conservation of expression and sequence of metabolic genes is reflected by activity across metabolic states. PLoS Computational Biology, 2(8):e106, 2006. [39] [17] A. Bj¨ orner, M. L. Vergnas, B. Sturmfels, N. White, and G. M. Ziegler. Oriented Matroids. Cambridge University Press, 1999. [6, 12, 13, 15] [18] R. Bland. Complementary orthogonal subspaces of Rn and orientability of matroids. PhD thesis, Cornell University, 1974. [11] [19] B. A. Boghigian, H. Shi, K. Lee, and B. A. Pfeifer. Utilizing elementary mode analysis, pathway thermodynamics, and a genetic algorithm for metabolic flux determination and optimal metabolic network design. BMC Systems Biology, 2010. [23] [20] A. Bordbar, N. E. Lewis, J. Schellenberger, B. O. Palsson, and N. Jamshidi. Insight into human alveolar macrophage and m. tuberculosis interactions via metabolic reconstructions. Molecular Systems Biology, 6:422, 2010. [29] [21] B. J. Bornstein, S. M. Keating, A. Jouraku, and H. M. LibSBML: An API library for SBML. Bioinformatics, 24(6):880–881, 2008. [66] [22] A. Brøndsted. An Introduction to Convex Polytopes. Graduate Texts in Mathematics. Springer, 1982. [31] [23] A. Burgard, P. Pharkya, and C. Maranas. Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnology and Bioengineering, 84:647–657, 2003. [1] [24] M. R. Bussieck and S. Vigerske. MINLP solver software. Internet, 2010. http://vs24. kobv.de/opus4-matheon/frontdoor/index/index/docId/684. [78] [25] R. Byrd, J. Nocedal, and R. Waltz. KNITRO: An integrated package for nonlinear optimization. Nonconvex Optimization and Its Applications, 83:35–59, 2006. [42] [26] G. Codato and M. Fischetti. Combinatorial Benders’ cuts for mixed-integer linear programming. Operations Research, 54(4):756–766, 2006. [52, 58] [27] G. Cogne, M. R¨ ugen, A. Bockmayr, M. Titical, C.-G. Dussap, J.-F. Cornet, and J. Legrand. A model-based method for investigating bioenergetic processes in autotrophically growing eukaryotic microalgae: Application to the green algae Chlamydomonas reinhardtii. Biotechnol Progress, 27(3):631–640, 2011. [5, 19, 30, 45, 66, 68, 77] 84

BIBLIOGRAPHY

[28] A. R. Conn, N. I. Gould, and P. L. Toint. Lancelot: A Fortran package for large-scale nonlinear-optimization (Release A), volume 17 of Springer Series in Computational Mathematics. Springer-Verlag, 1992. [42] [29] M. W. Covert, C. H. Schilling, and B. O. Palsson. Regulation of gene expression in flux balance models of metabolism. Journal of Theoretical Biology, 213(1):73–88, 2001. [3, 77] [30] K. Dettmer, P. A. Aronov, and B. D. Hammock. Mass spectrometry-based metabolomics. Mass Spectrom Rev, 26:51–78, 2007. [2] [31] R. G. Downey and M. Fellows. Parametrized Complexity. Monographs in Computer Science. Springer, 1999. [64] [32] H. Driouch, G. Melzer, and C. Wittmann. Integration of in vivo and in silico metabolic fluxes for improvement of recombinant protein production. Metabolic Engineering, 2011. doi:10.1016/j.ymben.2011.11.002. [31] [33] W. B. Dunn, N. J. Bailey, and H. E. Johnson. Measuring the metabolome: current analytical technologies. Analyst, 130:606–625, 2005. [2] [34] M. Durot, P.-Y. Bourguignon, and V. Schachter. Genome-scale models of bacterial metabolism: reconstruction and applications. FEMS Microbiology Reviews, 33:164–90, 2009. [1, 39] [35] J. S. Edwards, R. U. Ibarra, and B. O. Palsson. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nature Biotechnology, 19:125– 130, 2001. [i, iii, 5, 39] [36] A. M. Feist, C. S. Henry, J. L. Reed, M. Krummenacker, A. R. Joyce, P. D. Karp, L. J. Broadbelt, V. Hatzimanikatis, and B. O. Palsson. A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Molecular Systems Biology, 3(121), 2007. [1, 5, 19] [37] A. M. Feist and B. O. Palsson. The biomass objective function. Current Opinion in Microbiology, 13:344–349, 2010. [2, 39] [38] O. Fiehn. Metabolomics – the link between genotypes and phenotypes. Plant Molecular Biology, 2002. [2] [39] M. Fischetti, D. Salvagnin, and A. Zanette. Minimal infeasible subsystems and Benders cuts, 2008. [56, 57] [40] R. Fleming and I. Thiele. Mass conserved elementary kinetics is sufficient for the existence of a non-equilibrium steady state concentration. Arxiv preprint arXiv:1109.4498, 2011. [78] [41] R. Fleming, I. Thiele, and H. Nasheuer. Quantitative assignment of reaction directionality in constraint-based models of metabolism: Application to Escherichia coli. Biophysical Chemistry, 145:47–56, 2009. [2, 5, 19] [42] R. Fleming, I. Thiele, G. Provan, and H. Nasheuer. Integrated stoichiometric, thermodynamic and kinetic modelling of steady state metabolism. Journal of Theoretical Biology, 264:683–692, 2010. [43, 66] 85

BIBLIOGRAPHY

[43] R. M. Fleming, C. M. Maes, M. A. Saunders, Y. Ye, and B. O. Palsson. A variational principle for computing nonequilibrium fluxes and potentials in genome-scale biochemical networks. Journal of Theoretical Biology, 292:71–77, 2012. [44, 45] [44] R. Fletcher, S. Leyffler, and P. L. Toint. On the global convergence of a filter-SQP algorithm. SIAM Jorunal on Optimization, 2002. [42] [45] J. Flum and M. Grohe. Parametrized Complexity Theory. Texts in Theoretical Computer Science. An EATCS Series. Springer, 2006. [64] [46] S. Fortune, J. Hopcroft, and J. Wyllie. The directed subgraph homeomorphism problem. Theoretical Computer Science, 10(111-121), 1980. [28, 29] [47] C. Francke, R. J. Siezen, and B. Teusink. Reconstructing the metabolic network of a bacterium from its genome. Trends in Microbiology, 13:550–558, 2005. [1] [48] J. Gagneur and S. Klamt. Computation of elementary modes: a unifying framework and the new binary approach. BMC Bioinformatics, 5(175), 2004. [33] [49] R. Ganian, P. Hlinˇen´ y, J. Kneis, A. Langer, J. Obdrˇz´alek, and P. Rossmanith. On digraph width measures in parameterized algorithmics. Parameterized and Exact Computation, 5917/2009:185–197, 2009. [76] [50] A. Gevorgyan, M. G. Poolman, and D. A. Fell. Detection of stoichiometric inconsistencies in biomolecular models. Bioinformatics, 2008. [2] [51] P. E. Gill, W. Murray, and M. A. Saunders. SNOPT: An SQP algorithm for large-scale constrained optimization. SIAM Review, 47(1):99–131, 2005. [42] [52] J. Gleeson and J. Ryan. Identifying minimally infeasible subsystems of inequalities. ORSA Journal on Computing, 2(1):61–63, 1990. [56] [53] M. A. Goberna, V. Jornet, and M. M. Rodr´ıguez. On linear systems containing strict inequalities. Linear Algebra and its Applications, 360:151–171, 2003. [78] [54] M. A. Goberna and M. M. Rodr´ıguez. Analyzing linear systems containing strict inequalities via evenly convex hulls. European Journal of Operational Research, 169:1079–1095, 2006. [78] [55] J. L. Griffin. The Cinderella story of metabolic profiling: does metabolomics get to go to the functional genomics ball? Philosophical Transactions of the Royal Society B, 2006. [2] [56] I. E. Grossmann. Review of nonlinear mixed-integer and disjunctive programming techniques. Optimization and Engineering, pages 227–252, 2002. [42] [57] M. Gr¨ otschel, L. L., and S. A. Geometric Algorithms and Combinatorial Optimization. Springer, 1988. [50] [58] B. Gr¨ unbaum. Convex Polytopes. Graduate Texts in Mathematics. Springer, 2 edition, 2003. [31, 32] [59] S. Gudmundsson and I. Thiele. Computationally efficient flux variability analysis. BMC Bioinformatics, 11(489), 2010. [71] 86

BIBLIOGRAPHY

[60] V. Hatzimanikatis, C. Li, J. A. Ionita, C. S. Henry, M. D. Jankowski, and L. J. Broadbelt. Exploring the diversity of complex metabolic networks. Bioinformatics, (8):1603–1609, 2005. [9] [61] U.-U. Haus, S. Klamt, and T. Stephen. Computing knock-out strategies in metabolic networks. Journal of Computational Biology, 15(3):259–268, 2008. [31] [62] C. S. Henry, L. J. Broadbelt, and V. Hatzimanikatis. Thermodynamics-based metabolic flux analysis. Biophysical Journal, 92:1792–1805, 2007. [45] [63] W. J. Heuett and H. Qian. Combining flux and energy balance analysis to model large-scale biochemical networks. Journal of Bioinformatics and Computational Biology, 4(6):1227– 1243, 2006. [43, 44, 66, 77] [64] A. Hoppe, S. Hoffmann, and H.-G. Holzh¨ utter. Including metabolite concentrations into flux balance analysis: thermodynamic realizability as a constraint on flux distributions in metabolic networks. BMC Systems Biology, 1(23), 2007. [45, 66, 68, 77] [65] P. Hunter and S. Kreutzer. Digraph measures: Kelly decompositions, games, and orderings. Theoretical Computer Science, 399(3):206–219, 2008. [76] [66] L. R. J. Conventional wisdom and P=NP. Internet-Blog, G¨odel’s Lost Letter and P=NP, January 2012. https://rjlipton.wordpress.com/conventional-wisdom-and-pnp/. [76] [67] M. D. Jankowski, C. S. Henry, L. J. Broadbelt, and V. Hatzimanikatis. Group contribution method for thermodynamic analysis of complex metabolic networks. Biophysical Journal, 95:1487–1499, 2008. [19] [68] T. Johnson, N. Robertson, S. P.D., and R. Thomas. Directed tree-width. Journal of Combinatorial Theory, Series B, 82:138–155, 2001. [29, 76] [69] D. E. Kaufman and R. L. Smith. Direction choice for accelerated convergence in hit-and-run sampling. Operations Research, 46(1):84–95, 1998. [30] [70] L. Khachiyan, E. Boros, K. Borys, K. Elbassioni, and V. Gurvich. Generating all vertices of a polyhedron is hard. Discrete Computational Geometry, 39:174–190, 2008. [33] [71] C. Khannapho, H. Zhao, B. L. Bonde, A. M. Kierzek, C. A. Avignone-Rossa, and M. E. Bushell. Selection of objective function in genome scale flux balance analysis for process feed development in antibiotic production. Metabolic Engineering, 10(5):227–233, 2008. [i, iii, 39]

[72] P. Kharchenko, D. Vitkup, and G. M. Church. Filling gaps in a metabolic network using expression information. Bioinformatics, 20:178–185, 2004. [1] [73] S. Klamt and E. D. Gilles. Minimal cut sets in biochemical reaction networks. Bioinformatics, 20(2):226–234, 2004. [39, 58, 59] [74] M. V. Kritz, M. T. dos Santos, S. Urrutia, and J.-M. Schwartz. Organising metabolic networks: Cycles in flux distributions. Journal of Theoretical Biology, 265:250–260, 2010. [17]

[75] H. Kuhn. Solvability and consistency for linear equations and inequalities. The American Mathematical Monthly, 63(4):217–232, 1956. [78] 87

BIBLIOGRAPHY

[76] A. K¨ ummel, S. Panke, and M. Heinemann. Putative regulatory sites unraveled by networkembedded thermodynamic analysis of metabolome data. Molecular Systems Biology, (2006.0034), 2006. [1, 5, 19, 77] [77] A. K¨ ummel, S. Panke, and M. Heinemann. Systematic assignment of thermodynamic constraints in metabolic network models. BMC Bioinformatics, 7(512), 2006. [19] [78] A. Larhlimi and A. Bockmayr. Minimal metabolic behaviors and the reversible metabolic space. In Preprint 299, DFG Research Center Matheon, 2005. [32] [79] M. Las Vergnas. Matro¨ıdes orientables, 1974. papers/1974MO.ps. [11]

www.math.jussieu.fr/~mlv/maths/

[80] W. Liebermeister and E. Klipp. Bringing metabolite networks to life: convenience rate law and thermodynamic constraints. Theoretical Biology and Medical Modelling, 3(41), 2006. [78]

[81] R. Mahadevan and C. Schilling. The effects of alternate optimal solutions in constraintbased genome-scale metabolic models. Metabolic Engineering, 5:264–276, 2003. [i, iii, 39] [82] M. L. Mavrovouniotis. Group contributions for estimating standard Gibbs energies of formation of biochemical compounds in aqueous solution. Biotechnology and Bioengineering, 36:1070–1082, 1990. [19] [83] M. L. Mavrovouniotis. Duality theory for thermodynamic bottlenecks in bioreaction pathways. Chemical Engineering Science, 51(9):1495–1507, 1996. [20, 21] [84] R. Nigam and S. Liang. A pivoting algorithm for metabolic metworks in the presence of themodynamic constraints. In Computational Systems Bioinformatics Conference, 2005. Proceedings. 2005 IEEE, pages 259–267, 2005. [19] [85] R. Nigam and S. Liang. Algorithm for perturbing thermodynamically infeasible metabolic networks. Computers in Biology and Medicine, 37:126–133, 2007. [18, 19] [86] I. E. Nikerel, W. A. van Winden, P. J. Verheijen, and J. J. Heijnen. Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics. Metabolic Engineering, 11:20–30, 2009. [3] [87] R. P. Nolan, A. P. Fenley, and K. Lee. Identification of distributed metabolic objectives in the hypermetabolic liver by flux and energy balance analysis. Metabolic Engineering, 8:30–45, 2006. [23, 77] [88] R. A. Notebaart, F. H. van Enckevort, C. Francke, R. J. Siezen, and B. Teusink. Accelerating the reconstruction of genome-scale metabolic networks. BMC Bioinformatics, 7(296), 2006. [1] [89] J. S. Oliveira, C. G. Bailey, J. B. Jones-Oliveira, and D. A. Dixon. An algebraiccombinatorial model for the identification and mapping of biochemical pathways. Bulletin of Mathematical Biology, 63:1163–1196, 2001. [11, 35] [90] J. S. Oliveira, C. G. Bailey, J. B. Jones-Oliveira, D. A. Dixon, D. W. Gull, and M. L. Chandler. A computational model for the identification of biochemical pathways in the krebs cycle. Journal of Computational Biology, 10(1):57–82, 2003. [35] 88

BIBLIOGRAPHY

[91] J. S. Oliveira, J. B. Jones-Oliveira, D. A. Dixon, C. G. Bailey, and D. W. Gull. Hyperdigraph-theoretic analysis of the EGFR signaling network: Initial steps leading to GTP: Ras complex formation. Journal of Computational Biology, 11(5):812–842, 2004. [33, 35]

[92] M. A. Orman, I. P. Androulakis, F. Berthiaume, and M. G. Ierapetritou. Metabolic network analysis of perfused livers under fed and fasted states: Incorporating thermodynamic and futile-cycle-associated regulatory constraints. Journal of Theoretical Biology, 293:101–110, 2012. [31] [93] J. Orth, T. Conrad, J. Na, J. Lerman, H. Nam, A. Feist, and B. Palsson. A comprehensive genome-scale reconstruction of Escherichia coli metabolism—2011. Molecular Systems Biology, 7(535), 2011. [i] [94] J. D. Orth, I. Thiele, and B. O. Palsson. What is flux balance analysis. Nature Biotechnology, 28:245–248, 2010. [5] [95] G. Oster, A. Perelson, and A. Katchalsky. 234(5329):393–399, 1971. [4, 5]

Network thermodynamics.

Nature,

[96] A. J. Papin, J. Stelling, N. D. Price, S. Klamt, S. Schuster, and B. O. Palsson. Comparison of network-based pathway analysis methods. TRENDS in Biotechnology, 22(8):400–405, 2004. [31, 32] [97] J. L. Papin, Jason A. Reed and B. O. Palsson. Hierarchical thinking in network biology: the unbiased modularization of biochemical networks. TRENDS in Biochemical Sciences, 29(12):641–647. [29] [98] N. D. Price, I. Famili, D. A. Beard, and B. O. Palsson. Extreme pathways and Kirchhoff’s second law. Biophysical Journal, 83:2879–2882, 2002. [18] [99] N. D. Price, J. L. Reed, and B. Ø. Palsson. Genome-scale models of microbial cells: evaluating the consequences of constraints. Nature Reviews Microbiology, 2:886–897, 2004. [3, 4, 10, 39]

[100] N. D. Price, I. Thiele, and B. O. Palsson. Candidate states of Helicopacter pylori’s genomescale metabolic network upon application of ”loop law” thermodynamic constraints. Biophysical Journal, 90:3919–3928, 2006. [5, 19, 26, 30, 41] [101] H. Qian and D. A. Beard. Thermodynamics of stoichiometric biochemical networks in living systems far from equilibrium. Biophysical Chemistry, 114:213–220, 2005. [4, 5, 10] [102] H. Qian and D. A. Beard. Metabolic futile cycles and their functions: A systems analysis of energy and control, June 2006. [17] [103] H. Qian, D. A. Beard, and S.-d. Liang. Stoichiometric network theory for nonequilibrium biochemical systems. European Journal of Biochemistry, 270:415–421, 2003. [4] [104] V. N. Reddy, M. L. Mavrovouniotis, and M. N. Liebman. Petri net representations in metabolic pathways. In Proceedings International Conference on Intelligent Systems for Molecular Biology, volume 1, pages 328–336, 1993. [6] [105] M. A. Safari. D-width: A more natural measure for directed tree width. Mathematical Foundations of Computer Science, 3618/2005:745–756, 2005. [76] 89

BIBLIOGRAPHY

[106] B. Sarıyar, S. Perk, U. Akman, and A. Horta¸csu. Monte carlo sampling and principal component analysis of flux distributions yield topological and modular information on metabolic networks. Journal of Theoretical Biology, 242:389–400, 2006. [29] [107] U. Sauer. Metabolic networks in motion: Biology, 2(62), 2006. [2]

13

C-based flux analysis. Molecular Systems

[108] A. Saxena. On the set-covering polytope: I. all the facets with coefficients in {0, 1, 2, 3}. https://student-3k.tepper.cmu.edu/GSIADOC/wp/2004-E29.pdf, April 2004. [56] [109] J. Schellenberger, N. E. Lewis, and B. Ø. Palsson. Elimination of thermodynamically infeasible loops in steady-state metabolic models. Biophysical Journal, 100:544–553, 2011. [5, 26, 30, 45, 61, 66, 67, 70, 75]

[110] J. Schellenberger and B. O. Palsson. Use of randomized sampling for analysis of metabolic networks. The Journal of Biological Chemistry, 284(9):5457–5461, 2009. [29, 30] [111] J. Schellenberger, J. O. Park, T. M. Conrad, and B. O. Palsson. BiGG: a biochemical genetic and genomic knowledgebase of large scale metabolic reconstructions. BMC Bioinformatics, 11(213), 2010. [65, 70] [112] J. Schellenberger, R. Que, R. M. Fleming, I. Thiele, J. D. Orth, A. M. Feist, D. C. Zielinski, A. Bordbar, N. E. Lewis, S. Rahmanian, J. Kang, D. R. Hyduke, and B. O. Palsson. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA toolbox v2.0. Nature Protocols, 6(9):1290–1307, 2011. [4, 61, 70] [113] C. H. Schilling, D. Letscher, and B. O. Palsson. Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function form a pathwayoriented perspective. Journal of Theoretical Biology, 203:229–248, 2000. [17] [114] A. Schrijver. Combinatorial Optimization, Polyhedra and Efficiency, volume A of Algorithms and Combinatorics. Springer, 2003. [34] [115] S. Schuster, D. A. Fell, and T. Dandekar. A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks. Nature Biotechnology, 18:326–332, 2000. [31] [116] S. Schuster, T. Pfeiffer, and D. A. Fell. Is maximization of molar yield in metabolic networks favoured by evolution? Journal of Theoretical Biology, 252(3):497–504, 2007. [40] [117] J.-M. Schwartz and M. Kanehisa. Quantitative elementary mode analysis of metabolic pathways: the example of yeast glycolysis. BMC Bioinformatics, 7(186), 2006. [31] [118] D. Segre, D. Vitkup, and G. Church. Analysis of optimality in natural and perturbed metabolic networks. Proceedings of the National Academy of Sciences of the United States of America, 2001. [1] [119] Y. Shiloach. A polynomial solution to the undirected two paths problem. Journal of the Association for Computing Machinery, 27(3):445–456, 1980. [29] [120] R. Silva, M. Ulbrich, S. Ulbrich, and V. L.N. A globally convergent primal-dual interiorpoint filter method for nonlinear-programming: new filter optimality measures and computational results. Pr´e-Publica¸c˝oes DMUC. 08-49, 2008. http://hdl.handle.net/10316/ 11218. [42] 90

BIBLIOGRAPHY

[121] E. J. Smid, D. Molenaar, and d. V. W. M. T. B. Hugenholtz, Jeroen. Functional ingredient production: application of global metabolic models. Current Opinion in Biotechnology, 2005. [1] [122] J. Stelling, T. Klamt, K. Bettenbrock, S. Schuster, and E. D. Gilles. Metabolic network structure determines key aspects of functionality and regulation. Nature, 420(190-193), 2002. [1] [123] R. Steuer, T. Gross, J. Selbig, and B. Blasius. Structural kinetic modeling of metabolic networks. Proceedings of the National Academy of Sciences of the United States of America, 2006. [3] [124] M. Terzer. Large Scale Methods to Enumerate Extreme Rays and Elementary Modes. PhD thesis, Swiss Federal Institute of Technology, Zurich, 2009. [31, 33, 36, 37, 41] [125] M. Terzer, N. D. Maynard, M. W. Covert, and J. Stelling. Genome-scale metabolic networks. Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 1(3):285–297, 2009. [1, 4, 10, 39, 65]

[126] B. Teusink, A. Wiersma, L. Jacobs, R. Notebaart, and E. Smid. Understanding the adaptive growth strategy of Lactobacillus plantarum by in silico optimisation. PLoS Computational Biology, 2009. [40] [127] I. Thiele, N. D. Price, T. D. Vo, and B. O. Palsson. Candidate metabolic network states in human mitochondria. The Journal of Biological Chemistry, 280:11683–11695, 2005. [29] [128] G. H. Thomas, J. Zucker, S. Macdonald, A. Sorokin, I. Goryanin, and A. E. Douglas. A fragile metabolic network adapted for cooperation in the symbiotic bacterium Buchnera aphidicola. BMC Systems Biology, 3(24), 2009. [1] [129] S. G. Thorleifsson and I. Thiele. rBioNet: A COBRA toolbox extension for reconstructing high-quality biochemical networks. Bioinformatics, 27(14):2009–2010, 2011. [1] [130] R. J. van Berlo, D. de Ridder, J.-M. Daran, P. A. Daran-Lapujade, B. Teusink, and M. J. Reinders. Predicting metabolic fluxes using gene expression differences as constraints. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2011. [3, 77] [131] R. J. Vanderbei. LOQO: an interior point code for quadratic programming. Optimization Methods and Software, 1999. [42] [132] A. Varma and B. O. Palsson. Metabolic flux balancing: Basic concepts, scientific and practical use. Nature Biotechnology, 12:994–998, 1994. [39] [133] K. Voss, M. Heiner, and I. Koch. Steady state analysis of metabolic pathways using Petri nets. In Silico Biology, 3(31), 2003. [6] [134] W. Wiechert.

13

C metabolic flux analysis. Metabolic Engineering, 3:195–206, 2001.

[2]

[135] L. A. Wolsey. Integer Programming. Wiley-interscience series in discrete mathematics. John Wiley & Sons, 1998. [48, 51, 75] [136] J. Wright and A. Wagner. Exhaustive identification of steady state cycles in large stoichiometric networks. BMC Systems Biology, 2(61), 2008. [36, 65] 91

BIBLIOGRAPHY

[137] M. Xu, R. Smith, and J. Sadhukhan. Optimization of productivity and thermodynamic performance of metabolic pathways. Industrial & Engineering Chemistry Research, 2008. [23]

[138] F. Yang and D. A. Beard. Thermodynamically based profiling of drug metabolism and drugdrug metabolic interactions: A case study of acetaminophen and ethanol toxic interaction. Biophysical Chemistry, 120:121–134, 2006. [5, 19, 77] [139] F. Yang, F. Qi, and D. A. Beard. Directionality is an inherent property of biochemical networks, Feb. 2007. [27] [140] F. Yang, H. Qian, and D. A. Beard. Ab initio predicition of thermodynamically feasible reaction directions from biochemical network stroichiometry. Metabolic Engineering, 7:251– 259, 2005. [27] [141] G. M. Ziegler. Lectures on Polytopes, chapter 6.3, pages 160–162. Springer, 2007. 15, 16]

92

[12, 13,

Declaration I assert to have written this master’s thesis by my own, that I cited also the sources I used, and that I did not use any illegitimate tools. Ich versichere diese Arbeit selbst¨ andig, ohne unerlaubte Hilfsmittel verfasst und keine außer den angef¨ uhrten Quellen verwendet zu haben. Berlin, den 27. Februar 2012

Arne M¨ uller

93