HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Nielsen, J. E. Articles by McCammon, J. A. Search for Related Content PubMed PubMed Citation Articles by Nielsen, J. E. Articles by McCammon, J. A. Social Bookmarking                   What's this?

On the evaluation and optimization of protein X-ray structures for pKa calculations

Howard Hughes Medical Institute and Department of Chemistry and Biochemistry, Department of Pharmacology, University of California, San Diego, La Jolla, California 92093, USA

Reprint requests to: Jens Erik Nielsen, Department of Chemistry and Biochemistry M/C 0365, UCSD, La Jolla, CA 92093, USA; e-mail: jnielsen{at}mccammon.ucsd.edu; fax: (858) 534-7042.

(RECEIVED August 27, 2002; FINAL REVISION October 31, 2002; ACCEPTED November 4, 2002)

Article and publication are at http://www.proteinscience.org/cgi/doi/10.1110/ps.0229903.

	Abstract

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
The calculation of the physical properties of a protein from its X-ray structure is of importance in virtually every aspect of modern biology. Although computational algorithms have been developed for calculating everything from the dynamics of a protein to its binding specificity, only limited information is available on the ability of these methods to give accurate results when used with a particular X-ray structure. We examine the ability of a pKa calculation algorithm to predict the proton-donating residue in the catalytic mechanism of hen egg white lysozyme. We examine the correlation between the ability of the pKa calculation method to obtain the correct result and the overall characteristics of 41 X-ray structures such as crystallization conditions, resolution, and the output of structure validation software. We furthermore examine the ability of energy minimizations (EM), molecular dynamics (MD) simulations, and structure-perturbation methods to optimize the X-ray structures such that these give correct results with the pKa calculation algorithm. We propose a set of criteria for identifying the proton donor in a catalytic mechanism, and demonstrate that the application of these criteria give highly accurate prediction results when using unmodified X-ray structures. More specifically, we are able to successfully identify the proton donor in 85% of the X-ray structures when excluding structures with crystal contacts near the active site. Neither the use of the overall characteristics of the X-ray structures nor the optimization of the structure by EM, MD, or other methods improves the results of the pKa calculation algorithm. We discuss these results and their implications for the design of structure-based energy calculation algorithms in general.

Keywords: pKa calculations; crystal contacts; structural genomics; molecular dynamics; electrostatics; enzymes

Abbreviations: HEWL, hen egg white lysozyme • RMSD, root mean square difference • PDB, Protein Data Bank

	Introduction

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
Obtaining the X-ray structure of a protein has become a standard requirement in modern biology, not too unlike identifying the location of its gene on the chromosome and setting up an efficient expression system. Once available, the X-ray structure provides a wealth of information on how to interpret past experiments and how to design new experiments that will provide information on the function of the protein. The X-ray structure by itself, however, does not reveal much information regarding the physical characteristics of the protein. Extensive calculations are needed to determine the substrate or ligand specificity, the stability, the dynamics, and the electrostatic features of the protein, and it is often not possible to know whether the results of such calculations are trustworthy. With the vast number of protein X-ray structures being solved in various structural genomics projects (Heinemann et al. 2000; Yokoyama et al. 2000; Stevens et al. 2001), it is becoming more and more important to have access to fast and reliable algorithms that can tell us something about the physical characteristics of a protein from its X-ray structure. Presently algorithms are available for predicting everything from the pKa values of a protein (Bashford and Karplus 1990; Yang et al. 1993; Antosiewicz et al. 1994, 1996; Demchuk and Wade 1996; Alexov and Gunner 1997; Sham et al. 1997, 1998; Alexov and Gunner 1999; Mehler and Guarnieri 1999; Nielsen and Vriend 2001) to its large-scale motions (Amadei et al. 1993; de Groot et al. 1997, 1999) and binding characteristics (Goodsell et al. 1996; Kramer et al. 1999). However, although most of these algorithms give a good correlation with experimental data for a subset of existing protein structures (Kramer et al. 1997, 1999; Mehler and Guarnieri 1999; Gabdoulline and Wade 2001; Nielsen and Vriend 2001; Guerois et al. 2002), it is not trivial to know when an algorithm gives a reliable result when used with a novel X-ray structure. This is mainly due to the high sensitivity to the details of the protein structure displayed by many structure-based algorithms that calculate energies inside proteins. It is well known, for example, that small molecule docking programs are highly dependent on having the "correct" structure of the protein/receptor (Claussen et al. 2001). In addition, the results of molecular dynamics (MD) simulations (Braxenthaler et al. 1997) and electrostatic calculations (Nielsen et al. 1999; Nielsen and Vriend 2001) are known to be sensitive to structural details.

When calculating a property of a protein from its X-ray structure, which is dependent on structural details, it is therefore essential to know whether the structure is capable of giving accurate results with the algorithm in question. Information on the usefulness of the X-ray structure can often be inferred from a visual inspection of the structure. For example, if a ligand is present in the active site, then it is likely that the structure is well suited for drug design and docking studies, whereas an extended or open structure, as seen for example for some of the protein kinases (Cox et al. 1994), indicates that the structure is poorly suited for studies of the active form of the enzyme.

