Bioinformatics III: The Double Description Method

V17 The Double Description method:Theoretical framework behind EFM and EP in „Combinatorics and Computer Science Vol. 1120“ edited by Deza, Euler, Manoussakis, Springer, 1996:91 Bioinformatics III

Double Description Method (1953) All known algorithms for computing EMs are variants of the Double Description Method. • derive simple & efficient algorithm for extreme ray enumeration, the so-called Double Description Method. • show that it serves as a framework to the popular EM computation methods. Analogy with Computer Graphics problem: How can one efficiently describe the space in a dark room that is lighted by a torch shining through the open door? Bioinformatics III

Duality of Matrices This is the duality Left: all points above the dividing line (the shaded area) fulfill the condition x 0. Middle: the points in the grey area fulfill the conditions x1 0 and x2 0. But how could we describe the points in the grey area on the right side in a correspondingly simple manner? Obviously, we could define a new coordinate system (r1, r2) as a new set of generating vectors. But we could also try to transform this area back into the grey area of the middle panel and use the old axes x1 and x2. In 2D, this transformation can be obviously best performed by multiplying all vectors inside the grey area by a two-dimensional rotation matrix. Bioinformatics III

The Double Description Method A pair (A,R) of real matrices Aand R is said to be a double description pair or simply a DD pair if the relationship A x 0 if and only if x = R  for some   0 holds. Clearly, for a pair (A,R) to be a DD pair, the column size of Ahas to equal the row size of R, say d. For such a pair, the set P(A) represented by Aas is simultaneously represented by R as A subset P of d is called polyhedral cone if P = P(A) for some matrix A, and Ais called a representation matrix of the polyhedral cone P(A). Then, we say Ris a generating matrix for P. Clearly, each column vector of a generating matrix Rlies in the cone P and every vector in P is a nonnegative combination of some columns of R. Bioinformatics III

The Double Description Method Theorem 1 (Minkowski‘s Theorem for Polyhedral Cones) For any m  n real matrix A, there exists some d  m real matrix Rsuch that (A,R) is a DD pair, or in other words, the cone P(A) is generated by R. The theorem states that every polyhedral cone admits a generating matrix. The nontriviality comes from the fact that the row size of R is finite. If we allow an infinite size, there is a trivial generating matrix consisting of all vectors in the cone. Also the converse is true: Theorem 2 (Weyl‘s Theorem for Polyhedral Cones) For any d  n real matrix R, there exists some m  d real matrix Asuch that (A,R) is a DD pair, or in other words, the set generated by Ris the cone P(A). Bioinformatics III

The Double Description Method Task: how does one construct a matrix Rfrom a given matrix A, and the converse? These two problems are computationally equivalent. Farkas‘ Lemma shows that (A,R) is a DD pair if and only if (RT,AT) is a DD pair. A more appropriate formulation of the problem is to require the minimality of R: find a matrix R such that no proper submatrix is generating P(A). A minimal set of generators is unique up to positive scaling when we assume the regularity condition that the cone is pointed, i.e. the origin is an extreme point of P(A). Geometrically, the columns of a minimal generating matrix are in 1-to-1 correspondence with the extreme rays of P. Thus the problem is also known as the extreme ray enumeration problem. No efficient (polynomial) algorithm is known for the general problem. Bioinformatics III

Double Description Method: primitive form Suppose that the m  d matrix Ais given and let (This is equivalent to the situation at the beginning of constructing EPs or EFMs: we only know S.) The DD method is an incremental algorithm to construct a d  m matrix R such that (A,R) is a DD pair. Let us assume for simplicity that the cone P(A) is pointed. Let Kbe a subset of the row indices {1,2,...,m}of A and let AKdenote the submatrix of A consisting of rows indexed by K. Suppose we already found a generating matrix Rfor AK, or equivalently, (AK,R) is a DD pair. If A = AK ,we are done. Otherwise we select any row index i not in K and try to construct a DD pair (AK+i, R‘) using the information of the DD pair (AK,R). Once this basic procedure is described, we have an algorithm to construct a generating matrix Rfor P(A). Bioinformatics III

Geometric version of iteration step The procedure can be understood geometrically by looking at the cut-section C of the cone P(AK) with some appropriate hyperplane h in d which intersects with every extreme ray of P(AK) at a single point. Let us assume that the cone is pointed and thus C is bounded. Having a generating matrix R means that all extreme rays (i.e. extreme points of the cut-section) of the cone are represented by columns of R. Such a cutsection is illustrated in the Fig. Here, C is the cube abcdefgh. Bioinformatics III

