Electrostatic potential on human leukocyte antigen: implications for putative mechanism of chronic beryllium disease.

The pathobiology of chronic beryllium disease (CBD) involves the major histocompatibility complex class II human leukocyte antigen (HLA). Although occupational exposure to beryllium is the cause of CBD, molecular epidemiologic studies suggest that specific (Italic)HLA-DPB1(/Italic) alleles may be genetic susceptibility factors. We have studied three-dimensional structural models of HLA-DP proteins encoded by these genes. The extracellular domains of HLA-DPA1*0103/B1*1701, *1901, *0201, and *0401, and HLA-DPA1*0201/B1*1701, *1901, *0201, and *0401 were modeled from the X-ray coordinates of an HLA-DR template. Using these models, the electrostatic potential at the molecular surface of each HLA-DP was calculated and compared. These comparisons identify specific characteristics in the vicinity of the antigen-binding pocket that distinguish the different HLA-DP allotypes. Differences in electrostatics originate from the shape, specific disposition, and variation in the negatively charged groups around the pocket. The more negative the pocket potential, the greater the odds of developing CBD estimated from reported epidemiologic studies. Adverse impact is caused by charged substitutions in positions 55, 56, 69, 84, and 85, namely, the exact same loci identified as genetic markers of CBD susceptibility as well as cobalt-lung hard metal disease. These findings suggest that certain substitutions may promote an involuntary cation-binding site within a putatively metal-free peptide-binding pocket and therefore change the innate specificity of antigen recognition.

The recognition that CBD appeared to affect a subset of exposed workers and included an immunologic component suggested a genetic basis for susceptibility to BeS and CBD. Anti-human leukocyte antigen (HLA) antibodies blocked the Bespecific proliferation response in lymphocytes of workers sensitized to Be (Saltini et al. 1989). Subsequently, a molecular epidemiologic study implicated HLA-DPB1 encoding a glutamic acid in position 69 of the mature protein with susceptibility to CBD (Richeldi et al. 1993).
Recent studies have confirmed that genetic susceptibility to CBD involves HLA glycoproteins, T-cell receptor clonality, and have suggested the involvement of tumor necrosis factor polymorphisms (Tinkle et al. 2003). These HLA glycoproteins are heterodimers comprising an αchain and a β-chain, which form an antigen-binding moiety. Binding of the antigen occurs in an elongated groove on the surface of this complex (Figure 1). The floor of the groove is composed of β-sheets, and the walls of α-helices (Li et al. 2000).
To investigate the basis for the suggested interindividual genetic susceptibility to CBD, we constructed three-dimensional structural models of HLA-DP representing several distinct allotypes and calculated the electrostatic potential on the surface of these models. The allotypes were predicated upon two types of α-chain sequences and five types of β-chain sequences. The α-chain sequences encoded by HLA-DPA1*0103 and *0201 alleles, and β-chains encoded by HLA-DPB1*1701, *0201, *1901, *1301, and *0401 alleles were selected, allowing for 10 possible analytic models.

