HepatoNet1: a comprehensive metabolic reconstruction of the human hepatocyte for the analysis of liver physiology.

We present HepatoNet1, the first reconstruction of a comprehensive metabolic network of the human hepatocyte that is shown to accomplish a large canon of known metabolic liver functions. The network comprises 777 metabolites in six intracellular and two extracellular compartments and 2539 reactions, including 1466 transport reactions. It is based on the manual evaluation of >1500 original scientific research publications to warrant a high-quality evidence-based model. The final network is the result of an iterative process of data compilation and rigorous computational testing of network functionality by means of constraint-based modeling techniques. Taking the hepatic detoxification of ammonia as an example, we show how the availability of nutrients and oxygen may modulate the interplay of various metabolic pathways to allow an efficient response of the liver to perturbations of the homeostasis of blood compounds.

Introduction

The liver has a wide range of physiological functions, including detoxification of endo- and xenobiotic compounds, homeostatic regulation of the plasma concentration of a multitude of metabolites, synthesis of most plasma proteins, bile formation and hormone production. Hepatocytes account for about 60% of the liver in terms of cell number and for 80% of liver volume (Sasse et al, 1992). They are the principal site of the metabolic conversions underlying the diverse physiological functions of the liver. Accordingly, hepatocyte metabolism integrates a vast array of differentially regulated biochemical pathways and is highly responsive to changes in blood composition (Dardevet et al, 2006) evoked by special diets, overnutrition, starvation or enhanced physical activity. Overload of the liver's metabolic capacity as in drug abuse or inappropriate diets may give rise to long-term silent disease progression such as accumulation of triglycerides in the liver (non-alcoholic fatty liver), changes in the plasma level of lipoproteins, enhanced intracellular levels of reactive oxygen species and changes in the composition of bile. The complexity of the metabolic network and the numerous physiological functions to be achieved within a highly variable physiological environment necessitate an integrated approach with the aim of understanding liver metabolism at a systems level.

Over the past decades, a remarkable amount of biochemical, physiological and medical data on enzymes, membrane transporters, metabolic pathways and physiological functions of hepatocytes from human beings and other mammals have been compiled. In particular, two global reconstructions of the human metabolic network have been published (Duarte et al, 2007; Ma et al, 2007) that represent an excellent starting point for the reconstruction of tissue-specific metabolic networks. With the aim to predict metabolic networks of various human tissues solely based on large-scale transcriptional profiles, a computational method has been proposed (Shlomi et al, 2008), which relates the abundance of the protein transcript to the likelihood that an enzyme occurring in the global reconstruction of the human metabolic network carries metabolic flux. However, when testing the functional capacity of the predicted liver network to recapitulate typical metabolic objectives reported for the human hepatocyte, we failed in a substantial number of cases tested (see below). This negative finding underlines the need to choose the labor-intensive strategy to perform a thorough manual curation of the liver network in order to complement information on transcript and proteins levels of enzymes and membrane transporters with additional biochemical and physiological evidences available in the scientific literature and public databases. The curation process was tightly linked to rigorous computational testing of the functional capacity of the network by means of constraint-based optimization methods. To this end, we defined a large canon of known metabolic objectives of human hepatocytes ranging from the production of glucose from lactate (gluconeogenesis) to the degradation of endocytosed plasma proteins. We also performed negative tests by checking the incapability of the network to allow thermodynamically infeasible metabolic processes (e.g. the formation of energy-rich phosphate bonds without supply of external substrates) or metabolic conversions that have been shown not to exist in mammalian cells (e.g. the conversion of fatty acids with an even number of carbon atoms into glucose). As the result of this reconstruction process, we present HepatoNet1, a genome-scale metabolic network of human hepatocytes, which enables the application of constraint-based modeling techniques to discern allowable metabolic states in hepatocytes at a large variety of physiological conditions.

Results

Curation of the stoichiometric metabolic network of human hepatocytes

The primary aim of our work was to establish a stoichiometric model of human hepatocyte metabolism characterized by (i) comprehensive coverage of known biochemical activities of hepatocytes and (ii) due representation of the biochemical and physiological functions of hepatocytes as functional network states.

The initial list of enzymatic reactions of the network was assembled from two existing global reconstructions of the human metabolic network (Duarte et al, 2007; Ma et al, 2007) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al, 2008). A chemical reaction or membrane transport was included into the hepatocyte-specific sub-network if we could find biochemical evidences for the presence of the respective process (i.e. the chemical conversion or the intercompartmental transport) or other evidences (e.g. genetics, proteomics) for the presence of the respective protein either in human hepatocytes or—if testing of the network's functionality required postulation of enzyme or transporter activity for which biochemical evidence in human hepatocytes was lacking—in the liver of other mammals. The information on which the decision was based is indicated in Supplementary Data 1. Using the METANNOGEN software (Gille et al, 2007) for the integration of information from various sources, >1500 primary literature articles, reviews and biochemical textbooks have been manually evaluated. Transport reactions were compiled on the basis of a comprehensive survey of the scientific literature on physiologically characterized transport processes in the liver. The draft network was used to initiate an iterative process of rigorous testing network functionality and focused on data re-evaluation and curation.

Informational content of the network