Even much smaller changes in the structure of a protein are also likely to have a profound effect on the results of structure-based energy calculations, as illustrated by the large differences in the calculated pKa values of the hen egg white lysozyme (HEWL) active-site residues resulting from a 180° change in the ² angle of HEWL Asn 46 (Nielsen et al. 1999). Here we examine the sensitivity of a pKa calculation algorithm to the structural differences among 41 HEWL wild-type X-ray structures. This analysis provides us with information on the reliability of pKa calculations when used with a given X-ray structure. Because the desolvation energies and the electrostatic interaction energies that are calculated by the pKa calculation algorithm are essential components of most structure-based energy calculations, the conclusions that we present here are applicable to other types of structure-based energy calculation methods.

We investigate whether it is possible to select a more reliable subset of HEWL structures for pKa calculations by using properties of the X-ray structure such as the resolution, the crystallization conditions, and the output of structure validation software. We also investigate whether it is possible to standardize the HEWL structures by a computational protocol so that all give the correct result with the pKa calculation algorithm. Finally, we discuss the implications of the present results for the design and application of protein structure-based energy calculations in general.

	pKa calculation algorithms

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
pKa calculation algorithms are used mainly in studies of enzyme mechanisms (Raquet et al. 1997; Lamotte-Brasseur et al. 2000; Morikis et al. 2001b) and in the study of protein stability (Yang and Honig 1993, 1994; Lambeir et al. 2000; Morikis et al. 2001a). In the study of enzyme mechanisms, these algorithms aid by identifying the residues that are likely to be proton donors and proton acceptors, and in protein stability studies they are capable of predicting the origins of the pH-dependence of protein stability. Most pKa calculation algorithms rely on finite-difference solvers of the Poisson-Boltzmann equation (FDPB-solvers) to provide the electrostatic energies of a protein structure, although several alternative approaches to pKa calculations exist (Sham et al. 1997; Mehler and Guarnieri 1999; Sandberg and Edholm 1999). In the present paper we deal exclusively with an FDPB-based pKa calculation algorithm.

The major differences between FDPB-based pKa calculation algorithms lie in the way that they model protein flexibility. The treatment of the protein flexibility can be divided roughly into two classes: explicit treatment and implicit treatment. Methods that treat the flexibility of the protein explicitly employ MD simulations (Zhou and Vijayakumar 1997; van Vlijmen et al. 1998; Gorfe et al. 2002), proton optimization (Alexov and Gunner 1997), or rotamer optimization techniques (Alexov and Gunner 1999). Methods with implicit treatment of protein flexibility typically adjust the dielectric constant for the entire protein (Antosiewicz et al. 1994; Karshikoff 1995; Antosiewicz et al. 1996) to achieve better correlation with experimental results, although algorithms that use a residue-dependent value of the protein dielectric constant have also been developed (Demchuk and Wade 1996; Nielsen and Vriend 2001).

Generally, the methods that use an implicit description of the protein flexibility have been more successful in obtaining a good overall correlation with experimental data, whereas the methods that optimize the hydrogen-bond network have proven superior in calculating active-site pKa values (Nielsen and Vriend 2001), presumably because the details of the hydrogen-bond network are very important in active sites. It has also been reported that structural averaging (van Vlijmen et al. 1998; Gorfe et al. 2002) can improve the correlation between experimental and calculated pKa values, although the improvements in some cases seem to be insignificant (Koumanov et al. 2001). Attempts at incorporating pKa calculations in MD algorithms have also been made (Baptista et al. 1997), but presently these algorithms have not proven to give a significant improvement in the accuracy of the calculated pKa values.

	Calibrating pKa calculation methods

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
The accuracy of a pKa calculation result is evaluated by calculating the RMSD between the calculated and experimentally determined pKa values for all protein titratable groups that titrate within 5–6 pH units of physiological pH. Consequently, pKa calculation algorithms have been calibrated to give low RMSD values for a set of well-behaved model proteins with experimentally measured pKa values. However, most titratable groups are situated on the surface of proteins, and a significant fraction of surface groups are involved in crystal contacts (Carugo and Argos 1997; Valdar and Thornton 2001). Because crystal contacts perturb the details of the local structure, and restrict the mobility, they can induce the formation of salt bridges and charged hydrogen bonds in the crystal that are present only transiently in solution. The calibration of pKa calculation methods that use an implicit description of protein dynamics on a large unfiltered set of titratable residues is therefore bound to introduce a bias in the pKa calculation methods, such that the pKa values of surface residues are calculated correctly even though the conformations of these are different from the conformations that they occupy in solution (Nielsen and Vriend 2001).

It is tempting to speculate that this is a significant part of the reason why many pKa calculation methods give the best results with a relatively high protein dielectric constant (Antosiewicz et al. 1994; Demchuk and Wade 1996), which essentially smears out the effect of the surrounding protein environment. It is therefore our belief that more accurate pKa calculation methods can be constructed by calibrating pKa calculation methods on a set of experimentally measured pKa values which does not contain any titratable groups that are influenced by crystal contacts.