Geometric version of iteration step The newly introduced inequality Aix  0 partitions the space dinto three parts: Hi+ = {x d : Aix> 0 } Hi0 = {x d : Aix = 0 } Hi- = {x d : Aix < 0 } The intersection of Hi0 with P and the new extreme points i and j in the cut-section C are shown in bold in the Fig. Let J be the set of column indices of R. The rays rj (j J ) are then partitioned into three parts accordingly: J+ = {j J : rj Hi+ } J0 = {j J : rj Hi0} J- = {j J : rj Hi-} We call the rays indexed by J+, J0, J-the positive, zero, negative rays with respect to i, respectively. To construct a matrix R‘ from R, we generate new | J+|  | J-| rays lying on the ith hyperplane Hi0 by taking an appropriate positive combination of each positive ray rj and each negative ray rj‘ and by discarding all negative rays. Bioinformatics III

Geometric version of iteration step The following lemma ensures that we have a DD pair (AK+i ,R‘), and provides the key procedure for the most primitive version of the DD method. Lemma 3 Let (AK,R) be a DD pair and let i be a row index of A not in K. Then the pair (AK+i ,R‘) is a DD pair, where R‘ is the d  |J‘ | matrix with column vectors rj (j  J‘) defined by J‘ = J+ J0  (J+  J-), and rjj‘= (Airj)rj‘– (Airj‘)rjfor each (j,j‘) J+  J- Proof omitted. Bioinformatics III

Finding seed DD pair It is quite simple to find a DD pair (AK,R) when |K| = 1, which can serve as the initial DD pair. Another simple (and perhaps the most efficient) way to obtain an initial DD form of P is by selecting a maximal submatrix AK of A consisting of linearly independent rows of A. The vectors rj‘s are obtained by solving the system of equations AK R = I where I is the identity matrix of size |K|, R is a matrix of unknown column vectors rj, j J. As we have assumed rank(A) = d, i.e. R = AK-1 , the pair (AK,R) is clearly a DD pair, since AKx  0  x = AK-1,   0. Bioinformatics III

Primitive algorithm for DoubleDescriptionMethod This algorithm is very primitive, and the straightforward implementation will be quite useless, because the size of J increases extremely fast. This is because many vectors rjj‘ generated by the algorithm (defined in Lemma 3) are unnessary. We need to avoid generating redundant vectors. For this, we will use the zero set or active set Z(x) which is the set of inequality indices satisfied by x in P(A) with equality. Noting Ai• the ith row of A, Z(x) = {i : A i• x = 0} Bioinformatics III

Towards the standard implementation Two distinct extreme rays r and r‘ of P are adjacent if the minimal face of P containing both contains no other extreme rays. Proposition 7. Let r and r‘ be distinct rays of P. Then the following statements are equivalent (a) r and r‘ are adjacent extreme rays, (b) rand r‘ are extreme rays and the rank of the matrix AZ(r)  Z(r‘) is d – 2, (c) if r‘‘is a ray with Z(r‘‘)  Z(r)  Z(r‘) then either r‘‘ ≃ ror r‘‘ ≃r ‘. Lemma 8. Let (AK,R) be a DD pair such that rank(AK) = d and let i be a row index of A not in K. Then the pair (AK+i , R‘) is a DD pair, where R‘ is the d  |J‘| matrix with column vectors rj (j  J‘) defined by J‘ = J+ J0  Adj Adj = {(j,j‘)  J+ J- : rjand rj‘ are adjacent in P(AK)} and r = (Ai rj ) rj‘ – (Airj ) rjfor each (j,j‘) Adj. Furthermore, if R is a minimal generating matrix for P(AK) then R‘ is a minimal generating matrix for P(AK+i). Bioinformatics III

Algorithm for standard form of double description method Hence we can write a straightforward variation of the DD method which produces a minimal generating set for P: DDMethodStandard(A) such that R is minimal Lemma 8 To implement DDMethodStandard, we must check for each pair of extreme rays rand r‘of P(AK) with Air> 0 and Ai r‘< 0 whether they are adjacent in P(AK). Bioinformatics III

Remaining slides will not be used this year. Bioinformatics III

Tools for analyzing network states The two steps that are used to form a solution space — reconstruction and the imposition of governing constraints — are illustrated in the centre of the figure. Several methods are being developed at various laboratories to analyse the solution space. Ci and Cj concentrations of compounds i and j; EP, extreme pathway; vi and vj fluxes through reactions i and j; v1 –v3 flux through reactions 1-3; vnet, net flux through loop. Price et al. Nature Rev Microbiol 2, 886 (2004) Bioinformatics III