HepatoNet1 represents each reaction by chemically well-specified metabolites (given in Supplementary Data 1.1), each metabolite carrying a subscript that indicates its assignment to one of the six intracellular compartments (cytosol, endoplasmic reticulum and Golgi apparatus, lysosome, mitochondria, nucleus and peroxisome) and two extracellular compartments (bile canaliculus and sinusoidal space) (Table I). A total of 777 metabolites represent specific chemical entities and 69 describe pooled metabolites, that is composites of chemically similar molecular species (Supplementary information). A typical example is the pooled metabolite ‘triacylglycerols in VLDL', which is composed of different fatty-acyl residues in proportions as reported for a typical VLDL particle. Moreover, seven pseudo-metabolites were defined to take into account cofactor-using reactions that are not part of the network. For example, the pseudo-metabolite ‘NADH-redox potential' lumps together all cellular moieties that are not part of the network, but may either use or generate NADH, that is NADH<=>NADH-redox-potential+NAD+. Such a reaction permits to discriminate between the change in the redox state of NAD+ in which the sum of NAD+ and NADH remains constant from a net consumption or production of the pyridine/purine dinucleotide (e.g. consumption of NAD+ in ADP-ribosylations). For a definition of all seven pseudo-metabolites species, see Supplementary Table 1.

Table 1

Statistics of HepatoNet1

Item	Number
Reactions	2539
Membrane transporter	1466
Metabolites (total)	777
Cytosol	611
Mitochondrion	253
ER, Golgi apparatus	123
Lysosome	65
Nucleus	11
Peroxisome	79
Bile canaliculus	23
Sinusoidal space	255
Literature references	1527
Metabolic objectives
For validation	319
Physiological metabolic functions	123
For knock-out simulations

According to these types of metabolites and the level of mechanistic detail provided, the network comprises two main categories of reactions (Supplementary Data 1.2): (i) physical reactions, that is reactions and transport processes that interconvert chemically well-defined reactants (i.e. enzyme catalyzed reactions, free and facilitated diffusion) and (ii) abstract reactions, which as lumped processes describe the formation or degradation of large molecular complexes (e.g. lipoproteins). Physical reactions can be further classified into enzyme catalyzed reactions, spontaneous reactions, diffusion (sub-classes free and facilitated diffusion) and transport reactions (sub-class active transport). Whenever available, IUBMB EC numbers and TCDB-IDs are referenced for individual enzymatic and transport reactions, respectively. Abstract reactions comprise particle dis-/assembly, pool reactions, protein assembly, protein degradation, postulated transport reactions (i.e. uptake or secretion of compounds lacking stoichiometric detail), vesicular transport and virtual reactions used for definition of the aforementioned pseudo-compounds.

The network accounts for reactant and substrate stoichiometry and conservation of mass of each reaction and transport process. Charges are also balanced with the exception of protons because lacking knowledge on buffer capacities in the various cellular compartments supersedes the inclusion of bulk-phase protons in flux-balance analyses. It has to be noted, however, that protons involved in membrane transport processes are charge balanced as well.

HepatoNet1 is available in SBML (Supplementary Data 2) and tabular format (Supplementary Data 1). Definition, parameters and results of the validation process are contained in the Supplementary Data 3. For an overview of all supplements, see the ‘supplement guide'.

Functional validation and model performance

Our aim was to obtain a functionally competent network rather than a collection of all reactions reported for the hepatocyte. To this end, all metabolites known to be exchanged by the hepatocyte with bile and blood were functionally annotated (Supplementary Table 2). These exchange metabolites were used to prune the initially compiled hepatocyte-specific network to a reduced network that does not contain unreachable reactions, that is which cannot be forced to carry a non-zero flux in any of the metabolic objectives tested (Hoffmann et al, 2007).

We performed flux-balance analyses to validate whether for each of the 442 different metabolic objectives listed in Supplementary Data 3.2, a non-zero stationary flux distribution could be established in the network. In addition, we asserted that impossible tasks (e.g. albumin synthesis in the absence of phenylalanine and phosphorylation of ADP to ATP in the absence of any external energy source) are not accomplishable with the network.

Flux-balance computations were carried out by setting the target fluxes associated with the considered metabolic objective to a non-zero value and determining the stationary flux distribution that minimizes the sum of internal network fluxes (Holzhütter, 2004, 2006) and meets the criterion of thermodynamic realizability. The latter criterion assures concordance between the direction of fluxes and the Gibb's free energy of the reaction (Hoppe et al, 2007).

Functional testing of the network was performed at two different sets of metabolites that can be taken up or released by the network.

Minimal exchange set

In order to challenge the network to reveal its full functional capacity, we performed functional testing on a minimal set of importable nutrients (input set) containing besides the specific metabolite to be metabolized (e.g. lactate as substrate for the production of glucose) only glucose (in some simulations omitted, in which fatty acids or amino acids are forced to be the only nutrient), essential amino and fatty acids, vitamins, oxygen (in simulations at anoxic conditions omitted), sulfate, orthophosphate and water. In addition, we demanded that the metabolic objective tested has to be accomplished at strict homeostatic conditions, that is without release of other intermediates into the blood or bile other than the non-degradable end products of hepatic metabolism as, for example, CO2, urea or bilirubin diglucoronide. Compilation of the minimal input and output set was based on literature knowledge of essential nutrients and non-degradable metabolic end products of liver metabolism, respectively (see Supplementary Table 3.1).