Molecular Modeling
There are no experimentally determined coordinates for HLA-DP in the Protein Data Bank (PDB; http://www.rcsb.org/ pdb/). We modeled the extracellular part of HLA-DP by homology to a known HLA species. Overall, the primary, secondary, and tertiary structure in the class II family of major histocompatability complex (MHC) is well conserved. We have selected one HLA-DR, PDB entry 1FV1 (HLA-DRA*0101/B5*0101; Li et al. 2000) as a template. Our choice was based on a relatively good resolution of the template compared with other HLA entries (1.90 Å, R-factor = 0.233), the total number of crystallographically resolved residues, and the results of a sequence alignment using the GCG program (Accelrys, Inc., San Diego, CA).
The SYBYL 6.6 molecular modeling software (Tripos, Inc., St. Louis, MO) was used to make modifications to the HLA-DR template. It was first necessary to repair several residues in the template, in which coordinates of some side-chain (SC) atoms were missing. Next, specific SC substitutions were introduced, which converted the sequence of HLA-DR protein to that of HLA-DP. The deletion of two residues in the gap required reconstruction of one loop region. After addition of all hydrogen atoms, positions of the selected SC atoms and the backbone (BB) of the engineered loop were subjected to 50 iterations of the steepest descent minimization, followed by an exhaustive minimization using the Powell method (Powell 1977), until the full convergence by the energy gradient was reached. For minimization, the Kollman all-atom force field (Weiner et al. 1984) with distance-dependent dielectric and a cutoff radius of 10 Å was used. All BB and C β (except proline C β ) atoms were held fixed by treating them as a rigid body aggregate. Stereochemical quality of the HLA-DP model was tested using the WHATIF program (Vriend 1990). Modeled coordinate sets are available from the authors upon request.

Electrostatic Calculations
The electrostatic potential was calculated using the DelPhi module of the Insight II program (Accelrys). Insight II was also used to generate the isopotential surface maps. The electrostatic potential was calculated by numerically solving the linearized Poisson-Boltzmann equation using a finite difference method. The low-dielectric macromolecule (protein interior) was embedded in a high-dielectric continuum environment (water exterior). A solution with charged ions was simulated having an assigned ionic strength of 0.145, which typifies conditions at the physiologic pH. Using an ionic radius of 2.0 Å for the solvated ions, a Stern layer was constructed around the solvent-accessible surface, having a null ionic strength inside, which determines the maximum distance that an ion can approach the solvent-accessible surface. The assigned temperature was 298 K. The resulting system was discretized on a grid, and the potential at the grid points was solved iteratively starting from the Debye-Hückel boundary conditions. To improve the accuracy of calculated potentials, we used a method of grid focusing. This involved three additional calculations in which the molecule was allowed to occupy a successively larger fraction of the cubic grid volume so that values of the potential at the grid boundary could be calculated using the larger grid from each preceding calculation. The grid dimensions were selected to be 1,200 Å (spacing 4 Å/grid point), 600 Å (spacing 2 Å/grid point), 300 Å (spacing 1 Å/grid point), and 100 Å (spacing 0.75 Å/grid point). The interior and exterior dielectric constants were 4 and 80, respectively. The solute boundary was defined by a solvent-accessible surface generated by a rolling probe sphere of 1.4 Å radius. Partial atomic charges and radii were taken from the PARSE (Sitkoff et al. 1994) parameter set, which was originally optimized for proteins by fitting the solvation free energies of amino acid analogs.

Estimation of the Odds of Adverse Health Outcome Associated with Specific HLA-DPB1 Haplotypes in Beryllium Workers
It was suggested that susceptibility to CBD conferred by different HLA-DPB1*E69 haplotypes is not equal. In Be-exposed workers inheritance of certain HLA-DPB1*E69 haplotypes might convey greater susceptibility than others (Wang et al. 1999). Initially, we examined the data they reported and attempted to deduce allele-specific risk estimates (McCanlies et al. 2003;Tinkle et al. 2003). Because the number of cases was small (n = 20), we further pooled these data with data from several more-recent studies (Rossman et al. 2002;Saltini et al. 2001;Wang et al. 2001). A simple allelic hierarchy was developed based on the proportional allelic frequencies reported for CBD, Be-sensitized, and controls. Odds ratios (ORs) were calculated using the Mantel-Haenzel estimate (SAS Institute, Inc. 1999) of the common relative risk for each allele (HLA-DPB1*0201, *0401, *0601, *0901, *1001, *1301, *1601, *1701, and *1901) by analyzing each study population individually and by pooling data from all studies.

Results
In this study we investigated a putative link between electrostatic properties of HLA-DP allotypes and corresponding ORs  of epidemiologically observed susceptibility to CBD. The molecular electrostatic potential (MEP) represents the electrostatic properties of a molecule, quantummechanical effects notwithstanding. Stretching around the molecule, the MEP comprehensively describes intermolecular electrostatic interactions, which often have a significant impact on ligand binding. The electrostatic potential around the extracellular part of the HLA-DP protein was calculated using computer-generated homology models of HLA-DP allotypes. The calculations were performed in the absence of antigenic peptides, which for HLA-DP are not yet known. Further, these calculations did not take into account the phospholipid membrane responsible for anchoring the protein. Because the antigen-binding groove and the transmembrane domain of class II MHC are located a relatively large distance apart (approximately 60 Å) at the opposite ends of the extracellular part of the complex, the effect of the membrane on the electrostatic potential in the vicinity of the binding groove was assumed to be similar in different HLA-DP allotypes. Thus, the present analysis focuses on the distinctions in electrostatic potential on different HLA-DP allotypes.
The homology models of HLA-DP allotypes were generated from the X-ray coordinate set of the HLA-DR protein, Brookhaven code 1FV1. The amino acid sequence alignments between the corresponding chains of 1FV1 and HLA-DP are shown in Figure 2. The sequence alignment yielded approximately 64% identity in a 181 amino acid overlap between 1FV1A and the HLA-DPA1 sequences, and 65% identity in a 190 amino acid residue overlap between 1FV1B and the HLA-DPB1 sequences. There was a single two-residue gap in the alignment of β-chains, beginning at position 24, and there were no gaps at all in the alignment of the α-chains. Overall, the results of sequence analysis indicate a high homology between the three-dimensional structure of 1FV1 and HLA-DP. Usually, the high structural homology between two protein species implies an unequivocal reconstruction of unknown species from the coordinates of a known template.
Further validation of the derived structural models was conducted using the WHATIF program (Vriend 1990), and the coordinates of homology models were subjected to a number of quality control tests. To predict the probabilistic likelihood for the particular rotamer of a given SC in the model, the distribution of the χ 1 torsion angles of the model was compared with a distribution of rotamers in the WHATIF database. Scores between 0 and 1 were computed, in which a value close to 1 implies high likelihood for that rotamer. A histogram of the scores collected for all residues in both the αand β-chains is presented in Figure 3. In general, for most crystal structures subjected to this test, less than 10% of residues have scores lower than 0.4. Only three residues in the α-chain scored below 0.4, and there were no scores in any of the chains lower than 0.3. These results agree with expectations arising from a comprehensive analysis of Toxicogenomics | Electrostatic potential on HLA-DP   An additional WHATIF check involved looking at the fine-packing quality for both the homology model for HLA-DPA1*0103/B1*1701 and the 1FV1 protein. Quality control values for individual residues are intended to measure the fit of a residue in the particular part of the structure. If the z-score for the protein is greater than -2.5 (2.5 standard deviations less than the mean), then the structure is acceptable or good. A z-score less than -3 suggests that the structure likely is poor. For all BB-SC contacts, the z-scores were -2.78 for the 1FV1 protein, and -2.87 for the model structure. For all BB-BB contacts, the z-scores were 0.38 and 0.62 for 1FV1 and model, respectively. It appears that reconstruction of the loop region in the model produced an improvement in quality with respect to BB-BB contacts, as this was the only modification made to the BB of the 1FV1 template. The z-scores for all types of contacts (BB-BB, BB-SC, SC-BB, and SC-SC) were -1.52 and -1.27 for the 1FV1 and model, respectively, indicating an overall improvement in the packing quality for the homology model compared with the template.
As a final WHATIF check, a torsion angle evaluation was performed. The computed scores indicate, for each residue, how normal the torsion angles are, as described by normality scores using all torsion angles except ω. Average values and standard deviations are extracted from residues in the WHATIF database and used to compute zscores. A residue with a z-score less than -2.0 is poor, indicating that more than one torsional angle is in a highly unlikely position. For 1FV1, a total of 5 residues (3 in the α-chain and 2 in the β-chain) were listed as poor. For the HLA-DPA1*0103/ B1*1701 model, a total of seven residuals were listed as suspicious (2 in the α-chain and 5 in the β-chain). Residues having forbidden ϕ-ψ combinations are listed, along with residuals with unusual ω angles, which are angles that deviate by more than 3σ from the normal value. It is typical to find that approximately 5% of the residues in a structure have unusual ϕ-ψ combinations. There were 17 residues (out of 369 in the structure) listed with poor scores for the 1FV1 protein, and 22 residues (out of 367 in the structure) listed for the model. The χ 1 /χ 2 correlation z-score indicates the quality of correlation between χ 1 and χ 2 angles. Both values are considered acceptable. Thus, all quality control tests performed on the structures indicate a high likelihood of the derived homology model.
The initial assessment of electrostaticsrelated effects was carried out using the primary structures of HLA-DP allotypes. Figure 4 represents the amino acid sequence alignments of commonly found HLA-DP alleles. Epidemiologic studies have not provided convincing evidence for an association between and susceptibility to CBD of α-chain alleles of HLA-DP ( Figure 4A). Nevertheless, they are required for reconstruction of the three-dimensional structure of the extracellular part of MHC class II complex. Some β-chain alleles have also been surveyed. Figure 4B represents the β-chain alleles for which epidemiologic association data were available at the time of the study.
Except for the residue α50 ( Figure 4A), all substitutions involving ionizable residues are confined in a relatively short stretch of the β-chain alignment ( Figure 4B). These are residues β55-β56, β69, and β84-85. From the analysis of the three-dimensional model, it appears that all residues involved in the substitutions affecting the charge are spatially adjacent to the groove (Figure 1). Table 1 lists the total charge on the βchains encoded by the HLA-DPB1 alleles. The variation in total charge between the different alleles is determined by the presence of acidic (-1) or basic (+1) substitutions in the three aforementioned regions. For example, comparing *0401 and *1701 allotypes, the substitution of Asp β55, Glu β56, Asp β84, and Glu β85 for Ala β55, Ala β56, Gly β84, and Gly β85 produces a difference of -4e, and substitution of Glu β69 for Lys β69 contributes a -2e difference, making *1701 six units more negative than *0401. HLA-DPB1*1701, *1601, *1001, *0901, and *0601 contain similarly charged residues in the three regions. Thus, all these β-chains carry the same total charge, and the results for other alleles belonging to this group are expected to be similar. Only one of these β-chains, *1701, was used in this study.
In Table 1 the β-chain alleles are partitioned in four groups based on the total charge on the chain (column 2). Grouping β-chain alleles by charge in this way provides for a more robust statistical analysis. The ORs for either CBD or BeS  are reported in columns 6 and 7, respectively. The highest estimates of ORs for CBD were found consistently in β-chain alleles coding for isotypes that carry the greatest negative charge. Conversely, the lowest ORs for CBD were associated with alleles coding for the least negatively charged HLA-DPB1*0401 (K69) isotype. HLA-DPB1*1901 and *0201 and *1301 carry intermediate negative charge, and they also were associated with intermediate ORs for CBD. Overall, there seems to be an obvious link between the amount of negative charge on the β-chain and corresponding ORs for CBD. The highest ORs for BeS were also found in β-chain alleles that code for isotypes carrying the greatest negative charge. Comparison of these ORs with ORs associated with less negatively charged isotypes also indicates a likely possibility of connection between the charge and ORs. However, in the case of BeS, the picture is masked by the HLA-DPB1*1901 allele, whose OR disrupts the otherwise regular rank order in column 7. At present we do not have a comprehensive explanation for this observation, except that it may be the result of a small number of observations. The epidemiologic basis for the calculated *1901 OR in BeS is one reported case; this is reflected in large, overlapping confidence intervals (CIs) calculated for this allele. In addition, the deviation may be caused by interference of unknown α-chains (see below), which have not been reported in the source studies.
It is evident from the table that the ORs have an approximately log-linear relationship to the magnitude of negative charge on the β-chain of HLA-DP. Figure 5 shows fits of ungrouped log-ORs using the log-linear regression. The variance-covariance matrix analysis suggests that log-ORs for both CBD and BeS are reasonably well correlated with the charge. The correlation coefficients are -0.85 and -0.67, respectively. As the data are naturally stratified by charge, we further broke the matrix into lack-of-fit and pure-error terms. The lack-of-fit term represents the deviation of fitted values from group means, i.e., it characterizes the quality of the fit. The pure-error term specifies a scatter of the data around the means, which is the statistical noise that interferes with the regression. The lack-of-fit test indicated a statistically representative fitting with loglinear regression of both the CBD and BeS data. The calculated test-statistic numbers F 2,5 were 0.16 and 1.04, respectively, indicating that the null hypothesis for a poor fit could not be accepted even at the 60% level. A fairly large scatter of the data within the groups, which is seen in Figure 5, is due to the limited data available for individual alleles. The coefficients of determination, which measure the random scatter, were 0.71 and 0.45 for the CBD and BeS regressions, respectively. Despite this, both regressions were statistically significant according to the standard test for the zero slope of regression line (F 1,7 are 17.51 and 5.68, respectively). Thus, we conclude that the analyzed epidemiologic association data support the hypothesis of an association between the charge on an allele and observed ORs for both CBD and BeS.
To more effectively analyze how the charge distribution varies among HLA-DP proteins, we calculated the electrostatic potential on the solvent-accessible surface of each protein. Figure 6 illustrates the locations of the three regions, which according to Table 1 determine the variation of charge on the β-chain. Each of the three regions is located in proximity to the peptide-binding groove. Figure 6A and B compare the electrostatic potential at the solvent-accessible surface of HLA-DP allotypes involving the least and one of the most negatively charged β-chain alleles. They are HLA-DPA1*0103/B1*0401 and HLA-DPA1*0103/B1*1701, respectively.
These maps clearly show intensification of negative potential in the binding-groove region of the protein encoded by *1701 allele compared with that encoded by *0401. In addition to the solvent-accessible potential surface maps, which show how the potential varies along the surface drafted by the centers of the rolling probe sphere (e.g., the ion), it is also useful to analyze the MEP on the protein, which is given by the isopotential surface maps. Figure 7 illustrates the MEP at a contour level of -2 kT/e (red) and +2 kT/e (blue) for HLA-DPA1*0103/B1*0401, *0201, *1301, and *1701, respectively. The MEPs on HLA-DPA1*0103/B1*1301, *1901 and *1701 are compared in Figure 8, in which the view is rotated by 90°to analyze the effect of the substitution in position β56. Variations of MEP in other parts of the protein are fairly minor, as the substitutions of charged residue are confined to the binding groove only.
Comparison of MEPs in Figures 6-8 shows how the potential surface is altered by the amino acid substitutions, and hence, the charge in the regions β55-β56, β69, and β84-β85. In the case of the *0401 protein ( Figure 7A), in which the β-chain does not carry any negatively charged residues in positions β55, β56, β84, and Toxicogenomics | Electrostatic potential on HLA-DP Environmental Health Perspectives • VOLUME 111 | NUMBER 15 | November 2003 CI, confidence interval. a The charge on each allotype was determined under an assumption of standard pK a values for all ionizable residues. It implies that at the physiological pH, the residues adopt fully ionized states, i.e., aspartate and glutamate each carry a charge of -1e, and arginine and lysine carry a charge of +1e. Histidine and cysteine residues were assumed to be neutral. b ORs calculated from pooled data reported by Rossman et al. (2002), Saltini et al. (2001), and Wang et al. (1999Wang et al. ( , 2001. c Could not be deduced from Wang et al. (1999Wang et al. ( , 2001. β85, there is virtually no appearance of negative MEP in the peptide-binding groove region. Instead, the central part of the groove is occupied by a stretch of positive MEP, which is due to a positively charged Lys β69. Figure 7B illustrates how substitutions for negatively charged residues in positions β55, β56, and β69 induce negative MEP in the binding groove of the *0201 protein. Similarly, Figure 7C shows the emergence of negative MEP resulting from substitutions for negative residues in the positions β69, β84, and β85 of the *1301 protein. Comparison of Figure 7B and C indicates how much Glu β69 contributes to the negative MEP on allotypes involving the *0201 and *1301 alleles. Figure 7D demonstrates a situation in which negatively charged residues are substituted for all residues constituting the three genetic marker regions. When Figure 7B-D is compared, we can see that the substitutions in the three regions result in a roughly additive increase of the surface of negative MEP in and around the groove. In the future, in-depth quantitative analysis will be applied to determine the magnitude of negative MEP inside the groove. However, the shape of negative MEP around Glu β69 in Figure 7D is apparently enlarged compared with Figure 7B and C, probably originating from the increased negativity of MEP inside the groove. If this is the case, it may partially explain the link between charge and ORs reported in Table 1 and Figure 5. The effects of replacing Ala β55 in the *1301 protein with Glu β55 in the *1901 protein and replacing both Ala β55 and Ala β56 with negatively charged residues in the *1701 protein are illustrated in Figure 8. Similar to Figure 7, substitution for charged residues results in a corresponding appearance of increasingly more negative MEP. All studied HLA-DP β-chain proteins contain either an Asp or Glu residue in the position β57, which itself influences the shape and magnitude of MEP in this region. For example, there is a Glu in position β57 in the *0201 protein, whereas the *1701 protein has an Asp in this position. Evidently, the presence of Glu or Asp in this position affects the appearance of the negative MEP around the region of β56, β57, as can be seen by examining Figure 7B and D (or Figure 9A and B, described below).
The total charge on the α-chain encoded by HLA-DPA1*0201 is -17, and it is -18 on HLA-DPA1*0103. The replacement of Gln50 in the α-chain encoded by HLA-DPA1*0103 with Arg50 encoded by HLA-DPA1*0201 imparts one additional unit of positive charge. The MEP on two HLA-DP allotypes involving the *0201 allele of the α-chain is shown in Figure 9. The figure represents homology models of HLA-DPA1*0201/HLA-DPB1*1701 and HLA-DPA1*0201/HLA-DPB1*0201. These should be compared with Figures 7B  and 7D, which are for HLA-DPA1*0103/ HLA-DPB1*1701 and HLA-DPA1*0103/ HLA-DPB1*0201, respectively. Comparison of these two pairs shows that the effect of adding one additional unit of positive charge on the α-chain can be deduced in the absence (Figures 7B and 9A) and presence ( Figures 7D and 9B) of negative charge in the β84-β85 region of the β-chain. The Q50R substitution in the α-chain causes a profound effect on the MEP in the immediate vicinity of α50. However, it is not clear how much this substitution affects the MEP inside the groove. Conversely, unlike β69, this residue is located at the very edge of the peptide-binding groove, and thus its effect on the binding may be fairly small.  The surface is a solvent-accessible surface generated as described in the text. The wire frame denotes a putative location of the antigen as it is given by the 1FV1 template. Red and blue represent contours at -2 and +2 kT/e, respectively. The isopotential surface map is superimposed on the solvent-accessible surface (shown in white), so that only a part of the isopotential surface that lies on or above the solvent-accessible surface is visible. The wire frame denotes a putative location of the antigen as it is given by the 1FV1 template.
Alternatively, the other two important markers of the β-chain, β55-β56 and β84-β85, are similarly distant from the center of the groove as α50. From Figures  8 and 9, we already know that the effect of substituting different residues in one particular region can, to some extent, alter MEP in an adjacent region. Additional epidemiologic studies may help to clarify the role of α-chain substitutions in the etiology of CBD. From a structural point of view we do not see obstacles to interference with binding from yet unknown substitutions in the α-chain. Alternatively, if such substitutions indeed do not affect ligand binding, additional computational investigations are needed to explain the physiologic differences between substitutions in the αand β-chains.

Discussion
Occupational or environmental exposures to metals may cause severe health disorders related to modulation of immune homeostasis. Metal immunomodulation has been implicated in chronic inflammatory processes, autoimmune diseases, and related adverse effects (Cunningham-Rundles et al. 2000;Lawrence and McCabe 2002). Molecular mechanisms of immunological actions of metals are largely unknown. In this work we explored a possible link between molecular properties of genetic products conferring susceptibility to Be and the health effects of Be exposures. The results of this work are useful for interpreting the molecular recognition in MHC class II proteins and for investigating their interactions with Be 2+ . Previous studies have suggested Glu β69 as an important genetic marker for conferring susceptibility to CBD (Díaz et al. 1998;Fontenot et al. 2000;Lombardi et al. 2001;Potolicchio et al. 1999;Richeldi et al. 1993Richeldi et al. , 1997Rossman et al. 2002;Saltini et al. 2001;Wang et al. 1999Wang et al. , 2001. More recent studies have led some investigators to suggest that the charged residues Asp β55 and Glu β56 might also contribute to the recognition of HLA-DP by T cells (Díaz et al. 1998;Fontenot et al. 2000;Potolicchio et al. 1999). In a study of the binding of cobalt to HLA-DP, Potolicchio et al. (1999) compared the binding affinity for cobalt by HLA-DPB1*0201 (which has Asp β55, Glu β56, and Glu β69) and by HLA-DPB1*0402 (which has Asp β55, Glu β56, and Lys β69), and also by HLA-DPB1*0401 (which has Ala β55, Ala β56, and Lys β69). They found that the binding affinity for radioactive 57 Co 2+ decreased in the order HLA-DPB1*0201 > HLA-DPB1*0402 > HLA-DPB1*0401, confirming earlier evidence for involvement of Glu β69 but also suggesting that positions β55 and β56 may also be involved in the interaction. A similar study with HLA-DPB1*0101 and HLA-DPB1*0401, which both carry Ala β55, Ala β56, and Lys β69, failed to present Be 2+ to CD4 + T cells, whereas HLA-DPB1*0202 and other alleles with Glu β55, Asp β56, and Glu β69 were positive for presentation of Be 2+ to T cells (Fontenot et al. 2000).
In this study the electrostatic potential on model proteins encoded by several HLA-DPA1 and HLA-DPB1 alleles was analyzed. The genetic polymorphisms in different alleles are responsible for substitutions of charged residues (Glu, Asp, Arg, Lys) in positions 55-56, 69, and 84-85 of the amino acid sequence of the β-chain (Figure 1). We found a clear correlation between the total charge on protein variants and epidemiologic odds ratios conferring susceptibility to either BeS or CBD. The MEP contours in Figures 7 and 8 illustrate the distribution of charge in the three significant regions in proximity to the binding groove. Hence, HLA-DPA1*0103/B1*1701, which carries the most negative charge, also appears to confer the greatest susceptibility for BeS and CBD. We appreciate that in all cases of HLA-DPB1 alleles coding for glutamic acid in position 69, the 95% confidence intervals overlap. This stems from the lack of sufficient power of the studies, both independently and together, to answer the question of allele-specific risk hierarchy. However, given the consistency of the different studies, the data strongly suggest that such a hierarchy exists.
Originally, we examined data reported by Wang et al. (1999) and were able, to a limited extent, to account for homozygosity. The analyses presented here do not account for the impact of homozygosity versus heterozygosity. Moreover, we are unaware of the impact of any potential allele misclassification, which is suggested by the growing numbers of haplotypes reported since 1999. Regarding the ranking of alleles with respect to BeS (Table 1), we considered the possibility of bias, potentially introduced by the data reported by Saltini et al. (2001), as, unlike those of Wang et al. (2001) and Rossman et al. (2002), these data did not show an association between inheritance of Toxicogenomics | Electrostatic potential on HLA-DP Environmental Health Perspectives • VOLUME 111 | HLA-DPB1*E69 and BeS. However, elimination of the BeS data reported by Saltini et al. (2001) from our analyses did not change the overall ranking. All of the reports described were case-control studies, and we remain ignorant of many facets of the natural history of CBD, especially as it pertains to the progression from BeS to disease. In addition, diagnosis of BeS remains problematic. It may be that a relatively high OR associated with CBD but a relatively low allele-specific OR associated with BeS reflects a greater propensity for progression from BeS to CBD (e.g., as may be the case for HLA-DPB1*1001).
As Figure 5 illustrates, there appears to be a log-linear relationship between ORs and charge on the protein. Although an improved set of data may alter the fit, it is noteworthy that such a relationship exists, as it suggests that this kind of analysis may be useful for eventually predicting disease susceptibility for new allotypes for which epidemiologic data are lacking or unavailable.
As far as the molecular mechanism of CBD is concerned, a previous study has indicated that Be 2+ may bind in the vicinity of Glu β69, perhaps directly in the peptidebinding groove . In this connection the log-linear dependence of the odds ratios (essentially the probability expressed in the units of background) on the charge is curious because it resembles the Boltzmann distribution in a two-site charge-charge interaction model when all other parameters are fixed. At present, we are unaware whether it is a coincidence or fact reflecting a yet unknown peculiar mode of ion binding. This puzzle may be resolved in future studies, which will shed additional light on the molecular mechanisms of metal-dependent immunomodulation.

Conclusion
Our results complement recent findings, which implicate positions β55, β56, β69, β84, and β85 of HLA-DPB1 in susceptibility to CBD. This approach may prove useful in predicting risk associated with previously unknown alleles and may help lead to the elucidation of mechanistic issues associated with CBD. Further studies will be directed at identifying the beryllium binding site on the HLA-DP protein.