Application of elementary modesMetabolic network structure of E.coli determineskey aspects of functionality and regulation Compute EFMs for central metabolism of E.coli. Catabolic part: substrate uptake reactions, glycolysis, pentose phosphate pathway, TCA cycle, excretion of by-products (acetate, formate, lactate, ethanol) Anabolic part: conversions of precursors into building blocks like amino acids, to macromolecules, and to biomass. Stelling et al. Nature 420, 190 (2002) Bioinformatics III

Metabolic network topology and phenotype Idea: Can the total number of EFMs for given conditions be used as quantitative measure of metabolic flexibility? a, Relative number of EFMs N enabling deletion mutants of gene i (i) inE. coli to grow (abbreviated by µ) for 90 different combinations of mutation and carbon source. Shown are results for 90 deletions of different individual genes. Stelling et al. Nature 420, 190 (2002) Answer: Yes, the # of EFMs for mutant strain allows correct prediction of growth phenotype in more than 90% of the cases. Bioinformatics III

Robustness analysis The # of EFMs qualitatively indicates whether a mutant is viable or not, but does not describe quantitatively how well a mutant grows. Define maximal biomass yield Ymass as the optimum of: eiis the single reaction rate (growth and substrate uptake) in EFM i selected for utilization of substrate Sk. Stelling et al. Nature 420, 190 (2002) Bioinformatics III

Robustness Analysis Dependency of the mutants' maximal growth yield Ymax(Δi) (open circles) and the network diameter D(Δi) (open squares) on the share of elementary modes operational in the mutants. Stelling et al. Nature 420, 190 (2002) Central metabolism of E.coli behaves in a highly robust manner because mutants with significantly reduced metabolic flexibility show a growth yield similar to wild type. 17. Lecture WS 2013/14 Bioinformatics III

Current metabolomics Review: (1) recent work on metabolic networks required revising the picture of separate biochemical pathways into a densely-woven metabolic network (2) Connectivity of substrates in this network follows a power-law (Yeong&Barabasi). (3) Constraint-based modeling approaches (FBA) were successful in analyzing the capabilities of cellular metabolism including - its capacity to predict deletion phenotypes - the ability to calculate the relative flux values of metabolic reactions, and - the capability to identify properties of alternate optimal growth states in a wide range of simulated environmental conditions Open questions - what parts of metabolism are involved in adaptation to environmental conditions? - is there a central essential metabolic core? - what role does transcriptional regulation play? Bioinformatics III

Distribution of fluxes in E.coli Aim: understand principles that govern the use of individual reactions under different growth conditions. Nature 427, 839 (2004) Stoichiometric matrix for E.coli strain MG1655 containing 537 metabolites and 739 reactions taken from Palsson et al. Apply flux balance analysis to characterize solution space (all possible flux states under a given condition). vjis the flux of reaction j and Sijis the stoichiometric coefficient of reaction j. Bioinformatics III

Optimal states Denote the mass carried by reaction j producing (consuming) metabolite i by Observation: Fluxes vary widely: e.g. dimensionless flux of succinyl coenzyme A synthetase reaction is 0.185, whereas the flux of the aspartate oxidase reaction is 10.000 times smaller, 2.2  10-5. Using linear programming and adapting constraints for each reaction flux viof the form imin≤ vi ≤ imax, the flux states were calculated that optimize cell growth on various substrates. Plot the flux distribution for active (non-zero flux) reactions of E.coli grown in a glutamate- or succinate-rich substrate. Bioinformatics III