Customized exchange sets

Here, the minimal input and export set is complemented by further metabolites according to the specific metabolic objective tested. For example, testing the objective ‘conversion of fatty acids into cholesterol', the input set was complemented by palmitate, whereas glucose was omitted (Supplementary Table 3). The type of exchange set used is indicated at each simulation given in Supplementary Data 3.

The 442 metabolic objectives outlined in Supplementary Data 3 can be accomplished by HepatoNet1. The respective flux modes (Supplementary Data 4) have been visualized and manually validated for their biochemical and physiological feasibility. For a given metabolic objective, the computed functional flux mode comprises a larger number of non-zero fluxes when the minimal exchange set is used rather than a physiological exchange set as all intermediates have to be broken down into no further degradable end products.

Analyzing metabolic objectives with HepatoNet1: bile acid synthesis in starvation

Bile acid synthesis being one important and hepatocyte-specific physiological liver function is described here as one example of the constraint-based modeling results obtained with HepatoNet1. Applying the flux minimization principle constrained by the minimal input and output set, we investigated the functional flux mode accounting for the synthesis of glycochenodeoxycholate (see Figure 1). The minimal input set contained palmitate instead of glucose to mimic systemic starvation conditions characterized by an increase of non-esterified fatty acids in the blood and continuous hepatic glucose output. The flux mode comprises 173 reactions including 46 transport reactions across the basolateral and canalicular plasma membranes and membranes separating the cytosol from mitochondria, endoplasmic reticulum and peroxisome. Palmitate, alanine, serine and O2 are the sole substrates. The pathways directly linked to the synthesis of the bile acid are the mevalonate pathway and the de novo synthesis of cholesterol. Other pathways such as gluconeogenesis, the pentose phosphate pathway or the ornithine cycle are part of the flux mode because the calculation was performed on the minimal input and output set, that is forcing all reactants to be balanced and all exportable intermediates to be catabolized into non-degradable end products. This example shows how the challenges of limited substrate supply and catabolism of intermediates to non-degradable excretion products leads to numerous cross-links between metabolic pathways traditionally perceived as separate entities. For example, alanine is used as gluconeogenetic substrate to form glucose-6-phosphate, which is used in the pentose phosphate pathway to generate NADPH. The glycine moiety for bile acid conjugation is derived from serine.

Figure 1

Functional flux mode for the synthesis of glycocholate. Flux distribution obtained for canalicular glycocholate as target flux using the minimal exchange set. Metabolites are labeled and color coded according to their sub-cellular localization (s, sinusoidal; c, cytosol; r, endoplasmic reticulum; m, mitochondrion; b, bile). All internal metabolites are balanced, including cofactors, which may appear in several instances in the graph. Reaction arrows correspond to flux direction and magnitude predicted by TR-FBA and reaction identifiers are indicated. The flux graph was created with the program CytoScape (Shannon et al, 2003; Killcoyne et al, 2009) and our own CytoScape plug-in FluxViz (König and Holzhütter, 2010).

Only a minor part of cytosolic β-HMG-CoA (the precursor of cholesterol) is produced from cytosolic acetyl-CoA (flux=0.29 mol/s per formation of 1 mol/s glycochenodeoxycholate) that is derived from citrate by the lipogenic enzyme ATP-citrate-lyase (Hoffmann et al, 1980). Indeed, only a small flux of citrate carbon into lipids had been shown by isotope experiments (Watson and Lowenstein, 1970). The major portion of cytosolic β-HMG is contributed by the ketone body acetoacetyl-CoA formed in the mitochondrion in equally large proportions from condensations of two moieties of acetyl-CoA (flux=2.86) and incomplete β-oxidation of palmitate (flux=2.35). This finding is in contrast to textbook lines stating acetyl-CoA to be the exclusive precursor of cholesterol. Mitochondrial acetoacetate is exchanged against cytosolic puruvate (Kummel, 1983, 1987) and activated by a cytosolic acetoacetyl-CoA synthetase (Bergstrom and Edmond, 1985). It has been shown that this predicted major route for the formation of cytosolic β-HMG-CoA from mitochondrial ketone bodies may indeed account for 19–80% of cholesterol synthesis (Endemann et al, 1982).

Analyzing metabolic objectives with HepatoNet1: detoxification of ammonia

