Quantitative Characterization of the Human Metabolic Bibliome Reveals an Uneven Knowledge Landscape.

Bottom-up reconstruction of Recon 1 required extensive manual surveys of the primary literature to evaluate biological evidence associated with each gene, protein, and reaction. Viewing confidence scores for these individual components at the system level reveals a global knowledge landscape with specific “peaks” and “valleys” in our understanding of human metabolism (Fig. 1). Three categories of metabolic pathways were identified based on the degree of characterization of their corresponding reactions.

Open in a separate window

Fig. 1.

Human metabolic knowledge landscape. Colors represent the percentage of reactions within a pathway that have a confidence score of 3 (biochemical or genetic evidence), 2 (physiological data or evidence from a nonhuman mammalian cell), 1 (modeling evidence), or 0 (unevaluated). Metabolic pathways (primarily defined by the Kyoto Encyclopedia of Genes and Genomes LIGAND database) were classified into three categories based on their level of characterization as detailed in the text.

Category I pathways are those where extensive primary literature is available. Chondroitin sulfate is a common component of proteoglycans, which are important in cell adhesion, signaling, and connective tissue composition. Catabolism of chondroitin sulfate chains (SI Fig. 15) is a typical example of a Category I pathway in which nearly all of the enzymes have been biochemically characterized and their corresponding genes have been identified. However, the initial steps of chondroitin sulfate catabolism in endosomes (CSBPASEly) and final degradation of the core tetrasaccharide linkage (LINKDEG2ly) are not well known (23).

Category II pathways, such as glyoxylate metabolism (SI Fig. 16), have a roughly equal proportion of highly characterized enzymes and those with moderate biological evidence (see Materials and Methods). For instance, although the peroxisomal and mitochondrial degradation of glyoxylate to l-glycine (reactions AGTix and AGTim, respectively) has been extensively studied, the presence of the glycerate kinase reaction (GLYCK2) was inferred based on the observation that individuals with d-glycericaciduria (who lack the enzyme catalyzing this reaction) cannot further metabolize d-glycerate and excrete gram amounts of it in their urine (24).

Category III pathways exhibit a wide range of confidence scores and gene coverage. That some of these pathways have not been completely elucidated is somewhat surprising, and arguably these knowledge deficits may not have been comprehensively identified without a systems approach. For example, the mechanism that cycles the end products of vitamin C degradation back to the glycolytic pathway appear to be poorly understood (SI Fig. 17) despite evidence in human erythrocytes that it may be used as an energy source (25). Furthermore, although most of ubiquinone 10 biosynthesis was inferred from physiological evidence, there are a few well studied enzymes interspersed in the pathway (SI Fig. 18). A large number of intracellular transport reactions are also included in this category, indicating that as a whole they require considerably more investigation to elucidate precise mechanistic reactions. Thus, the reconstruction of H. sapiens Recon 1 has resulted in a comprehensive appraisal of our knowledge of human metabolism and has led to direct suggestions where further experimental studies are needed (SI Tables 9 and 10).

Singular Value Decomposition (SVD) of the Stoichiometric Matrix Highlights the Importance of Compartmentalization.

Network capabilities are constrained in part by the overall structure of the stoichiometric matrix (S), which can be analyzed by using computational approaches (21). SVD (26) can identify the most influential components of a network. We used SVD to calculate the effective dimensionality of the human network to assess the number of network components needed to account for a given percentage of its structure. The cumulative normalized singular value spectrum is shown in Fig. 2 and compared with that of previous reconstructions (20, 27). The human and yeast networks are considered with and without metabolite compartmentalization to evaluate its effect on network complexity. We observed that the compartmentalized networks have a significantly larger effective dimensionality, requiring a larger number of independent modes to fully describe their contents. For a functional metabolic model, this result may be interpreted as the expansion of nonredundant metabolic functionality rather than simply a division of linear metabolic pathways across various cellular compartments.

Open in a separate window

Fig. 2.

Normalized cumulative singular value spectra for H. sapiens, S. cerevisiae, and E. coli and dominant metabolite modes. (A) Compartmentalized networks have a greater effective dimensionality than their noncompartmentalized counterparts, requiring a larger number of singular values to completely reconstruct the network. Each spectrum shows the number of decomposed modes (x axis) required to reconstruct a given fraction (y axis) of the S matrix's content. (B) The first five modes of the human metabolite coupling matrix (38) highlight the importance of the production and exchange of energy equivalents and the potentially significant impact of osmotic regulation. c, cytoplasmic; e, extracellular; g, Golgi apparatus; m, mitochondrial.

SVD was further applied to demonstrate unbiased systemic links between the metabolites of Recon 1. To minimize any bias arising from differing stoichiometric coefficients and lumped reactions, the metabolite-coupling matrix M (28) was constructed and decomposed. The five most dominant modes of the SVD of M (Fig. 2) are interpreted as independent groups of interacting metabolites, ordered with monotonically decreasing importance. These have biologically meaningful interpretations as currency exchanges, including high-energy phosphate group transfer (modes 1 and 4), reducing equivalent exchange (mode 2), and sugar transfer (mode 3). Mode 5 is dominated by ions and water, linking biochemical and mechanical relationships in the cell by osmotic force balance. Interestingly, these modes are compartment-specific, highlighting known relationships between compartments and cellular functions. Collectively, these observations highlight some of the implications of compartmentalization in the human metabolic reconstruction and support the notion that compartments may function as independent reaction sets to achieve specific metabolic objectives (29). This finding reinforces the idea that intracellular compartmentalization has functional implications for various metabolites, including glutathione (30), amino acids (31–33), cholesterol and bile acids (34, 35), and sphingolipids (36).