Overall flux organization of E.coli metabolic network a, Flux distribution for optimized biomass production on succinate (black) and glutamate (red) substrates. The solid line corresponds to the power-law fit that a reaction has flux v P(v)  (v + v0)- , with v0 = 0.0003 and  = 1.5. d, The distribution of experimentally determined fluxes from the central metabolism of E. coli shows power-law behaviour as well, with a best fit to P(v) v- with  = 1. Both computed and experimental flux distribution show wide spectrum of fluxes. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Response to different environmental conditions Is the flux distribution independent of environmental conditions? b, Flux distribution for optimized biomass on succinate substrate (black) with an additional 10% (red), 50% (green) and 80% (blue) randomly chosen subsets of the 96 input channels (substrates) turned on. The flux distribution was averaged over 5,000 independent random choices of uptake metabolites.  Yes, the flux distribution is independent of the external conditions. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Use scaling behavior to determine local connectivity The observed flux distribution is compatible with two different potential local flux structures: (a) a homogenous local organization would imply that all reactions producing (consuming) a given metabolite have comparable fluxes (b) a more delocalized „high-flux backbone (HFB)“ is expected if the local flux organisation is heterogenous such that each metabolite has a dominant source (consuming) reaction. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Characterizing the local inhomogeneity of the flux net a, Measured kY(k) shown as a function of k for incoming and outgoing reactions, averaged over all metabolites, indicates that k×Y(k) k0.73. Inset shows non-zero mass flows, v^ij, producing (consuming) FAD on a glutamate-rich substrate.  an intermediate behavior is found between the two extreme cases.  the large-scale inhomogeneity observed in the overall flux distribution is also increasingly valid at the level of the individual metabolites. The more reactions that consume (produce) a given metabolite, the more likely it is that a single reaction carries most of the flux, see FAD. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Clean up metabolic network Use simple algorithm that removes for each metabolite systematically all reactions but the one providing the largest incoming (outgoing) flux distribution. The algorithm uncovers the „high-flux-backbone“ of the metabolism, a distinct structure of linked reactions that form a giant component with a star-like topology. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Maximal flow networks glutamate rich succinate rich substrates Directed links: Two metabolites (e.g. A and B) are connected with a directed link pointing from A to B only if the reaction with maximal flux consuming A is the reaction with maximal flux producing B. Shown are all metabolites that have at least one neighbour after completing this procedure. The background colours denote different known biochemical pathways. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

FBA-optimized network on glutamate-rich substrate High-flux backbone for FBA-optimized metabolic network of E. coli on a glutamate-rich substrate. Metabolites (vertices) coloured blue have at least one neighbour in common in glutamate- and succinate-rich substrates, and those coloured red have none. Reactions (lines) are coloured blue if they are identical in glutamate- and succinate-rich substrates, green if a different reaction connects the same neighbour pair, and red if this is a new neighbour pair. Black dotted lines indicate where the disconnected pathways, for example, folate biosynthesis, would connect to the cluster through a link that is not part of the HFB. Thus, the red nodes and links highlight the predicted changes in the HFB when shifting E. coli from glutamate- to succinate-rich media. Dashed lines indicate links to the biomass growth reaction. (1) Pentose Phospate (11) Respiration (2) Purine Biosynthesis (12) Glutamate Biosynthesis (20) Histidine Biosynthesis (3) Aromatic Amino Acids (13) NAD Biosynthesis (21) Pyrimidine Biosynthesis (4) Folate Biosynthesis (14) Threonine, Lysine and Methionine Biosynthesis (5) Serine Biosynthesis (15) Branched Chain Amino Acid Biosynthesis (6) Cysteine Biosynthesis (16) Spermidine Biosynthesis (22) Membrane Lipid Biosynthesis (7) Riboflavin Biosynthesis (17) Salvage Pathways (23) Arginine Biosynthesis (8) Vitamin B6 Biosynthesis (18) Murein Biosynthesis (24) Pyruvate Metabolism (9) Coenzyme A Biosynthesis (19) Cell Envelope Biosynthesis (25) Glycolysis (10) TCA Cycle Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Interpretation Only a few pathways appear disconnected indicating that although these pathways are part of the HFB, their end product is only the second-most important source for another HFB metabolite. Groups of individual HFB reactions largely overlap with traditional biochemical partitioning of cellular metabolism. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

How sensitive is the HFB to changes in the environment? b, Fluxes of individual reactions for glutamate-rich and succinate-rich conditions. Reactions with negligible flux changes follow the diagonal (solid line). Some reactions are turned off in only one of the conditions (shown close to the coordinate axes). Reactions belonging to the HFB are indicated by black squares, the rest are indicated by blue dots. Reactions in which the direction of the flux is reversed are coloured green. Only reactions in the high-flux territory undergo noticeable differences! Type I: reactions turned on in one conditions and off in the other (symbols). Type II: reactions remain active but show an orders-in-magnitude shift in flux under the two different growth conditions. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Flux distributions for individual reactions Shown is the flux distribution for four selected E. coli reactions in a 50% random environment. a Triosphosphate isomerase; b carbon dioxide transport; c NAD kinase; d guanosine kinase. Reactions on the   vcurve (small fluxes) have unimodal/gaussian distributions (a and c). Shifts in growth-conditions only lead to small changes of their flux values. Reactions off this curve have multimodal distributions (b and d), showing several discrete flux values under diverse conditions. Under different growth conditions they show several discrete and distinct flux values. Almaar et al., Nature 427, 839 (2004) Bioinformatics III