Conversion of ammonia into non-toxic nitrogen compounds is one central homeostatic function of hepatocytes. Using the HepatoNet1 model, we investigated how the consumption of oxygen, glucose and palmitate is affected when an external nitrogen load is channeled in varying proportions to the non-toxic nitrogen compounds: urea, glutamine and alanine. We performed a high-resolution series of 20 301 flux minimization computations by varying the relative proportions of urea, glutamine and alanine in 200 steps between 0 and 100% and determined functional flux modes allowing oxygen, glucose and palmitate as substrates. A similar phase plane analysis has been used by Resendis-Antonio et al (2007) to analyze the physiological capabilities of the bacterium Rhibobium etli during different stages of nitrogen fixation. Oxygen demand is highest (molar ratio=0.407) if nitrogen is exclusively transformed into urea. At lower fluxes into urea, an intriguing pattern for oxygen demand is predicted: oxygen demand attains a minimum (0.105) if the nitrogen load is directed to urea, glutamine and alanine with relative fluxes of 0.17, 0.43 and 0.40, respectively (Figure 2A). Oxygen demand in this flux distribution is four times lower than for the maximum (100% urea) and still 77 and 33% lower than using alanine (0.186) and glutamine (0.140) as exclusive nitrogen compounds, respectively. This computationally predicted tendency is consistent with the notion that the zonation of ammonia detoxification, that is the preferential conversion of ammonia to urea in periportal hepatocytes and to glutamine in perivenous hepatocytes, is dictated by the availability of oxygen (Gebhardt, 1992; Jungermann and Kietzmann, 2000).

Figure 2

Simulation of oxygen and substrate demands for ammonia detoxification. Analysis of oxygen demand for ammonia detoxification is based on the simulation ‘NH3 degradation' (Supplementary Data 3.1) and the following settings: weight of the palmitate and glucose import is set to 10 and 100, respectively. The import flux of ammonia is set to 1. Export flux for urea iterates from 0 to 0.5 in 200 steps. Correspondingly, the export flux of glutamine ranges from 0 to 0.5 minus the export flux of urea. The export flux of alanine is set to the remaining nitrogen atoms: 1–2 × (export flux of urea+export flux of glutamine). Thus, the export of every atom imported with ammonium is fixed into predefined shares. (A) Oxygen demand as a function of the relative proportion of urea, glutamine and alanine as nitrogen compounds (arbitrary units). (B) Substrate demand, sum of glucose and palmitate fluxes. Fluxes are weighted by the number of carbon atoms in the substrate molecules, that is 6 (glucose) and 16 (palmitate). Plot axes as in (A).

The decreased oxygen demand in flux distributions using higher proportions of glutamine and alanine is accompanied by increased uptake of the substrates glucose and palmitate (Figure 2B). This is due to an increased demand of energy and carbon for the amidation and transamination of glutamate and pyruvate to discharge nitrogen in the form of glutamine and alanine, respectively. A detailed analysis of the flux modes obtained at high oxygen (89.3% urea, 5.3% glutamine, 5.3% alanine) and low oxygen (17.3% urea, 42.7% glutamine, 40.0% alanine) demand settings (Figures 3 and 4, respectively) revealed that both solutions rely on glycolysis, the oxidative pentose pathway, fatty-acid β-oxidation and the Krebs cycle for substrate entry and breakdown. In both simulations, non-zero fluxes through a complex part of the central metabolism occur consisting of (i) reactions and intermediates attributable to a single pathway; for instance, the urea cycle or Krebs cycle, (ii) reactions and intermediates interconnecting these traditionally perceived pathways (e.g. transaminase reactions, pyruvate) and (iii) transport processes for the flow and balancing of substrates, energy and redox equivalents between compartments. However, there are considerable differences between the two scenarios: the flux solution for high urea production, which is characterized by high oxygen demand, carries a high flux through the urea cycle and the recovery pathway of aspartate as a precursor for argininosuccinate synthesis. A sequence of reactions transports fumarate into mitochondria, produces and transaminates oxaloacetate (using Krebs cycle enzymes and aspartate transaminase) and transports aspartate to the cytosol. These reactions carry a considerably lower flux in the low oxygen demand solution, which is, conversely, characterized by a marked increase in glycolytic flux and exclusive use of pyruvate for transamination to alanine and formation of malate. In the high urea setting, pyruvate is primarily transported to the mitochondria and converted to acetyl-CoA and fed into the Krebs cycle. Intriguingly, in the low oxygen setting, not only is the flux through the complexes of the respiratory chain diminished, but the F0/F1-ATPase reaction does not carry any flux at all, indicating that in this solution the ATP demand can be satisfied entirely by substrate-level phosphorylation and that the respiratory chain and oxygen uptake is needed solely to drive secondary active transport processes. This finding is consistent with the observations that the proton motive force in hepatocytes is used only partly for ATP synthesis as approximately 20% of respiratory oxygen consumption can be attributed to proton-driven mitochondrial transport (Brand et al, 1991). With regard to the inversely correlated demand for oxygen and organic substrates, we observe several combinations of the chosen nitrogen compounds, which appear equally optimal with respect to effective substrate usage: high proportions of urea require high oxygen uptake, but lower quantities of glucose and palmitate, whereas flux distributions with proportions of nitrogen compounds as 17% urea, 43% glutamine and 40% alanine use considerably less oxygen at the expense of a higher substrate intake. For other nitrogen compound ratios, notably settings with very high flux into alanine, the reduction in oxygen consumption is accompanied by an even higher increase in substrate demand. Taken together, this analysis reveals strong dependencies between the available level of oxygen and variations in the substrate demand of hepatocytes required for effective ammonia detoxification by the liver.

Figure 3

Functional flux mode for the detoxification of NH3 with high oxygen demand. Flux distribution obtained for ammonia detoxification into the nitrogen compounds urea (89.3%), alanine (5.3%) and glutamine (5.3%). The setting is detailed in the result section (detoxification of ammonia). Presentation is similar to Figure 1.