In the present study we chose to focus exclusively on the calculated pKa values for the two key active-site residues of HEWL. We did so because one of the most important uses for pKa calculation algorithms is to identify the proton donor from a set of titratable residues in the active site of an enzyme (Raquet et al. 1997; Lamotte-Brasseur et al. 1999, 2000), and the pKa values of surface residues are not essential when answering such a question.

	Identifying the proton donor in a catalytic mechanism

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
Many enzymes have bell-shaped pH activity profiles, and this naturally leads to the assumption that catalysis at low pH is limited by the protonation of an active-site residue, and similarly that catalysis at high pH is limited by the deprotonation of another active-site residue. Enzymatic pH activity profiles can generally be decomposed into a pH-k_cat profile and a pH-K_m profile. From the pH-k_cat profile, one can extract the pKa values of the active-site groups when the substrate is bound, whereas the pH-k_cat/pH-K_m profile will give the pKa values of these two groups in the apo-form of the enzyme (Kyte 1995). Here we will assume that the same two groups are responsible for the shapes of both the pH-k_cat profile and the pH-k_cat/pH-K_m profile, and, furthermore, that these two groups are the catalytic nucleophile (the group that limits activity at low pH) and the proton donor (the group that limits activity at high pH) in the catalytic mechanism. It is not obvious that this is the case for most enzymes because experimental data on the subject are very scarce. However, for Bacillus circulans xylanase, convincing experimental data (McIntosh et al. 1996; Joshi et al. 2001) have been presented that justify the assumption that the catalytic nucleophile and the proton donor are indeed the residues that govern the shape of the pH-k_cat/pH-K_m profile. Because the catalytic mechanism of the xylanases is identical to that of HEWL, and because it is well known that the proton donor in HEWL (Glu 35) indeed has an elevated pKa value in the apo-form of the enzyme (Demchuk and Wade 1996), we will identify the proton donor in the HEWL catalytic mechanism by examining pKa values calculated from apo-crystal structures of HEWL.

Ideally, one should identify the proton donor in a catalytic mechanism as the residue which has a pKa value identical to the pKa value for the proton donor determined from kinetic data. Unfortunately, pKa calculation methods are not yet accurate enough to match kinetically measured pKa values directly, and a better strategy is therefore to identify the proton donor as the acidic group in the active site predicted to have the highest pKa value. In many cases, the choice is between two or three acidic residues, such as is the case for lysozyme and most other glycosyl hydrolases (Davies and Henrissat 1995), and in the present work we propose a set of criteria that, to the best of our judgement, enables us to confidently identify the proton donor in a catalytic mechanism from calculated pKa values for two acidic residues.