Coupled Reaction Sets Suggest Potential Alternatives to Known Drug Targets.

Recon 1 has also enabled in silico characterization of the known human metabolic map by using well established constraint-based methods that have been applied extensively to microbial metabolism (37). Coupled reaction sets consist of reactions that are active together in functional states of a network (29, 38). Thus, a flux through one reaction results in a directionally equivalent flux in other reactions in the set. Flux coupling analysis (39) was used to identify coupled reaction sets under aerobic glucose metabolic conditions. One of the largest of the >250 coupled reaction sets identified under these metabolic conditions involves two branches of cholesterol biosynthesis (Fig. 3A and SI Fig. 19). 3-Hydroxy-3-methylglutaryl-CoA reductase (Entrez Gene ID 3156), a primary metabolic target of the antilipidemic class of statin drugs, is in this coupled reaction set; other members of the set are thus identified as potential alternative drug targets for treating hyperlipidemia. It has been proposed that deficiencies in enzymes belonging to the same functionally coupled reaction set may have similar phenotypes, and examples of this have been reported for the human mitochondria (38). This hypothesis is supported by the example in Fig. 3B (SI Fig. 20), which depicts a coupled reaction set involved in the production and transport of glutathione. The Online Mendelian Inheritance in Man (OMIM) (40, 41) identification tags associated with the genes encoding glutathione synthetase (MIM no. 231900; Entrez Gene ID, 2937) and glutamate-cysteine ligase (OMIM no. 230450; Entrez Gene ID, 2729) indicate that both deficiencies result in hemolytic anemia. These examples demonstrate that in silico experiments with these models provide an analytical approach to studying the causes and consequences of disease states, which can potentially lead to insights into new drug treatment targets. The use of functionally grouped reactions, such as coupled (39) and correlated reaction sets (42), present a promising approach for the functional analysis of complex networks with applications in elucidating causal relationships in the diseases and potentially identifying new treatment strategies and drug targets (29).

Open in a separate window

Fig. 3.

Coupled reaction sets involving cholesterol biosynthesis and glutathione production and transport. (A) The cholesterol biosynthesis coupled set includes all reactions except those shaded in gray. Note that the groups of reactions are not all directly connected and could not be identified by visual inspection alone. (B) Reactions in the glutathione reaction set were mapped to disease associations by using Mendelian Inheritance in Man identification tags. Deficiencies in glutathione synthetase (GTHS) or glutamate-cysteine ligase (GLUCYS) both result in hemolytic anemia, supporting the notion that enzyme deficiencies in the same coupled set may have similar phenotypes. Interference with or decreases in GLUCYS activity is associated with an increased risk of myocardial infarctions (MI). A high-resolution version of this figure is available in SI Figs. 19 and 20.

Functional Validation and Gap Analysis.

The reconstruction was assembled in SimPheny (Genomatica, San Diego, CA), and the stoichiometric matrix was formulated as described (48). Exchange reactions (SI Table 6) were added to enable uptake and secretion of extracellular metabolites for the purpose of simulations. Functional validation was performed by using flux balance analysis (39), allowing recycled cofactor pairs to enter and leave the system as needed (SI Tables 2–4). Comprehensive gap analysis of the stoichiometric matrix was performed after each round of functional validation. Every “dead-end” metabolite that could not be produced or consumed was manually reexamined by returning to the literature to identify possible reactions describing its degradation, production, or transport. A final round of gap analysis was performed upon completion of H. sapiens Recon 1, and a description of unresolved gaps is provided in SI Table 9.

Network Analysis.

The singular value spectra (26) were computed for H. sapiens Recon 1, Saccharomyces cerevisiae iND750 (20), and Escherichia coli iJR904 (48) as the normalized cumulative sum of the singular values by using Matlab Ver. 6.5 (MathWorks, Natick, MA). The stoichiometric matrix was not altered or scaled before SVD. The human and yeast networks were decompartmentalized by reassigning all intracellular metabolites to the cytoplasm. Modes were defined as the largest elements (each >0.25×max column value) of the columns of U (identical to the columns of V), where UΣV^T = ^S^S^T and Ŝ is the binary form of S. Noncompartmentalized versions of the human and yeast model were created by summing up the rows of the S matrix corresponding to the same metabolite in the different intracellular compartments.

Coupled reaction sets were calculated with a bilinear optimization algorithm described by Burgard et al. (39) by using LINDO r6.1 (Lindo Systems, Chicago, IL) called from within Matlab Ver. 6.5 for aerobic glucose uptake conditions. Briefly, for all of the fluxes in the network under specified input conditions, if two fluxes always respond in the same manner (both increase or both decrease) to any perturbation or alternative flux distribution, they are said to be coupled.

We thank Dr. Markus Herrgård for performing the flux coupling computations, Mr. Jan Schellenberger for developing the BiGG Database, and Dr. Evelyn Travnik and Mr. Thomas Conrad for assistance in preparing the supporting information. This work was supported by Training Grants from the National Institutes of Health.