Figure 4

Functional flux mode for the detoxification of NH3 with low oxygen demand. Flux distribution obtained for ammonia detoxification into the nitrogen compounds urea (17.3%), alanine (40%) and glutamine (42.7%). The setting is detailed in the result section (detoxification of ammonia). Presentation is similar to Figure 1.

Comparing the functional capacity of HepatoNet1 with that of other network reconstructions

In order to compare the functional capacity of HepatoNet1 with that of other available network reconstructions, we selected from the full list of tested metabolic objectives a comprehensive collection of 30 metabolic objectives that are either highly liver specific or known to be accomplishable in many different human cell types (Table II). Most of these objectives were evaluated using several simulations representing different physiological variants (e.g. gluconeogenesis from various precursors). All biochemical objectives listed in Table II are successfully accomplished by the HepatoNet1 model with the associated simulations yielding biochemically feasible functional modes (see Supplementary Table 4). Figure 5 illustrates the number of transporters and compartmentalized chemical reactions involved in the 123 functional flux modes. The smallest flux mode (formation of glucose from glycogen) comprises only 11 reactions and 2 compartments (cytosol and endoplasmic reticulum). The largest flux mode was obtained for the objective ‘VLDL synthesis' (247 reactions including 82 transporters). The flux modes also differ substantially in the number of involved compartments. Anaerobic re-phosphorylation of ATP takes place exclusively in the cytosol, whereas five different compartments are involved in the accomplishment of the objective ‘LDL degradation'.

Table 2

Representative set of metabolic objectives that can be accomplished by HepatoNet1

Metabolic function	Physiological significance (examples)	Simulations (see Supplement Table 4)
Re-phosphorylation of nucleoside triphosphates (ATP, GTP, CTP, TTP, UTP)	Cellular energy metabolism, principal cofactors in metabolism	(1)–(9)
Regeneration and de-novo synthesis of NADH and NADPH	Energy metabolism, reductive biosynthesis, redox homeostasis, principal cofactors in metabolism	(21)–(22)
Glycogenesis	Systemic glucose homeostasis; short-term regulation of blood glucose	(29)
Glycogenolysis	Systemic glucose homeostasis; short-term regulation of blood glucose	(30)
Gluconeogenesis	Systemic glucose homeostasis: homeostasis of blood glucose in post-prandial state; cori cycle, alanine cycle	(23)–(28)
Galactose and fructose metabolism	Breakdown of dietary carbohydrates	(31)–(32)
Formation of nucleotide-activated sugars	Glycoprotein and proteoglycan biosynthesis; biotransformation (conjugation)	(33)–(39)
Fatty-acid biosynthesis	Precursor of triglycerides and phospholipids; lipoprotein assembly and remodeling	(97)–(99)
Triglyceride biosynthesis	Systemic lipid homeostasis: cellular lipid store; lipoprotein assembly	(102)
Ketogenesis	Preferred substrate for the brain function after prolonged fasting	(105)–(106)
Farnesylpyrophosphate biosynthesis	Prenylation of proteins; ubiquinone biosynthesis; cholesterol precursor	(104)
Cholesterol biosynthesis	Membrane synthesis; lipoprotein assembly; bile formation; systemic lipid homeostasis	(103)
VLDL formation	Systemic lipid homeostasis: provision of endogenously derived triacylglycerides and cholesterol to peripheral tissue	(107)
LDL catabolism	Systemic lipid homeostasis: lipoprotein remodeling	(108)
Phospholipid biosynthesis	Membrane synthesis; lipoprotein and bile formation	(82)–(86)
Sphingolipid biosynthesis	Membrane synthesis	(87)–(89)
Salvage of purine and pyrimidine nucleotides	RNA and DNA biosynthesis; liver regeneration after injury and intoxication; energy metabolism	(18)–(20)
De novo synthesis of purine and pyrimidine nucleotides	RNA and DNA biosynthesis; liver regeneration after injury and intoxication; energy metabolism	(10)–(17)
Formation of non-essential amino acids	Systemic amino-acid homeostasis; precursors of endogenous proteins, plasmaproteins and proteoglycan; creatine and glutathione synthesis	(42)–(53)
Ureogenesis	Ammonia detoxification; systemic pH regulation, cori cycle, alanine cycle, systemic nitrogen and amino-acid homeostasis	(78)–(79)
Creatine biosynthesis	Central metabolite for energy metabolism of muscle cells	(80)
Heme biosynthesis	Biosynthesis of heme proteins, in particular cytochrome P450 (biotransformation and metabolic detoxification)	(81)
Protein biosynthesis	Enzymes, cytoskeleton, membrane proteins, plasma proteins (e.g. acute-phase response, pro-hormone synthesis)	(77)
Amino-acid breakdown	Systemic amino-acid homeostasis; carbon sources for gluconeogenesis, ketogenesis, alanine cycle	(54)–(76)
Biosynthesis of cofactors	Principal components of intermediate metabolism, biotransformation and detoxification	(117)–(118)
Bile formation	Cholesterol homeostasis; intestinal uptake and digestion of dietary lipids, transport of lipophilic vitamins, intestinal motility	(109)–(116)
Formation of glutathione	Intracellular and systemic redox homeostasis, systemic GSH supply; anti-oxidant function; cofactor in biotransformation	(119)
Bilirubin catabolism	Heme catabolism	(121)
Detoxification of reactive oxygen species	Defense against oxidative stress	(120)
Detoxification of xenobiotics	Metabolic detoxification	(122)–(123)

