INTRODUCTION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

The Human Genome Project (HGP) is now essentially complete. One result of the HGP is the definition of single nucleotide polymorphisms (SNPs) and their effects on the development of human disease. Although the number of SNPs in the human genome is expected to be a few million, it is estimated that only 100,000 to 200,000 will effectively define a unique human genotype. A subset of these SNPs is believed to be "informative" with respect to human disease (Syvanen 2001). Many of these "informative" SNPs will fall into coding regions while others will be found in regulatory regions. The human genotype-phenotype relationship is very complex and it will be hard to determine the causal relationship between sequence variation and physiological function. One way to deal with this intricate relationship is to build large-scale in silico models of complex biological processes (Fig. 1). Defects or alterations in the properties of a single component in complex biological processes can be put into context of the rest by using an in silico model. In this work, recent data on SNPs in key red blood cell enzymes (Fig. 1a) and corresponding alterations in their kinetic properties (Fig. 1b) were used in an in silico red blood cell model (Fig. 1c) to calculate the overall effect of SNPs on whole cell function (Fig. 1d).

View larger version (32K): [in this window] [in a new window]   Figure 1   A schematic depiction of the genotype-phenotype relationship; from sequence variation (a) to changes in biochemical function (b) to network properties (c) to overall physiological function (d). The upper path signifies a normal pathological situation. The DNA sequence codes a nondeficient enzyme that catalyzes the glucose-6-phosphate dehydrogenase or pyruvate kinase reactions normally resulting in homeostatic red blood cell metabolism with the cell assuming its traditional biconcave shape. The lower path represents a pathological situation, a defect in the DNA sequence that leads to a metabolic enzyme with altered catalytic activity. The result is a loss in the cell's ability to respond to oxidative and/or energy loads, hence an inability to maintain osmotic balance and stability that ultimately results in lysis-the phenotypic expression of the genetic defect.

Erythrocyte metabolism has been extensively studied. Over the years, a series of mathematical models of increasing scope and detail of erythrocyte metabolism have been formulated and they have been used to describe the red blood cell's integrated metabolic functions (Rapoport and Heinrich 1975; Heinrich et al. 1977; Schauer et al. 1981; Joshi and Palsson 1989; Mulquiney and Kuchel 1999). Such models can be used to describe the consequences of altered kinetic constants on red blood cell function (Heinrich et al. 1977; Holzhütter et al. 1985; Schuster 1995; Jacobasch and Rapoport 1996; Martinov et al. 2000). The study of variations in the kinetic properties of red blood cell enzymes is not merely an academic study of the quality of a mathematical model, but has real utility in the clinical diagnosis and treatment of enzymopathies and can provide a link to the underlying sequence variation (Fig. 1). Here, an in silico model is used to study SNPs in two of the most frequent red blood cell enzymopathies: glucose-6-phosphate dehydrogenase (G6PD) and pyruvate kinase (PK).

	RESULTS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

For both enzyme deficiencies, clinical data was obtained from the published literature to determine measured values for the various kinetic parameters (V_maxs, K_ms, K_is) associated with each clinically diagnosed variant. These numerical values were then used in the in silico model (Jamshidi et al. 2001) and sensitivities to various oxidative and energy loads (above normal, baseline values) were simulated. The results are interpreted with respect to the genetic basis of the enzymopathy in an attempt to establish a direct link between the genotype and phenotype (Fig. 1).

G6PD Deficiency

G6PD catalyzes the first step in the oxidative branch of the pentose pathway (Fig. 1c) and is thus of critical importance in maintaining the red blood cell's resistance to oxidative stresses. G6PD has been studied extensively since its initial discovery and characterization in the early 1960's (Kirkman 1960; Kirkman 1962; Thorburn and Kuchel 1985; Kirkman and Gaetani 1986). It is the most common erythrocyte enzymopathy, affecting ~400 million people worldwide. The World Health Organization (WHO) has standardized methods by which to describe these deficiencies, and over 400 different biochemical variants have been described based on the biochemical properties (e.g., electrophoretic mobility, V_max, etc.) of the enzymes (Beutler 1990). Four classes have been defined ranging from class I (chronic nonspherocytic hemolytic anemia) to class IV (no deficiency) grading the level of hemolytic anemia observed in patients having G6PD deficiencies (Beutler and Yoshida 1988; Beutler 1990; Jacobasch and Rapoport 1996).