Summary Metabolic network use is highly uneven (power-law distribution) at the global level and at the level of the individual metabolites. Whereas most metabolic reactions have low fluxes, the overall activity of the metabolism is dominated by several reactions with very high fluxes. E. coli responds to changes in growth conditions by reorganizing the rates of selected fluxes predominantly within this high-flux backbone. Apart from minor changes, the use of the other pathways remains unaltered. These reorganizations result in large, discrete changes in the fluxes of the HFB reactions. Bioinformatics III

The same authors as before used Flux Balance Analysis to examine utilization and relative flux rate of each metabolite in a wide range of simulated environmental conditions for E.coli, H. pylori and S. cerevisae: For each system they considered 30.000 randomly chosen combinations where each uptake reaction is a assigned a random value between 0 and 20 mmol/g/h.  adaptation to different conditions occurs by 2 mechanisms: (a) flux plasticity: changes in the fluxes of already active reactions. E.g. changing from glucose- to succinate-rich conditions alters the flux of 264 E.coli reactions by more than 20% (b) less commonly, adaptation includes structural plasticity, turning on previously zero-flux reactions or switching off active pathways. Bioinformatics III

Emergence of the Metabolic Core The two adaptation method mechanisms allow for the possibility of a group of reactions not subject to structural plasticity being active under all environmental conditions. Assume that active reactions were randomly distributed. If typically a q fraction of the metabolic reactions are active under a specific growth condition, we expect for n distinct conditions an overlap of at least qn reactions. This converges quickly to 0. Bioinformatics III

Emergence of the Metabolic Core (a–c) The average relative size of the number of reactions that are always active as a function of the number of sampled conditions (black line). (d and e) The number of metabolic reactions (d) and the number of metabolic core reactions (e) in the three studied organisms. In a-c, as the number of conditions increases, the curve converges to a constant enoted by the dashed line, identifying the metabolic core of an organism. Red line : number of reactions that are always active if activity is randomly distributed in the metabolic network. The fact that it converges to zero indicates that the real core represents a collective network effect, forcing a group of reactions to be active in all conditions. Bioinformatics III

Metabolic Core of E.coli: The constantly active reactions form a tightly connected cluster! Shown are all reactions that are found to be active in each of the 30,000 investigated external conditions. Blue: Metabolites that contribute directly to biomass formation, Red (green): core reactions (links) catalyzed by essential (or nonessential) enzymes. Black-colored links: enzymes with unknown deletion phenotype. Blue dashed lines: multiple appearances of a metabolite, Links with arrows: unidirectional reactions. Note that 20 out of the 51 metabolites necessary for biomass synthesis are not present in the core, indicating that they are produced (or consumed) in a growth-condition-specific manner. Blue and brown shading: folate and peptidoglycan biosynthesis pathways White numbered arrows denote current antibiotic targets inhibited by: (1) sulfonamides, (2) trimethoprim, (3) cycloserine, and (4) fosfomycin. A few reactions appear disconnected since we have omitted the drawing of cofactors. Bioinformatics III

Metabolic Core Reactions The metabolic cores contain 2 types of reactions: (a) reactions that are essential for biomass production under all environment conditions (81 of 90 in E.coli) (b) reactions that assure optimal metabolic performance. Bioinformatics III

Characterizing the Metabolic Cores (A) The number of overlapping metabolic reactions in the metabolic core of H. pylori, E. coli, and S. cerevisiae. The metabolic cores of simple organisms (H. pylori and E.coli) overlap to a large extent. The largest organism (S.cerevisae) has a much larger reaction network that allows more flexbility  the relative size of the metabolic core is much lower. (B) The fraction of metabolic reactions catalyzed by essential enzymes in the cores (black) and outside the core in E. coli and S. cerevisiae.  Reactions of the metabolic core are mostly essential ones. (C) One could assume that the core represents a subset of high-flux reactions. This is apparently not the case. The distributions of average metabolic fluxes for the core and the noncore reactions in E. coli are very similar. Bioinformatics III

Summary - Adaptation to environmental conditions occurs via structural plasticity and/or flux plasticity. Here: identification of a surprisingly stable metabolic core of reactions that are tightly connected to eachother. - the reactions belonging to this core represent potential targets for antimicrobial intervention. Bioinformatics III

Bioinformatics III: The Double Description Method