Figure 5

Composition of functional flux modes related to the metabolic objectives in Table II. The heights of the colored bars indicate the number of membrane transporters and compartmentalized chemical reactions constituting the complete flux mode. For the metabolic objectives underlying the 123 flux modes, see Supplementary Table 4.

We then tested how the biochemical objectives shown in Table II are represented in the global human metabolic network reconstruction, Recon1 (Duarte et al, 2007). When comparing the HepatoNet1 model to Recon1, 799 localized metabolites and 731 reactions can be found in both metabolic networks (Supplementary Data 1), whereas HepatoNet1 contains 1709 enzymatic and transport processes for which no equivalent to Recon1 could be found. This is partly a consequence of differences in the representation of enzymatic processes. For instance, the fatty-acyl residue distribution of lipids in HepatoNet1 is tailored for the analysis of liver-specific questions. Thus, uptake, formation and excretion of lipoprotein particles or bile-specific phospholipids and the corresponding reactions do not coincide between the networks, although the same processes are covered. Accounting for such differences in compound representation, 119 of the 123 simulations in Table II could be directly compared. As expected, for most objectives Recon1 allowed feasible flux modes. The remaining 14 cases in which Recon1 failed are mostly because of the lack of liver-specific reactions. For example, a considerable number of hepatic transport processes, either intracellular or at the basolateral and canalicular plasma membranes, are unique to HepatoNet1.

Next, we tested a liver-specific prediction of metabolic activity in Recon1 based on transcription profiles (Shlomi et al, 2008). Taking into account all reactions of Recon1 and omitting those reactions that according to this prediction should not be active in the liver (indicated by a negative confidence level), only in 71 of 119 simulations, a stationary flux distribution could be obtained (see Supplementary Table 4). This result shows the current limitations inherent in the prediction of species- or tissue-specific networks derived from transcript data alone. The global reconstructions of human metabolism assisted by computational methods are a useful starting point to establish tissue-specific networks. However, currently, there is no alternative to the additional manual reconstruction and validation in order to attain a physiologically functional model.

Using HepatoNet1 to explore the robustness of metabolic liver functions against enzyme deficiencies

An intriguing question that can be addressed with large-scale networks pertains the essentiality of enzymes and transporters (Behre et al, 2008; Suthers et al, 2009; Plaimas et al, 2010). Using HepatoNet1, we have analyzed how the knock-out of a single enzyme or transporter may compromise the metabolic objectives accomplished by the human hepatocyte. To this end, we have performed 123 computational knock-out studies in which we disabled one after the other each single reaction of the original 123 functional flux modes related to the metabolic objectives (Table II) and counted how often a computational knock-out could be compensated by an alternative flux mode (=non-essential reaction). As the essentiality of an enzyme or transporter depends on the availability of external substrates that can be used to bypass the invalidated reaction, these knock-out studies were performed on less restrictive customized exchange sets specified in Supplementary Data 3. Hence, the number of essential enzymes detected in these knock-out simulations will be lower than at constraints imposed by the minimal exchange set, but probably higher than at comfortable conditions in which the hepatocyte has full access to all plasma metabolites.

Ranking the essential reactions determined in these 123 knock-out studies in descending order according to the frequency of their occurrence allows to discriminate between reactions with cardinal essentiality that are indispensable for almost all metabolic objectives tested (e.g. electron carriers of the respiratory chain) and reactions that are essential for only a limited (usually closely related) number of metabolic objectives (Supplementary Table 5). For 80 enzymes and transporters that turned out to be essential in at least one knock-out simulation, enzymopathies with clinical symptoms have been reported. To our surprise, clinically manifested enzymopathies have been reported even for those enzymes and transporters, which in case of a complete knock-out are predicted to impair a larger set of metabolic objectives. This may point to the existence of protein isoforms that are not affected by the specific type of deficiency. Alternatively, one may speculate that non-lethal but clinically manifested enzymopathies are restricted to moderate defects that still allow significant fluxes.

Discussion

Human physiology relies on the integration of specialized metabolic functions in different tissues and cell types. Here, a system-scale stoichiometric model of human hepatocyte metabolism is presented that accounts for a large variety of biochemical functions of this cell type. The model is based on a comprehensive evaluation of the currently available knowledge on hepatocyte metabolism combined from several information sources and lines of evidence. This metabolic network is taken as the starting point of a computational protocol and validation strategy that establishes a network with a functional scope determined by a comprehensive list of hepatic physiological functions and documented transport activities at the cell boundaries. The model is confined to reactions, which can actually carry flux, exchanging only those metabolites that are contained in the set of exchangeable metabolites across the system boundary.