We require that the proton donor has a pKa value of at least 5.0, and that the difference between the pKa value of the proton donor and that of the other acid is at least 1.5 units, with the proton donor having the higher pKa value of the two. In the following we refer to these criteria as the "local identification criteria" or local ID criteria.

	Hen egg white lysozyme

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
HEWL is a 129-residue enzyme which serves as one of the paradigms for investigating the effect of crystallization conditions. The PDB contains more than 100 structures of wild-type HEWL, and HEWL thus provides an excellent model system for studying the effect of structural variation on the results of pKa calculation methods. HEWL is a monomeric single-domain enzyme, which consists of an all- region and a ß-rich region. The active site is situated in a cleft between the two regions, and the two key active-site residues are Glu 35 and Asp 52 (Fig. 1). HEWL is a retaining glycosyl hydrolase (Family 22 in the CaZy database http://afmb.cnrs-mrs.fr/CAZY/; Coutinho and Henrissat 1999), and recently it was elegantly proven that hydrolysis proceeds via a covalent enzyme-substrate intermediate (Vocadlo et al. 2001) with Glu 35 being the proton donor and Asp 52 the nucleophile in the catalytic mechanism. The initial step of the catalytic mechanism (Fig. 2) is the donation of a proton by Glu 35 to the glycosidic oxygen of the substrate. Subsequently, Asp 52 performs a nucleophilic attack on the anomeric carbon atom of the substrate, thus forming a covalent bond with the substrate. In the final step, the covalent enzyme-substrate intermediate is hydrolyzed by a water molecule, and the initial protonation states are regenerated.

View larger version (149K): [in this window] [in a new window]   Figure 1. Hen egg white lysozyme (PDB Id: 7LYZ). The active-site residues Glu 35 and Asp 52 are shown in red and cyan, respectively. The figure was prepared with the MOLSCRIPT (Kraulis 1998) and Raster3D (Merrit and Bacon 1997) programs.

View larger version (9K): [in this window] [in a new window]   Figure 2. The general catalytic mechanism for retaining glycosyl hydrolases. (I) Protonation of the glycosidic oxygen by the proton donor (Glu 35) and attack on the glucose C1 by the nucleophile (Asp 52). Departure of the reducing end of the substrate. (II) Activation of a water molecule, cleavage of C1-Asp 52 covalent bond. (III) Regeneration of the initial protonation states.

	Objective

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

	Results

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
We examined the ability of a protein pKa calculation algorithm to identify the proton donor in the catalytic mechanism of HEWL. Initially we evaluated the feasibility of identifying the proton donor using a large number of unmodified X-ray structures. We then continued to investigate the possibility of using crystallization criteria, the resolution of the X-ray structure, and structure validation software to select a better subset of HEWL X-ray structures, and we examined the reasons for correctly and incorrectly calculated pKa values using three representative HEWL structures as an example. Next we explored several ways of "correcting" HEWL X-ray structure coordinates to improve the correlation with experimental data, and finally we discuss the implications for protein structure-based energy calculations in the light of the results presented here.

Calculating pKa values for 41 HEWL X-ray structures
pKa values were calculated for all residues of 41 HEWL wild-type X-ray structures to identify the proton donor in the catalytic mechanism (Table 1). None of the structures contain any inhibitors or substrate molecules in the active site, and the structures therefore present a set of X-ray structures that could be obtained for the apo-form of any given enzyme. RMSD values for C positions between 2LZT and the rest of the structures are low (maximum RMSD, 1.62; mean value, 0.66), thus demonstrating that the structures are indeed very similar as measured by C-RMSD values.

Resolution
Low resolution can be responsible for incorrect positioning of atoms in X-ray structures, and we therefore expect a correlation between the resolution of the X-ray structure used for the pKa calculation and the accuracy of the pKa calculation results. Figure 3 shows the correlation between the percentage of structures for which the local ID criteria correctly identify Glu 35 and the cutoff for the resolution of the structures. For the set of X-ray structures presented here, there is no strong evidence of a correlation between the resolution and the quality of the pKa calculation, although the two structures with a resolution less than 1 Å both give a correct result.

View larger version (19K): [in this window] [in a new window]   Figure 3. Percentage of structures where the pKa calculations correctly identify Glu 35 as the proton donor () or identify Asp 52 () as judged by the local ID method. () The total number of structures within the cutoff. Only structures that have a resolution that is equal to or less than the cutoff value (x-axis) are considered.

Similarly there is no correlation between the temperature and the pKa calculation results, as evidenced by comparing the results for 1LSA, 1LSB, 1LSC, 1LSD, 1LSE, and 1LYS.

Space groups
The present set of HEWL structures are crystallized in five different space groups: P 43 21 2: 20 structures, P 21 21 21: 6 structures, P 1 21 1: 3 structures, P 21: 7 structures, and P 1: 5 structures. Only structures in the P 43 21 2 space group give rise to the identification of Asp 52 as the proton donor, and additionally seven more structures in this space group give pKa values that are inconclusive (i.e., neither Glu 35 nor Asp 52 can be identified according to the local ID criteria). In the four other space groups (21 of the 41 structures), only three structures give an inconclusive result, and in all remaining cases, Glu 35 is identified as the proton donor. In all but one of the P 43 21 2 structures, Asn 44 forms a crystal contact with either Arg 45 or Arg 68 from a symmetry-related molecule. Asn 44 forms a hydrogen-bond with Asp 52 in several structures, and is thus of critical importance for the protonation state of the active-site residues, as will be illustrated later.

Structure validation tools
It is possible to get a correct prediction of the catalytic proton donor for almost three-quarters of the unmodified HEWL X-ray structures, and although this number seems encouraging it also means that for one-quarter of all crystal structures, we are likely to get an inconclusive (or even wrong) answer when we apply a structure-based energy calculation method to an X-ray structure. Wrong and inconclusive answers do not present a major obstacle in themselves; the real problem is that it is not possible to distinguish "bad" from good results based only on the resolution and the crystallization conditions of the X-ray structure. Several tools have been constructed for the validation of protein X-ray structures. The more well known of these tools are WHAT_CHECK (Hooft et al. 1996b) and PROCHECK (Laskowski et al. 1993). We used WHAT_CHECK with all of the 41 HEWL X-structures, and analyzed the correlation of WHAT_CHECK Z-scores with the ability of the local ID criteria to correctly identify the proton donor in the catalytic mechanism. Figure 4 shows the percentage of correct and wrong predictions by the local ID criteria versus the Z-score cutoff. Only structures with WHAT_CHECK Z-scores above or equal to the cutoff were included in the analysis using the local ID criteria. Figure 4 shows a weak correlation between increasing Z-score and the percentage of structures that give correct results, but the significance of this correlation is too low to warrant any conclusions as to whether Z-scores are able to discriminate between "good" and "bad" structures for pKa calculations.

View larger version (52K): [in this window] [in a new window]   Figure 4. Correlation between WHAT_CHECK Z-score cutoffs for structures and the percentage of these structures that give correct results () or wrong results () when the local ID criteria are used to identify the proton donor. () The total number of structures that fulfill the cutoff criteria. (a) 1st generation packing quality, (b) 2nd generation packing quality, (c) Rotamer normality score, (d) Ramachandran plot appearance, (e) Backbone conformation score, (f) the sum of all Z-scores.

View larger version (103K): [in this window] [in a new window]   Figure 5. Glu 35 and Asp 52 from 2LZT (red), 4LYT (green), and 7LYZ (yellow). Only very small differences between the positions of these two side chains are observed for these three structures. The figure was prepared with the MOLSCRIPT (Kraulis 1998) and Raster3D (Merrit and Bacon 1997) programs.

View larger version (16K): [in this window] [in a new window]   Figure 6. The circular hydrogen-bond network around Asp 52 in 2LZT. Only side chains and their hydrogen bonds are shown. The notation used in the text for this hydrogen-bond network is Asp 52 – Asn 44 – Asn 46 – Ser 50 – Asn 59 – (Asp 52). Figure prepared with WHAT IF (Vriend 1990).

The solution structure of HEWL
The fact that 2LZT and 7LYZ give good results with the WHAT IF pKa calculation package for the two active-site residues does not necessarily mean that these two structures provide a more accurate description of the environment of Glu 35 and Asp 52 than does 4LYT. It is possible that 2LZT and 7LYZ give good results merely because many pKa calculation packages (including the one we use here) were calibrated on a set of X-ray structures that includes 2LZT. This would provide an example of "getting the right result for the wrong reason", and in order to exclude this possibility it is important to verify that the hydrogen-bond network around Asp 52 is formed in solution in the same way as it is in the two crystal structures 2LZT and 7LYZ. The PDB contains a single NMR structure of HEWL (PDBID: 1E8L; Schwalbe et al. 2001) which contains 50 models. Considering the resolution of NMR structures, we find Asn 44, Asn 46, and Asn 59 to be within hydrogen-bonding distance of Asp 52. Further evidence for the donation of hydrogen bonds to both Os of Asp 52 comes from an X-ray structure of HEWL D52S (Hadfield et al. 1994). The X-ray structure (PDBID: 1LSY) shows clear differences in the position of residues 44 and 46 when Asp 52 is mutated to serine, thus indicating that these residues are dependent on the Asp 52 side chain for proper orientation.

Calculating pKa values from the average titration curves of the 50 structures in 1E8L gives pKa values of 9.00 and 4.50 for Glu 35 and Asp 52, respectively. Although Glu 35 thus would be identified by the local ID criteria by using the NMR ensemble, it is clear that the pKa values of the two active-site residues are quite different from the experimentally measured values, and we therefore conclude that NMR structures should be used with caution for pKa calculations.

Calculating the pKa values of Glu 35 and Asp 52 correctly is thus highly dependent on having the correct hydrogen-bond network around Asp 52 in 2LZT, 4LYT, and 7LYZ. The generality of this conclusion is confirmed by examining the number of hydrogen bonds to Asp 52 for all of the structures in Table 1. In only six of the structures where both Asp 52 Os accept hydrogen bonds is Glu 35 not identified, and in four of these cases, a change in one of the pKa values of less than 0.10 units would have allowed identification of Glu 35 according to the local ID criteria. Glu 35 is identified in several structures where Asp 52 accepts only a single hydrogen bond. In these cases the environment provides several interactions that favor a negative charge on Asp 52, although the interactions cannot be described as hydrogen-bonds, because the angular criteria and distance criteria that normally define hydrogen bonds are not fulfilled.

Optimizing X-ray structures
We investigate the effect of using EMs, MD simulations, and CONCOORD (de Groot et al. 1997) analysis to optimize each of the 41 HEWL X-ray structures. EM, MD, and CONCOORD analyses were performed as described in Materials and Methods, and pKa values were computed for the final EM structure, for the average structure and for the minimized average structure for both the short MD simulations (200 psec) and the CONCOORD analysis. After the first 100 psec of each simulation, snapshot structures were recorded every 10 psec from the long MD simulations and used directly for pKa calculations. The final pKa values from the snapshot calculations were arrived at either by taking the average of the calculated pKa values for each frame or by averaging the titration curve for each group over all snapshots and determining the pKa value from this average titration curve (van Vlijmen et al. 1998).

Table 5 shows the results from the MD simulation protocols when the 41 X-ray structures were submitted to the above analyses. It is seen that none of the optimization methods are able to improve the frequency with which Glu 35 is identified as the proton donor. Instead, all methods used here make the predictions more incorrect, with the MD-based methods being worst of all.

Neither when using the average structure from an MD simulation nor when using the energy minimized average structure for MD simulations does the overall RMSD of the calculations improve. Similarly, a significant decrease is seen in the accuracy of the predictions for Glu 35 and Asp 52, with both residues predicted to have a pKa value below 5.

Setting the initial protonation states
The MD simulation methods used in Table 5 were performed with all titratable groups in their standard ionization state at pH 7.0. To examine whether the charges on Glu 35 and Asp 52 force HEWL to adopt a conformation in the simulations that favors negative charges on both these residues, we performed three additional MD simulations where we protonated either Glu 35 or Asp 52 or both residues. We note that HEWL contains only one histidine residue and that this His is more than 15 Å removed from the active-site residues, and the MD simulations with alternative protonation states can therefore be carried out without adjusting any other protonation states.

If the pKa values calculated from the snapshots of these MD simulations correspond to the protonation states used in the simulations, it means that the HEWL structure adjusts to stabilize the protonation state imposed by the set-up for the MD simulation. Table 6 shows that this is indeed the case. Protonating either Glu 35 or Asp 52 leads to a strong shift in the calculated pKa value for the protonated residue compared to the pKa value calculated from the trajectory where both groups are in their unprotonated form.

The protonation states used in the MD simulations thus determine the HEWL active-site pKa values that are calculated from the trajectory. More importantly, however, a comparison of the calculated pKa values for the HEWL X-ray structures crystallized at different pH values with the results of the simulations in Table 6 reveals that the two HEWL X-ray structures at pH 9.5 (1HSX and 1HSW) give significantly different results from the MD simulation with both Glu 35 and Asp 52 in their deprotonated form. Because the experimental pKa values of both Glu 35 and Asp 52 are significantly lower than 9.5 (Glu 35 pKa =6.20, Asp 52 pKa =3.68), this strongly suggests that the structural changes that occur in the Glu 35 -, Asp 52 - simulation are artefacts of the MD simulation, and not an accurate description of the actual structural changes that occur upon deprotonation of Glu 35 and Asp 52. Consequently, incorporating protonation and deprotonation reactions in MD simulations therefore cannot be expected to lead to structures that are more accurate for pKa calculations. Similarly, coupling pKa calculations and MD simulations is not likely to give any improvement in the realism of the MD simulation or in the accuracy of the calculated pKa values, as long as MD simulations are not capable of reproducing the structural changes (or lack of structural changes) that occur upon deprotonation of a residue.

	Discussion

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
The calculation of pKa values for active-site residues in enzymes is of importance in the study of enzyme mechanisms. In the present study, we investigated the best strategy for identifying the proton donor in the catalytic mechanism of an enzyme using pKa calculations on a set of X-ray structures. The results show that the best strategy is to use the unmodified X-ray structures directly for pKa calculations and identify the proton donor as the residue that fulfills the local ID criteria in the majority of the structures. We find that an improvement in the results can be achieved by excluding structures that make crystal contacts near the active-site residues. Excluding the HEWL structures from the P43 21 2 space group improves the prediction accuracy from 29/41 (70.7%) to 18/21 (85.7%), and eliminates the only two structures that identify Asp 52 as the proton donor.

The resolution of the X-ray structures, the crystallization conditions, and the quality of the X-ray structures did not show any significant correlation with the ability to correctly identify the proton donor, and none of these criteria should therefore be used to select a more trustworthy subset of X-ray structures. It should be noted, however, that the two high-resolution structures that we examined both gave the correct prediction, and we therefore cannot exclude the possibility that high-resolution structures will give more accurate results than low-resolution structures. We would indeed expect this to be true, but the present data do not warrant such a conclusion.

The use of energy minimizations, molecular dynamics simulations, and CONCOORD analysis marginally improves the overall agreement with experimental pKa values, but significantly degrades the ability to correctly identify the proton donor in HEWL. In our hands, MD simulations are particularly ill-suited for producing both average and snapshot structures for pKa calculations, because the calculated pKa values are highly dependent on the protonation state used in the MD simulation. Furthermore, it is evident that the structure produced by the MD simulations for a particular protonation state of Glu 35 and Asp 52 (Table 6) is different from the X-ray structure solved at a pH value where these two residues populate the same protonation state (Table 1), as seen by the large differences in the calculated pKa values. Here we have shown this to be the case only for a particular MD simulation package and for a particular MD simulation setup, but we expect similar results for other MD simulation packages, because the GROMACS standard force field and the simulation protocol used here do not differ significantly from other MD force fields and MD protocols. These conclusions are in contrast to the those reported by Van Vlijmen et al. (1998) and Gorfe et al. (2002), who argued that the use of MD improves the accuracy of pKa calculations, but our conclusions agree well with the observations of Koumanov et al. (2001) and Wlodek et al. (1997). We note that Van Vlijmen et al. (1998) and Gorfe et al. (2002) focused on the overall RMSD between experimental and calculated pKa values, and we believe that this is the main reason that they arrived at the conclusion opposite of ours. In our opinion, pKa calculations are interesting mostly for the very low number of active-site groups and buried groups that have highly shifted pKa values. The overall RMSD between calculated and experimentally measured pKa values is dominated by the contribution from surface-group pKa values that typically are almost unperturbed by the protein environment, and observing changes in the overall RMSD will therefore give a misleading picture of the effect of MD simulations and other protein structure improvement tools if one is primarily interested in the pKa values of functionally important groups.

The results presented here rely on the assumption that the WHAT IF pKa calculation method (Nielsen and Vriend 2001) will give the correct result when used with the "correct" structure. We define the correct structure as the structure of HEWL which is predominant in solution at the temperature and at the concentration that was used in the NMR experiments for determining the experimental pKa values. We have carefully examined the crystal contacts made in each crystal, and only for structures in the P 43 21 2 space group do we find crystal contacts in the vicinity of the two active acids. Because we thus examine only the pKa values of the two active-site acids Glu 35 and Asp 52, and because we get much worse results with the structure in the P 43 21 2, we are confident that the structures of the HEWL active site is a good representation of the solution structure of the active site, and consequently that the WHAT IF pKa calculation package calculates correct pKa values when used with the correct solution-like structure. This is furthermore corroborated by comparing the X-ray structures of HEWL to the NMR structure of the enzyme. Specifically it is important that the critical hydrogen-bond network around Asp 52 is formed in the NMR structure, as mentioned earlier. We have furthermore assumed that the WHAT IF pKa calculation routines are not biased by the choice of protein dielectric constant towards giving good results with only a certain subset of X-ray structures. We did this by recalculating all pKa values for the structures in Table 1 with a dielectric constant of both 4 and 16 and found that the correlations between the "correctness" of the results and the resolution, WHAT_CHECK scores, and crystallization conditions did not improve. The conclusions that we arrive at here are therefore not biased by our choice of dielectric constant. We were able to identify Glu 35 in only one of the 41 structures when using a dielectric constant of 16 but we observed a large drop in the overall RMSD, thus illustrating that a high dielectric constant will give poor pKa values for active-site residues, but accurate pKa values for surface residues. Using a dielectric constant of 4 gives a higher overall RMSD, and only a slightly worse performance for the active-site residues (Table 5).

In summary, we have proposed a set of criteria for identifying the proton donor from two candidate acidic residues in an enzyme. We have applied these criteria with the WHAT IF pKa calculation routines (Nielsen and Vriend 2001) to a set of 41 HEWL X-ray structures and showed that we successfully identify the correct residue in 85% of the structures if we exclude structures with crystal contacts near the active site. We furthermore find that EM, MD, and CONCOORD analyses are not currently capable of successfully optimizing protein X-ray structures for pKa calculations, and given the generality of the energies calculated in pKa calculation algorithms, we expect this conclusion to be true also for other structure-based energy calculation methods such as drug docking algorithms, protein design algorithms, and protein structure analysis tools.

In the light of this we consider it essential that the development and optimization of protein structure-based energy calculation methods are concerned not only with the construction of the algorithm for a few test cases, but also with the development of a specific protocol for preparing protein structures for the algorithm in question. Such protein preparation protocols should take into account the sensitivity of the algorithm in question, the source of the protein structure (NMR, X-ray, or homology model), and the desired level of detail of the results. Preferably the protein preparation protocol should be an integrated part of the energy-calculation procedure such that the X-ray structure is optimized simultaneously with the calculation of the desired energetic quantity. We note that the pKa calculation method developed by Alexov and Gunner (1999) presents an example of such a method. Unfortunately this method has not been benchmarked, and we are therefore unable to comment on the prediction accuracy of the method. Work on integrating structure optimization tools with structure-based energy calculation methods is on-going in our lab.

	Materials and methods

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
Selection of X-ray structures
The PDB was searched for structures which had a 100% sequence identity to the HEWL sequence in 2LZT. All structures that contained a significant number of ions or cofactors were excluded, and we arrived at a set of 64 wild-type HEWL structures. From this set we manually selected 41 structures that cover a wide range of crystallization conditions and resolutions.

pKa calculations
pKa calculations were carried out with the WHAT IF (Vriend 1990) pKa calculation routines as described (Nielsen and Vriend 2001), with the exception that the protein dielectric constant was set to 8 for all titratable groups. The WHAT IF pKa calculation routines perform a global optimization of the hydrogen-bond network for every single protonation state needed in FDPB-based pKa calculations. We employ the hydrogen-bond optimization algorithm developed by Hooft et al. (1996a) to produce the optimal hydrogen-bond network for every protonation state. The algorithm by Hooft et al. (1996a) does not change heavy atom positions, except in the cases where a better hydrogen-bond network can be produced by flipping the 2, 2, or 3 angles of His, Asn, and Gln, respectively. The WHAT IF pKa calculation routines employ DelPhi II (Nicholls and Honig 1991) for solving the Poisson-Boltzmann equation and use the OPLS forcefield as source of charges and radii (Jorgensen and Tirado-Rives 1988).

Retrieval of WHAT_CHECK scores
WHAT_CHECK (Hooft et al. 1996b) scores were retrieved from the compilation of WHAT_CHECK reports on all PDB files, which can be accessed at http://www.cmbi.kun.nl/pdbreport/. The Z-scores reported for "Users of a structure" were used.

Energy minimizations and molecular dynamics simulations
The GROMACS (Lindahl et al. 2001) molecular dynamics package (version 3.0) was used for all energy minimizations (EMs) and molecular dynamics (MD) simulations. The standard GROMACS force field was used in all calculations.

All EMs and MD simulations were performed in a box of water with a minimum distance between the edge of the box and the protein of 5 Å.

EMs were carried out with a steepest descents algorithm until the largest force was less than 2000 kJ mol^-1 nm^-1.

Two types of MD simulations were carried out: 200 psec (short) simulations and 1 nsec (long) simulations. Both simulations were preceded by steepest descents minimization as described above and by a 0.5 psec MD simulation where all protein atoms were kept fixed. The MD simulations were carried out with a Berendsen temperature coupling (Berendsen et al. 1984) to a bath at 300K. The MD stepsize was 2 fs, and a 10 Å cut-off was used for coulombic interactions which were calculated with a dielectric constant of one. All remaining MD parameters were set as described at http://www.gromacs.org/documentation/reference_3.0/online/getting_started.html#full.

Average structures
Average structures of the short MD simulations were calculated as the average position of all protein heavy atoms during the last 100 psec of the run using the program g_covar of the Gromacs 3.0 package.

Minimized average structures
Minimized average structures were produced by performing a steepest descents energy minimization of the average structures obtained from the short MD runs. The energy minimization was carried out until the largest force was smaller than 2000 kJ mol^-1 nm^-1.

Snapshot structures
The first 100 psec of the long MD simulations were discarded, and snapshot structures were taken every 10 psec. Each frame was used directly for pKa calculations. The final pKa value for each residue in the trajectory was calculated in two ways: (1) by calculating the average pKa value from all the individual snapshot pKa values, and (2) by calculating the average titration curve from all snapshots and thereafter calculating a pKa value from the average titration curve.

CONCOORD analysis
Coordinate sets were produced by the CONCOORD method (de Groot et al. 1997). Two thousand structures were generated. The average structure of the 2000 structures was calculated by the g_covar program of the GROMACS package.

The energy minimized structure of the average structure was calculated by performing a steepest descents energy minimization as described above.

	Acknowledgments

 
J.E.N. acknowledges financial support from the Danish Natural Science Council; Dr. Stewart Adcock and Dr. Wibke Sudholt for helpful discussions and for thoughtful comments on the manuscript; and Julie Schames for suggesting MD simulations with alternative protonation states. This work was supported in part by grants from the NIH, NSF, NBCR, and the W.M. Keck Foundation.

The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 USC section 1734 solely to indicate this fact.

	References

TOP Abstract Introduction pKa calculation algorithms Calibrating pKa calculation... Identifying the proton donor... Hen egg white lysozyme Objective Results Discussion Materials and methods References

 
Alexov, E.G. and Gunner, M.R. 1997. Incorporating protein conformational flexibility into the calculation of pH-dependent protein properties. Biophys. J. 72: 2075–2093.[Abstract/Free Full Text]

———. 1999. Calculated protein and proton motions coupled to electron transfer: Electron transfer from QA- to QB in bacterial photosynthetic reaction centers. Biochemistry 38: 8253–8270.[CrossRef][Medline]

Amadei, A., Linssen, A.B., and Berendsen, H.J. 1993. Essential dynamics of proteins. Proteins 17: 412–425.[CrossRef][Medline]

Antosiewicz, J., McCammon, J.A., and Gilson, M.K. 1994. Prediction of pH-dependent properties of proteins. J. Mol. Biol. 238: 415–436.[CrossRef][Medline]

———. 1996. The determinants of pKas in proteins. Biochemistry 35: 7819–7833.[CrossRef][Medline]

Baptista, A.M., Martel, P.J., and Petersen, S.B. 1997. Simulation of protein conformational freedom as a function of pH: Constant-pH molecular dynamics using implicit titration. Proteins 27: 523–544.[CrossRef][Medline]

Bashford, D. and Karplus, M. 1990. pKa’s of ionizable groups in proteins: Atomic detail from a continuum electrostatic model. Biochemistry 29: 10219–10225.[CrossRef][Medline]

Berendsen, H.J., Postma, J.P.M., DiNola, A., and Haak, J.R. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81: 3684–3690.[CrossRef]

Biswal, B.K., Sukumar, N., and Vijayan, M. 2000. Hydration, mobility and accessibility of lysozyme: Structures of a pH 6.5 orthorhombic form and its low-humidity variant and a comparative study involving 20 crystallographically independent molecules. Acta Crystallogr. D Biol. Crystallogr. 56: 1110–1119.[CrossRef][Medline]

Braxenthaler, M., Unger, R., Auerbach, D., Given, J.A., and Moult, J. 1997. Chaos in protein dynamics. Proteins 29: 417–425.[CrossRef][Medline]

Carugo, O. and Argos, P. 1997. Protein–protein crystal-packing contacts. Protein Sci. 6: 2261–2263.[Abstract]

Claussen, H., Buning, C., Rarey, M., and Lengauer, T. 2001. FlexE: Efficient molecular docking considering protein structure variations. J. Mol. Biol. 308: 377–395.[CrossRef][Medline]

Coutinho, P.M. and Henrissat, B. 1999. Carbohydrate-active enzymes: An integrated database approach. In Recent advances in carbohydrate bioengineering (eds. H.J. Gilbert, G. Davies,