G6PD deficiency is a sex-linked trait. The G6PD gene is located in position Xq28 in the long arm of the X chromosome, with 13 exons and 12 introns spanning 20 kb (Yoshida 1986; Beutler 1990; Jacobasch and Rapoport 1996). As a monomer, G6PD is a 59-kD protein, 515 amino acids long. It is active in the tetrameric (in acid) or dimeric (in base) forms (Adediran 1996). The monomer has two domains, an N-terminal domain (amino acids 27-200) and a + domain. The N-terminal domain contains a -- dinucleotide (NADP) binding site (amino acids 38-44). The -helix links the two domains and contains a highly conserved substrate-binding site (amino acids 198-206). The quaternary structure of the enzyme is important for its stability and its ability to maintain normal activity. Normal G6PD requires low concentrations of NADP to maintain its integrity (Beutler 1986). Thus, NADP serves two functions; it stabilizes the dimer and serves as a substrate. The same domain is believed to be responsible for structural stability and catalytic activity (Hirono et al. 1989).

G6PD from normal patients and patients with hemolytic anemia have been characterized on the molecular level. A total of 61 G6PD class I variants have been described at the molecular level. Of the 61 class I chronic variants, 55 are the result of SNPs involving amino-acid changes, five result from frame deletions and one results from a splicing defect (Fiorelli et al. 2000).

Clinically diagnosed SNPs cluster around important, active regions of G6PD enzyme including the dimer interface and substrate-binding sites (Fig. 2a). Numerical values of G6PD kinetic parameters were varied in silico to determine the sensitivity of red blood cell metabolic functions to these changes in enzyme function. The most sensitive parameters were found to be V_max and K_i-NADPH. The NADPH/NADP ratio proved to be the most informative indicator of metabolic status as it was the most sensitive to changes in these two parameters and it gives an indication as to the oxidative state of the cell (Kirkman et al. 1975).

View larger version (92K): [in this window] [in a new window]   Figure 2   The location of SNPs in the glucose-6-phosphate dehydrogenase (G6PD) protein. (a) The G6PD variants for which there is full molecular characterization available, are shown with their corresponding parameter values and amino-acid location of the single nucleotide polymorphism (SNP). Note that the SNPs are located in critical regions of the protein. (b) Shows that there is no obvious correlation between the altered numerical values for the two kinetic parameters in either the chronic or nonchronic cases.

For each documented variant, there appears to be no direct correlation between V_max and K_i-NADPH (Fig. 2b). Clinically, G6PD deficiencies are broken down into two main categories: chronic and nonchronic hemolytic anemia. Chronic cases show clinical symptoms and are very sensitive to the environment. Nonchronic cases appear normal under homeostatic conditions but can experience problems when subjected to large oxidative stresses (Jacobasch and Rapoport 1996). For this study, kinetic data for 12 chronic and 8 nonchronic cases from Yoshida and 19 chronic cases from Fiorelli were used (Yoshida 1986; Fiorelli et al. 2000). The G6PD variant data used from the references were characterized according to the WHO standards defined in Betke et al. (1967).

Under normal conditions (i.e., oxidative load, v_ox = 0), there are differences between the chronic and nonchronic groups, with the chronic group having a somewhat lower homeostatic steady-state NADPH/NADP ratio than the nonchronic group (data not shown). When subjected to an oxidative load (v_ox > 0), noticeable differences between the two groups (chronic and nonchronic) appear (Fig. 3). The NADPH/NADP ratio at the maximum tolerated oxidative load (v_ox = max value) correlates with this ratio in the unstressed situation (v_ox = 0). The group of chronic hemolytic anemia patients is clearly separated from the normal and nonchronic group. A number of the chronic cases can only withstand a very modest oxidative load. Of the variant cases studied, a handful has been characterized at the molecular (amino acid) level (Table 1). Of the cases considered, most of the single-base changes in the chronic (class I) variants occur at or near the dimer interface (exons 10,11 and 6,7) or near the NADP binding site, leading to an impaired ability to respond to systemic oxidative challenges.