Metabolic objectives representing synthetic capabilities of the model can be fulfilled using the minimal input set that contains only compounds known to be essential nutrients for human beings. The implications of the restrictive nature of this model boundary exchange imposed in our model calculations are two-fold. First, consistent with biochemical knowledge, it shows that the model faithfully represents the anabolic and catabolic pathways and metabolic versatility of hepatocytes needed to synthesize and degrade various components of intermediary metabolism starting from or leading to simple input and output metabolites. Second, the feasibility of steady state flux distributions under very restrictive exchange conditions with granting or prohibiting exchange of individual compounds suggests that the capacity of hepatocyte metabolism for homeostatic regulation of blood composition may allow for the adjustment of perturbations for individual compounds without major interference among blood constituents.

Computational testing of the network's functional capabilities was performed by using a variant of flux-balance analysis that assures thermodynamically feasible flux directions (Hoppe et al, 2007) and at the same time obeys the principle of flux minimization (Holzhütter, 2004, 2006), that is minimizes the sum of internal fluxes. Hitherto, constraint-based modeling approaches for study of metabolism is much more advanced in unicellular systems (Feist and Palsson, 2008) than in cell types with more complex architecture. In microorganisms, simple objective functions (e.g. biomass production) can often successfully predict cellular behavior (Edwards et al, 2001). However, for the metabolism of differentiated somatic cells as hepatocytes that meet many different objectives in parallel, one may think of a variety of criteria (Nagrath et al, 2007; Uygun et al, 2007), which may govern the regulation of the network and thus determine the magnitude of network fluxes. The minimization of internal fluxes is only one of these criteria. It relies on the plausible assumption that the cell aims to accomplish its metabolic tasks at a minimum expense of enzyme capacity and external resources. Our findings with respect to the oxygen demand for ammonia detoxification indicate that the principle of flux minimization may indeed yield predictions that are in line with experimental observations on the sinusoidal zonation of metabolic pathways (Gebhardt, 1992; Haussinger et al, 1992). Moreover, functional flux modes by definition represent a minimal set of non-zero fluxes required to meet a single metabolic objective and thus provide a glance at the minimal complexity of the flux-carrying sub-network involved. Nevertheless, critical evaluation of the feasibility of the flux minimization criterion used has yet to be performed by confronting calculated flux distributions with experimentally determined flux rates.

In terms of both scope and specificity, our model bridges the scale between models constructed specifically to examine distinct metabolic processes of the liver and modeling based on a global representation of human metabolism. The former includes models for the interdependence of gluconeogenesis and fatty-acid catabolism (Chalhoub et al, 2007), impairment of glucose production in von Gierke's and Hers' diseases (Beard and Qian, 2005) and other processes (Calik and Akbay, 2000; Stucki and Urbanczik, 2005; Ohno et al, 2008). The hallmark of these models is that each of them focuses on a small number of reactions pertinent to the metabolic function of interest embedded in a customized representation of the principal pathways of central metabolism. In contrast to either of the aforementioned modeling scales, HepatoNet1 provides the combination of a system-scale representation of metabolic activities and representation of the cell type-specific physical boundaries and their specific transport capacities. This allows for a highly versatile use of the model for the analysis of various tissue-specific physiological functions. Conceptually, from a biological system perspective, this type of model offers a large degree of comprehensiveness while retaining tissue specificity, a fundamental design principle of mammalian metabolism.

Recently, a computational method for network-based prediction of human tissue-specific metabolism using tissue-specific transcript and protein abundance data has been introduced (Shlomi et al, 2008). The method relies on the maximization of the number of reactions in the metabolic network whose predicted flux activity is consistent with the expression of the genes of the corresponding enzymes. We were able to establish only 71 of the 119 metabolic functions of human hepatocytes accomplishable with Recon1 when we excluded those reactions for which inactivity was predicted by Shlomi et al (2008), that is reactions carrying a score of less than zero. Several reasons may explain this negative outcome. One reason is certainly the poor coverage of the large-scale network with expression data. We have used gene expression data as one of the several lines of evidence in the course of curation of HepatoNet1. Only 55% of all reactions could be associated with gene expression data. Thus, for a substantial number of reactions in our network, gene expression data is not assigned or available. Generally, expression profiles were consistent with biochemical evidence, although we identified a small but important number of exceptions in which clear biochemical evidence for enzyme activity was not reflected in liver-specific transcript abundance (Supplementary Data 5; Supplementary Table 6, respectively). Another reason restraining the straightforward definition of tissue-specific networks on the sole basis of expression profiles consists of the poor overall correlation of transcript abundance and metabolic flux strength (Ovacik and Androulakis, 2008; Ponten et al, 2009). For example, in the set of expression data used by Shlomi et al (2008), many glycolytic enzymes are marked with ‘not expressed', although the glycolytic flux in human hepatocytes is large (Cherrington, 1999). Glycolytic enzymes in hepatocytes have to be continuously present, either for the purpose to convert extra glucose into triglycerides or to produce glucose from gluconeogenic precursors in fasting conditions.

HepatoNet1 incorporates different sources of information, which provides a greater degree of confidence. This avoids problems, which may occur if just one type of biological evidence is used, especially as the use of transcript abundance as a measure for enzyme activity relies on a number of simplifying assumptions, which, in sum, may significantly alter the structure and behavior of the model. The use of different types of biological evidence and manual integration and validation of the model components in a functional network context have been vital for obtaining a model that provides adequate resolution of hepatocyte biochemical functions and alignment with physiologically determined input and output parameters.

