Global reconstruction of the human metabolic network based on genomic and bibliomic data

Reconstruction and Validation of H. sapiens Recon 1.

A well annotated genome sequence is vital for bottom-up reconstruction because it enables the rapid identification of candidate network components (4) and the assembly of a preliminary network (13) that can be used as a starting point for manual curation [supporting information (SI) Fig. 5]. We used Enzyme Commission numbers (14) and Gene Ontologies (15) to identify an initial set of 1,865 human metabolic genes from the November 2004 annotations (Build 35) of Kyoto Encyclopedia of Genes and Genomes (KEGG) (16), National Center for Biotechnology Information's LocusLink (17) [now EntrezGene (18)], and the H-Invitational Database (19). These genes were mapped to a rudimentary network of 3,623 metabolic enzymes and 3,673 reactions from KEGG's LIGAND database and the compartmentalized yeast metabolic reconstruction (20). In addition to establishing initial network scope, LIGAND's pathway-based organizational structure also facilitated parallel network assembly. A team of researchers simultaneously curated network components by evaluation of >50 years of biological evidence from >1,500 primary literature articles, reviews, and biochemical textbooks. Strict quality control/quality assurance methods were used throughout the reconstruction (see Materials and Methods). Manual literature-based reconstruction ensured that the network components and their interactions were based on direct physical evidence and reflected the current knowledge of human metabolism.

Bottom-up reconstructions can be represented mathematically, enabling the computational interrogation of their properties (4, 21). We validated the basic functionality of the human metabolic network by simulating 288 known metabolic functions in silico (SI Tables 2–4). Like genome sequence assembly and annotation, network reconstruction is an iterative process, and consequently several rounds of iterative gap analysis (i.e., targeted literature searches to identify metabolic reactions that were missed initially) and comprehensive revalidation were required to achieve a BiGG database and a high-quality network reconstruction (SI Fig. 6).

The result of five iterative rounds of reconstruction and validation is H. sapiens Recon 1 and the first human cellular process to be comprehensively modeled at this level of detail and accuracy (Table 1). It was almost entirely constructed from human-specific data and includes many reactions directly extracted from the literature that are not described in any chart or database. Furthermore, Recon 1 represents several hierarchical levels of detail, namely:

• Carefully formulated metabolites and reactions, which account for known reaction stoichiometry, substrate/cofactor specificity, and directionality, as well as overall conservation of mass and charge-based metabolite ionization states at pH 7.2.

• Full compartmentalization of metabolites in and their exchange between seven intracellular locations (cytoplasm, mitochondria, nucleus, endoplasmic reticulum, Golgi apparatus, lysosome, and peroxisome) and the extracellular environment.

• Precise Boolean descriptions of gene–protein relationships such as alternatively spliced variants, protein complexes, and isozymes.

• Confidence scores and literature references based on known biological evidence associated with each gene, protein, and reaction.

The entire contents of H. sapiens Recon 1 is freely available in several formats [searchable database, metabolite and reaction lists, human-specific metabolic maps, stoichiometric matrix, and Systems Biology Markup Language (22)] at http://bigg.ucsd.edu and in SI Figs. 7–16 and SI Tables 5–10. We now describe three applications of Recon 1.

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.

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.

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).

Integrated Analysis of Gene Expression Data Reveals the Effects of Gastric Bypass Surgery on Skeletal Muscle Metabolism.

Recon 1 was used as a context for interpreting the effects of gastric bypass surgery on skeletal muscle metabolism. In this study, gene expression data were acquired from the same patients before and after bariatric surgery. Consequently, many genetic factors that might contribute to interindividual differences are not a concern in this analysis. Published gene expression measurements (43) were mapped to the reconstructed network by using gene–transcript–protein–reaction associations (see Materials and Methods) and visualized on maps of central metabolism and the electron transport chain (Fig. 4; SI Fig. 21). We observed a general trend of up-regulated anaerobic metabolism and down-regulated oxidative phosphorylation after surgery, with many genes in glycolysis, pentose phosphate pathway, methylglyoxal metabolism, and oxidative phosphorylation showing subtle but consistent overall patterns of expression change (SI Tables 11–15). The relative decrease in mitochondrial bioenergetics is also evident in terms of the smaller number of mitochondrial reactions in the down-regulated network (SI Fig. 22) and is consistent with gene expression changes observed in the skeletal muscle of rhesus monkeys subjected to long-term caloric restriction (44). Comparison of gene expression data in the context of Recon 1 suggests that 1 year after surgery, patients may still be feeling the effects of calorie restriction even after weight stabilization. Thus, Recon 1 represents a versatile and effective integration tool, enabling visualization and analysis of genome-scale context in the context of a highly curated metabolic network.

Reconstruction Procedure.

An initial component list was assembled as described in the text. This list was then divided into eight metabolic subsets (amino acids, carbohydrates, energy, glycans, lipids, nucleotides, secondary metabolites/xenobiotics, vitamins, and cofactors) for independent curation by a team of researchers. Putative gene assignments were verified based on evidence collected from genome annotation databases, namely EntrezGene (18), Gene Cards (46), and the scientific literature. Alternative transcripts were identified based on known RefSeq (17) mRNA transcripts for each locus. Substrate and cofactor preferences were identified from the literature and BRENDA (47). Metabolite formula and charge were calculated based on their ionization state at pH 7.2, which for simplicity was presumed to be constant across all compartments. Reaction directionality was determined from thermodynamic data or inferred from legacy data and textbooks. Compartmentalization was determined from protein localization data, sequence targeting signals, and indirect physiological evidence. If these data were unavailable, reactions were modeled as cytoplasmic. The intermembrane space of double-membrane organelles was also modeled as cytoplasmic. Gene–transcript–protein–reaction relationships (5, 6) were manually identified from the literature and formulated as Boolean logic statements. Isozymes (an “or” relationship) were defined as distinct proteins that catalyze the same substrate- and compartment-specific reaction and could arise from one gene due to alternative splicing or could be encoded by independent genes. Cases in which a reaction depended on the presence of more than one gene/protein (an “and” relationship, e.g., proteins with multiple subunits/chains or complexes composed of multiple enzymes) were classified as protein complexes. Confidence scores were assigned based on biological evidence associated with each reaction. Evidence from classical biochemical or genetic experiments, such as gene cloning and protein characterization, was given the highest confidence score (3). Midlevel scores (2) were assigned to reactions based on physiological data or biochemical/genetic evidence from a nonhuman mammalian cell (typically mouse, rat, or rabbit). Reactions with the lowest confidence score (1) were included solely based on in silico modeling because, during the process of model validation, they were deemed mandatory for a particular metabolic function. Transport reactions were entirely reconstructed based on literature reports and biochemistry textbooks because the current annotation of transporters is not sufficiently specific with regard to substrates and mechanisms.

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.

Gene Expression Analysis.

Gene expression data for the gastric bypass study were downloaded from Gene Expression Omnibus (GEO) [GSE5109 (49)]. Expression signals were normalized by the average value for each chip. Log10 ratios of post/pregastric surgery probe expression signals were calculated for all three patients in the gastric bypass data set and mapped to Entrez Gene and RefSeq mRNA identifications (18) in H. sapiens Recon 1 based on database identifiers in the Affymetrix (Santa Clara, CA) U133A Plus 2.0 annotation file (50). Probes whose expression ratio was qualitatively inconsistent across all three patients (i.e., not all up or all down) or that conflicted at the gene level were removed from the data set. Average gene expression ratios were then mapped to the reaction network by using Recon 1's gene–transcript–protein–reaction associations.