View larger version (84K): [in this window] [in a new window]   Figure 3   The ability of glucose-6-phosphate dehydrogenase (G6PD) variants to respond to oxidative load based on in silico red blood cell model results. The NADPH/NADP is plotted for maximum load (vox = max value tolerated) versus no oxidative load (vox = 0) for each G6PD variant. The nonchronic cases cluster near the normal case and can sustain a maximum oxidative load of near normal (3 mM/h). The chronic cases can be split into two categories, those that can sustain a near-normal load (chronic high load) and those who are well below normal (chronic low load). The horizontal arrows demarcate two regions: an oxidized state (Log [NADPH/NADP] < 0) where the majority of the cofactor is in the form of NADP and a reduced state (Log [NADPH/NADP] > 0) where the majority of the cofactor is in the form of NADPH. For G6PD enzymopathies, the oxidized state is prevalent, while the reduced state is found in normal individuals (Kirkman et al. 1975). The insert shows no clear correlation between Vmax, Ki-NADPH, and the cell's ability to withstand oxidative loads.

PK Deficiency

Pyruvate kinase (PK) is a key glycolytic regulatory enzyme. There have only been ~400 documented variants since PKs first description in 1961 (Tanaka and Zerez 1990; Jacobasch and Rapoport 1996; Zanella and Bianchi 2000). PK accounts for 90% of the enzyme deficiencies found in red blood cell glycolysis (Jacobasch and Rapoport 1996). It is autosomal recessive where clinical manifestations appear only in compound heterozygotes (two mutant alleles). There are four isozymes: L, R, M₁, and M₂, with the R type being exclusive to red blood cells. PK is encoded by the PK-LR gene on chromosome 1q21. The kinetics of the enzyme have been extensively studied (Otto et al. 1974). PK activity is regulated by F6P, ATP, Mg, and MgATP. Anemic heterozygotes have 5-40% of normal PK activity.

Active PK exists as a tetramer of four 60 kD subunits. There are four domains: A [a ()₈ barrel], B (a small  barrel), the C-terminal, and the N-terminal with a helix-turn-helix motif. The amino-acid sequence of PK is highly conserved. Based on PK studies of the rabbit, Escherichia coli, and yeast, the PEP and cation binding sites are in the cleft between the A and B domains, close to the sixth and eighth loops of the ()₈ barrel. PK exhibits 222 symmetry, and there are two states for the tetramer, the low-affinity, tight (T) conformation and the high-affinity, relaxed (R) conformation. The PEP and Mg²⁺ binding sites are in the cleft between the A and B domains close to the sixth and eighth loops of the ()₈ barrel. The FDP binding site is in the C domain. The enzyme is stabilized by the interaction between the A and C domains of opposing subunits (Zanella and Bianchi 2000). While the crystal structure of human R-type PK has not yet been explicitly determined, it shares sequence homology with feline PK, whose crystal structure has been determined. Given this information along with computer simulation studies, it has been determined that some of the PK SNPs leading to severe hemolytic anemia are not in the catalytic sites. Some of these mutations are believed to occur at subunit interfaces (Solinge et al. 1997).

Unlike for G6PD, the characterized PK SNPs are scattered throughout the protein coding region and do not appear to cluster near the corresponding active site of the enzyme (Zanella and Bianchi 2000). A summary of the PK variants is presented in Table 2. The biochemical data used for the PK variants from these references were characterized according to the protocols described by Miwa et al. (1979). The documented kinetic values for the main kinetic parameters V_max and K_PEP are shown (Fig. 4a). Similar to the G6PD variants, there is not a clear correlation between changes in the numerical V_max and K_PEP amongst the PK variants (Fig. 4b). Although changes in K_ADP are also documented for each variant and accounted for in the simulations, increases or decreases in its value did not significantly affect the red blood cell's steady-state metabolite concentrations or its ability to withstand energy loads (data not shown). Changes in K_PEP and V_max influence the concentration of ATP and 2,3DPG most significantly. When increased energy loads (v_e >0) are applied in silico, differences between the variants are observed. The ratio between the ATP concentration at maximum tolerated load (v_e  = max value) and the ATP concentration in the unchallenged state (v_e = 0) varies approximately linearly with the maximum tolerated load when all the variants are evaluated (Fig. 4c). Thus, the variants that tolerated the lowest maximum load have an [ATP]_max/[ATP]_{no load} ratio close to unity indicating their sharply diminished ability to deviate from the nominal homeostatic state. Interestingly, the computed energy charge [EC = (ATP + 1/2ADP)/(ATP+ADP+AMP)] (Atkinson 1977) stays relatively constant (Fig. 4d). This result indicates that red blood cell metabolism strives to maintain its EC within the tolerated load range, thus allowing for an energetically consistent metabolic function.