HepatoNet1 is expected to provide a structural platform for computational studies on liver function. The results presented herein highlight how internal fluxes of hepatocyte metabolism and the interplay with systemic physiological parameters can be analyzed with constraint-based modeling techniques. At the same time, the framework may serve as a scaffold for complementation of kinetic and regulatory properties of enzymes and transporters for analysis of sub-networks with topological or kinetic modeling methods.

Materials and methods

Network reconstruction

Biochemical reactions and their assignment to metabolic pathways in the global reconstructions of the human metabolic network (Duarte et al, 2007; Ma et al, 2007) and the KEGG (Kanehisa et al, 2008) were used as a starting point for the initial candidate list of network components. For each reaction, various enzyme and protein databases such as the Braunschweig Enzyme database (BRENDA) (Chang et al, 2009), Reactome (Matthews et al, 2009) and UniProtKB (The UniProt Consortium, 2008) have been used. Relevant content was analyzed to evaluate the available physiological, biochemical and genetic evidence with special regard to tissue specificity and sub-cellular localization. A reaction was included in the network if (i) experimental evidence was available for the reaction to occur in human hepatocytes or liver tissue in general or (ii) similar experimental evidence from other mammalian species and consideration of human orthologous genes allowed inference of the reaction. Biochemical (enzyme assay, transport assay, protein expression and localization) and genetic evidence (genes, transcript expression and genetic diseases) were the principal types of evidence used. Although generally, enzymatic and transport reactions represent mechanistic biochemical processes, the network contains a number of lumped reactions, which represent processes that are mechanistically beyond the scope of HepatoNet1 and describe its metabolic inputs and outputs. Glycogen is represented on the basis of the corresponding species in Recon1 (Duarte et al, 2007). An exhaustive search of the Gene Expression Omnibus (NCBI-GEO; Barrett et al, 2009) revealed a total of 12 datasets containing gene expression data from normal human liver tissue samples. Eight of these, GDS181 (Su et al, 2002), GDS422 to 426 (Yanai et al, 2005), GDS1096 (Ge et al, 2005) and GDS1209 (Yoon et al, 2006) include affymetrix detection calls and were used as additional lines of evidence for network reconstruction. Only samples of normal human liver tissue/hepatocytes were considered. Affymetrix probe set IDs were mapped to reactions using the Ensembl Homo sapiens database (Hubbard et al, 2009) and KEGG orthology records.

Reaction directionality

Principally, reactions are treated as reversible and constraints on reaction directionality are only imposed systematically in the course of computational analysis of the network. Empirical constraints of directionality have been manually set for a limited set of reactions to prevent unfeasible flux distributions. See Supplementary Data 1.2 for details.

Compartmentalization

Endoplasmic reticulum, microsomes and Golgi apparatus were modeled as one compartment because in many cases, experimental data do not allow to discriminate between these compartments. Mitochondrial intermembrane space was assigned to the cytosolic compartment. Sub-cellular localization of reactions was determined based on direct experimental evidence (protein localization, targeting sequences and subcellular fractionation) and indirect physiological or biochemical evidence. In the absence of information, reactions were modeled as cytosolic.

Structural network analysis

Atomic mass balance was analyzed by computation of net stoichiometry (Gevorgyan et al, 2008). Protons involved in transport processes across membranes occur as separate species in the model and are mass balanced, which is not the case for protons exchanged with the aqueous bulk phase.

Pruning

Reconstructed network data was subjected to functional pruning (Hoffmann et al, 2007) by confining the model to reactions that can potentially carry a non-zero flux. The hepatocyte-specific set of exchange processes used for pruning comprises substrates and end products of hepatocyte biosynthetic and catabolic activities and was assembled based on liver physiological and biochemical functions (Supplementary Table 2).

FBA

Pruned network models were subjected to flux-balance simulations setting the target fluxes according to the demanded metabolic objective (Supplementary Data 3) to non-zero positive or negative values and determining the stationary flux distribution. The optimization objective has been the minimization of internal fluxes (Holzhütter, 2004). Fluxes are weighted equally by default unless modified in selected cases to reflect differential activity of enzymes with alternative substrates or cofactors (e.g. affinity of hexokinase toward various hexoses) (Supplementary Data 6.1). Thermodynamic constraints on reaction directionality have been imposed (Hoppe et al, 2007). For this purpose, standard Gibb's energies listed in Supplementary Data 1.2 have been obtained from a prediction method (Jankowski et al, 2008) and physiological metabolite concentration ranges listed in Supplementary Data 6.2 from the Human Metabolome database (Wishart et al, 2007). Computation was performed with the aid of CPLEX 10.1 (ILOG, Gentilly, France).

In an iterative process, simulation results that were not in agreement with the expected outcome were used to spot inconsistencies in the underlying network structure. These inconsistencies were the subject of curation, followed by new simulations to revalidate HepatoNet1. Flux modes were visualized with BiNA 1.3.1 (Kuentzer et al, 2007) or CytoScape (Shannon et al, 2003) in combination with FluxViz (König and Holzhütter, 2010).