View larger version (67K): [in this window] [in a new window]   Figure 4   The ability of pyruvate kinase (PK) variants to respond to energy load. (a) The numerical values for Vmax and KPEP for the documented PK cases. (b) There is no obvious correlation between the two kinetic parameters amongst the variants studied. (c) As the maximum tolerated ATP load decreases, the ratio of the cells [ATP] at maximum load versus at no load increases resulting in a near linear relationship. (d) The energy charge remains near constant despite the significant drop in the each variant's ability to withstand an ATP load.

	DISCUSSION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

With the HGP effectively complete and the amount of sequence variation information growing, the need to define the consequences of such variations on physiological function is brought to the forefront. Here, we use the human red blood cell to show that an in silico model can be used to aid in this process.

Sequence variations in coding regions for metabolic enzymes can lead to altered kinetic properties. The kinetic properties of enzymes are described by many parameters, and a single SNP can alter one or many of these parameters. For the variants of G6PD and PK considered here, there appears to be no clear relationship between their kinetic parameters as a function of sequence variation. Thus, consequences of sequence variations on the function of a gene product must be fully evaluated to get a comprehensive assessment of the altered biochemical function.

The consequences of many simultaneously altered enzyme properties must in turn be evaluated in terms of the function of the enzyme in the context of the reaction network in which it participates. The assessment of sequence variation on biochemical and kinetic properties of enzymes may seem difficult and this challenge is currently being addressed (Yamada et al. 2001), but the assessment of sequence variation on entire network function is even more complicated. This highly complex and intricate relationship between sequence variation and network function can be studied through the use of a computer model. Here, we have shown that a large number of variants in red blood cell G6PD and PK can be systematically analyzed using an in silico model of the red blood cell. Correlations between sequence variation and predicted overall cell behavior is established, and in the case of G6PD, it in turn correlates with the severity of the clinical conditions.

The red blood cell in silico model has been developed through successive iterations over the past 30 years (Rapoport and Heinrich 1975; Heinrich et al. 1977; Schauer et al. 1981; Holzhütter et al. 1985; Joshi and Palsson 1989; Mulquiney and Kuchel 1999). Building in silico models of complex biological processes has therefore proven to be a difficult task. This challenge may now be more manageable given the development of data-driven constraints-based models (Schilling et al. 1999; Palsson 2000). We thus may be entering an era where in silico biology will develop computer representations of complex biological processes, and these in silico models may then, in turn, aid in the analysis, interpretation, and even prediction of intricate genotype-phenotype relationships, both for normal physiological and pathophysiological conditions.

	METHODS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

The metabolic network of the mature erythrocyte is composed of glycolysis, the pentose phosphate shunt, the 2,3-diphosphoglycerate (2,3DPG) shunt, and the nucleotide salvage pathway. Glycolysis produces ATP by substrate-level phosphorylation to drive the Na/K pump that enables the cell to maintain osmotic balance and its biconcave shape. The pentose phosphate pathway produces NADPH to keep glutathione in the reduced state and allows the cell to counteract oxidative loads (via oxygen radicals, superoxides, etc.). The two enzyme deficiencies considered here involve key fluxes in these two key metabolic pathways.

Simulation studies with a model of red blood cell metabolism were performed on a PC using Mathematica (Wolfram Research; Jamshidi et al. 2001). The model includes glycolysis, the pentose-phosphate pathway, the Rapoport-Leubering (2-3DPG) shunt, nucleotide synthesis and degradation reactions, the ATPase pump, equilibrium expressions accounting for magnesium complexing, and osmotic and electroneutrality balances. This in silico model is composed of 41 differential and algebraic equations that must be solved simultaneously (Joshi and Palsson 1989; Joshi and Palsson 1990; Lee and Palsson 1991).

Although the model can compute dynamic metabolic responses, here we focus solely on the steady-state solutions. The numerical value for any kinetic parameter in the model can be varied and the resulting steady state calculated. If the numerical value of a kinetic parameter is varied too far from its normal value, a point will be reached where there is no solution to the system of equations and the in silico cell is not able to reach a steady state. This point corresponds to uncontrolled cell swelling and lysis.

Loads on energy (v_e) and redox (v_ox) metabolism are simulated in silico by assigning increased fluxes to the following nonspecific conversions of cofactors to simulate elevated metabolic loads (Fig. 1c),

which are zero in the normal physiological state. When v_e or v_ox have reached values that lead to cell lysis, they are scored as the maximum flux loads sustainable by the cell. This maximum will vary with the numerical values of the kinetic parameters that in turn may be influenced by SNPs.