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 Warwicker, J. Search for Related Content PubMed PubMed Citation Articles by Warwicker, J. Social Bookmarking                   What's this?

Improved pK_a calculations through flexibility based sampling of a water-dominated interaction scheme

Department of Biomolecular Sciences, University of Manchester Institute of Science and Technology (UMIST), Manchester M60 1QD, United Kingdom

Reprint requests to: Jim Warwicker, Department of Biomolecular Sciences, UMIST, P.O. Box 88, Manchester M60 1QD, UK; e-mail: jim. warwicker{at}umist.ac.uk; fax: +44-(0)161-236-0409.

(RECEIVED April 2, 2004; FINAL REVISION June 30, 2004; ACCEPTED July 6, 2004)

	Abstract

TOP Abstract Introduction Results and Discussion Materials and methods References

 
Ionizable groups play critical roles in biological processes. Computation of pK_as is complicated by model approximations and multiple conformations. Calculated and experimental pK_as are compared for relatively inflexible active-site side chains, to develop an empirical model for hydration entropy changes upon charge burial. The modification is found to be generally small, but large for cysteine, consistent with small molecule ionization data and with partial charge distributions in ionized and neutral forms. The hydration model predicts significant entropic contributions for ionizable residue burial, demonstrated for components in the pyruvate dehydrogenase complex. Conformational relaxation in a pH-titration is estimated with a mean-field assessment of maximal side chain solvent accessibility. All ionizable residues interact within a low protein dielectric finite difference (FD) scheme, and more flexible groups also access water-mediated Debye-Hückel (DH) interactions. The DH method tends to match overall pH-dependent stability, while FD can be more accurate for active-site groups. Tolerance for side chain rotamer packing is varied, defining access to DH interactions, and the best fit with experimental pK_as obtained. The new (FD/DH) method provides a fast computational framework for making the distinction between buried and solvent-accessible groups that has been qualitatively apparent from previous work, and pK_a calculations are significantly improved for a mixed set of ionizable residues. Its effectiveness is also demonstrated with computation of the pH-dependence of electrostatic energy, recovering favorable contributions to folded state stability and, in relation to structural genomics, with substantial improvement (reduction of false positives) in active-site identification by electrostatic strain.

Keywords: protein electrostatics; pK_as; ionization entropy; side-chain packing; active-site identification; structural genomics

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

	Introduction

TOP Abstract Introduction Results and Discussion Materials and methods References

 
Ionizable group interactions are important factors in various biological processes (Warshel 1981, 2003; Honig and Nicholls 1995; Warshel and Papazyan 1998; Simonson 2001), and are a focus of attempts to identify functional sites for structural genomics (Elcock 2001; Ondrechen et al. 2001). Grid-based continuum electrostatics methods, such as Finite Difference Poisson-Boltzmann (FDPB; Warwicker and Watson 1982; Klapper et al. 1986; Warwicker 1986) using a low protein-relative dielectric (typically _p = 4; Gilson and Honig 1986) have been applied to pK_a calculations (Bashford and Karplus 1990). These can be useful when applied to regions with limited solvent accessibility (SA; Demchuk and Wade 1996; Warwicker 1998), but _p = 4 calculations based on a single conformer have been largely unreliable for overall pK_a analysis. Generally _p = 20 performs better (Antosiewicz et al. 1994, 1996), presumably accounting to some degree for other factors (Schutz and Warshel 2001) such as conformational variation (Simonson and Perahia 1995), proton/hydrogen-bond network relaxation (Nielsen et al. 1999), or specific internal water binding (Fitch et al. 2002). Indeed, a Debye-Hückel (DH) model with water dielectric also gives reasonable agreement over-all for pK_as, and for the pH dependence of folding energy when combined with a simple model for ionizable group interactions in the unfolded state (Warwicker 1999). A Tanford-Kirkwood model gives reasonable estimates for carboxylate pK_as in ubiquitin (Sundd et al. 2002). As with FDPB/_p = 20, the DH and Tanford-Kirkwood models fail to generate the large pK_as often associated with functional groups (Warwicker 1998).

The problem of accounting for large and small pK_as in one scheme is being addressed with inclusion of conformational and proton configurational relaxation in _p = 4 calculations. Multiple conformers are generally sampled from molecular dynamics simulations (You and Bashford 1995; Zhou and Vijayakumar 1997; van Vlijmen et al. 1998; Koumanov et al. 2001; Gorfe et al. 2002), with some improvement in match to experiment and significant increase in computational requirement. A key issue is to address those conformational adjustments of most relevance to a pH titration. Multiconformation continuum electrostatics (MCCE) samples side-chain ionization and conformation (Alexov and Gunner 1997; Georgescu et al. 2002; Alexov 2003), yielding a pK_a root-mean-square (RMS) error of 0.83. An approximation in this model is that all conformers are combined to produce a single dielectric boundary. A method that assigns higher or lower electrostatic screening functions, (applied to Coulomb potentials), according to the hydrophobicity/hydrophilicity of residue microenvironments performs well, giving an RMS error of 0.5 for a large pK_a set (Mehler and Guarnieri 1999).

The current work also pursues two interaction schemes, the DH model that is effective for flexible groups, and FDPB/_p = 4, which can give large pK_as of functional interest, with combination in a framework of conformational relaxation. FD (for FDPB/_p = 4) interactions are always sampled for the experimental conformer, and DH interactions are selectively introduced as a mimic for relaxation from this conformer. Two examples illustrate this idea. First, a buried lysine, neutral at pH 7 due to dehydration (lower residue in Fig. 1). The DH model omits the Born term and underestimates such a pK_a, whereas FD with the Born term provides the required destabilization. Because DH would give a higher statistical weight than FD, a combined scheme must exclude DH and reflect the restricted conformation. Secondly (top of Fig. 1), in a relatively flexible salt bridge the Born contributions are balanced by favorable charge–charge interactions. For a pH shift that neutralizes one of the partners, without conformational relaxation FD will record a Born term for the remaining ionized group that is not balanced by charge–charge interactions (other than potential hydrogen bonds). The ionized partner is likely to seek a more solvent accessible conformation, reducing the Born energy penalty. Interactions in this water-dominated situation could be estimated by DH modeling of the salt bridge, as an alternative to FD calculations on a range of generated conformers.

View larger version (18K): [in this window] [in a new window]   Figure 1. Schematic diagram of ionizable group relaxation with pH in two different environments. An upper salt bridge can alter hydration upon pH titration, while the lower basic group cannot (e.g., a lysine with a reduced pKa).

An earlier empirical analysis of hydration entropy and pK_a calculations (Warwicker 1997) has been extended with study of relatively buried ionizable groups in the FD/_p = 4 model, and comparison made to small molecule ionization entropies. It is concluded that pK_a adjustments due to hydration entropy are relatively small, except in the case of cysteine. Changes in hydration entropy that can be important in binding energetics (Jung et al. 2002) are discussed in the framework of the empirical analysis.

	Results and Discussion

TOP Abstract Introduction Results and Discussion Materials and methods References

 
Estimation of hydration entropy changes upon ionization and charge burial
A subset of ionizable groups with large pK_as and surrounding ionizations that can be reasonably assigned at a pH around the pK_a of interest are studied, partitioning the analysis from the wider prediction of pH dependence (Table 1). Figure 2 shows that these groups are relatively buried in terms of subsequent FD/DH analysis (Fig. 3), and estimates of changes in hydration entropy are made with the single conformer FD/_p = 4 model. Polar hydrogen optimization is applied for the hydroxyl groups of serine and threonine. Tyrosine hydroxyls are included where the unionized form is used, which for entropic term modeling is in all calculations other than those of Figure 2C. A previous study (Warwicker 1997) looked at carboxylate and thiolate groups in a range of SA environments, deriving differential hydration numbers upon ionization (N_s) of about 2 (carboxylate) and 6 (thiolate). Current work extends the analysis for groups with large pK_as, and includes polar hydrogen optimization.

View larger version (80K): [in this window] [in a new window]   Figure 2. Active-site groups, with conformational restriction, used to derive Ns values for hydration entropy modification. The Born and Q–Q (interaction with background and other ionizable group charge), pKa contributions are calculated components, which together with the experimental pKa (Expt), and the fraction of the hydration shell buried in the protein relative to isolated amino acid (Vf), are used to derive Ns (see Materials and Methods). (A) Glutamic acid; (B) Cysteine; (C) Tyrosine; (D) Lysine.

View larger version (26K): [in this window] [in a new window]   Figure 3. Schematic diagram for determination of access to DH interactions. (A) Solvent probing of a side chain gives two SA arcs according to whether probed in the context of the drawn rotamer set (SAconfig) or fixed atoms only (SAfixed-atoms). (B) A different rotamer set gives higher solvent accessibility, such that [SAmax/SAfixed-atoms] 0.75 and access to DH interactions is allowed. (C) In this case, the rotamer set of B is not possible, SA remains small, and DH interactions are disallowed.

Hen egg white lysozyme is a well-known model for electrostatic calculations, particularly E35 (pK_a = 6.1; Kuramitsu and Hamaguchi 1980). Figure 2A shows the restricted environment of E35 and the elements that give N_s = 1.5, again with qualitative agreement to the experiment without N_s modification (Warwicker 1997). If discrepancies in the carboxylate group calculations of Figure 2A are approximated as hydration entropy within the FD framework, then a related single N_s value is relatively small and positive (Table 2; Warwicker 1997).

Calculations (Fig. 2B) with a reduced DsbA crystal structure (Guddat et al. 1998) give N_s = 4.2, where H60, E24, and E38 were assigned neutral (Warwicker 1998), and the measured C30 pKa is 3.4 (Grauschopf et al. 1995). A reduced model of the D26N mutant was made from oxidized, wild-type thioredoxin, uncoupling the similar pK_as of C32 and D26, with C32 pK_a = 7.5 (Chivers et al. 1997) and N_s = 4.9. Calculated N_s of 4.5, 4.2, and 4.9 imply significant discrepancies in FD/_p = 4 calculations of pK_as for relatively buried cysteines, without the N_s term.

For tyrosine, the active-site residue Y9 of glutathione S-transferase A1–1 has a low pK_a (8.1) in the substrate-free enzyme (Björnestedt et al. 1995). Calculations with a monomer from the dimeric crystal structure (Cameron et al. 1995) give N_s = 1.4 (Fig. 2C). A pK_a of 6.1 has been measured for Y149 in the active site of UDP-galactose-4-epimerase with NAD⁺ bound (Liu et al. 1997). Removing UDP, but not NAD⁺, from the crystal structure (Thoden et al. 1996) gives N_s = 0.8 (Fig. 2C), with a strong interaction to the positively charged side chain of K153 and a hydrogen bond between the deprotonated Y149 side chain and a ribose hydroxyl group.

Hydration entropy and pK_as: Lysine, histidine, and arginine
Figure 2D shows the environment of K229 in rabbit muscle aldolase (Blom and Sygusch 1997), implicated in Schiff base formation, with pK_a probably matching the decline of activity around pH 6.5 (Morris and Tolan 1994). A monomer from the tetramer (with product removed) gives N_s = 0.2. The network of ionizable groups around K229 (Fig. 2D) were assigned unperturbed charges for neutral pH calculations. A catalytic antibody with aldolase activity has an active-site lysine (H chain K93) with a pK_a of 5.5 estimated from the pH dependence of enamine formation (Barbas et al. 1997). Figure 2D shows the enclosed environment of K93(H), with W103(H) adjacent and H27(L) distant, but displayed to facilitate the view. NZ atom locations are shown from 1axt [PDB] and modeled with more solvent exposure, because electron density for K93(H) does not extend beyond CE (Barbas et al. 1997). These give N_s = 2.0 (1axt [PDB] , quoted in Fig. 2D) and N_s = 0.8 (model), with the low pK_a dominated by the Born term.

Histidine 48 of porcine pancreatic phospholipase A₂ is catalytic (Thunnissen et al. 1990), with a pK_a of 5.5 in the calcium-bound form (Verheij et al. 1980). Calculation gives N_s = 2.0 (Table 2). In chymotrypsin, interactions between partially charged/protonated or partially charged/neutral forms of H57, and other catalytic residues D102 and S195 have been evaluated, with S195 hydroxyl directed toward neutral H57 and away from protonated H57. Comparison with the measured pK_a of 6.8 (Fersht and Renard 1974) gives N_s = 0.9 (Table 2). For T4 lysozyme, the noncatalytic residue H31 (pK_a = 9.1) is significantly buried in a salt bridge with D70 (Anderson et al. 1990), giving N_s = 3.1 (Table 2). This N_s-based study of buried side chains with known pK_as has excluded histidine residues with complicating factors such as strongly coupled titrations and ligand involvement. Variability in N_s values for histidine (prior to averaging) probably reflects relatively mixed environments and the omission of local relaxations (Edgcomb and Murphy 2002).

There is some evidence for reduced arginine pK_as (Lehoux and Mitra 1999; Morillas et al. 1999), but not a clear example coupled to an atomic structure. Arginine side chain and N-terminal groups are assigned (N_s) as the lysine side-chain average, and C-terminal groups follow the average of the carboxylate containing side chains.

Comparison with ionization data for amino acids and analogs
The empirical modification for hydration entropy is generally small in pK_a terms compared with the Born and charge–charge contributions, with the exception of cysteine (Fig. 2; Table 2). The sense of the modification is in line with water molecule release upon ionizable group burial being greater (albeit slightly in most cases) for ionized versus unionized forms. Table 3 shows the enthalpies and entropies of ionization for amino acid side chains and molecular analogs (Izatt and Christensen 1976). Hydration entropies make the major contributions to ionization entropies (Alberty 1983), which show a wide variation for acidic and basic side chains in amino acids, partially due to overlapping proton binding equilibria. In contrast, the analog data show similar ionization S within acidic and within basic groups, distinguishing between dissociation from one neutral to two charged species (acidic, larger S), and from one charged to one charged and one neutral species (basic, smaller S). Differencing between average S values for each grouping (–103.5 J/deg/mole for AH A^– B + H⁺ and –17.0 J/deg/mole for BH⁺ B + H⁺), cancelling the proton term, gives S =–86.5 J/deg/mole for AH A^– summed with B BH⁺. Approximating equal S for each of these unionized to ionized form transitions, S = –43.3 J/deg/ mole is obtained, equivalent to N_s = 1.7 (300 K) in the hydration shell model and consistent overall with the pK_a calculations (other than cysteine). This favorable comparison supports use of the empirical model.

View this table: [in this window] [in a new window]   Table 4. Comparison of ionization S values for carboxylates and thiolates

View larger version (114K): [in this window] [in a new window]   Figure 4. PSBD (blue backbone):E3 dimer (green and orange surfaces) interface, showing mutated charge network.

Combination of FDPB and DH models: The FD/DH interaction scheme
The empirical modification of the FD/_p = 4 single conformer pK_a calculations for hydration entropy is applied to the FD component of the FD/DH scheme. Figure 3 schematizes side-chain flexibility in different environments, and the derivation of SA_max and access to DH interactions. A VdW clash tolerance (VdW_tol) is included in mean-field calculations that determine allowed rotamers for SA analysis. This parameter is varied to approximate conformational relaxation not explicitly included (nonlibrary rotamer packing and main-chain movement), and the best value determined with respect to experimental pK_as (Fig. 5B). VdW_tol of around 0.8 Å is typically required to give a packing solution for united atom VdW radii, so that VdW_tol > 0.8 Å represents additional flexibility that will lead to greater SA_max and successive entry of more ionizable groups to the DH calculation regime.

View larger version (27K): [in this window] [in a new window]   Figure 5. VdWtol variation, SAmax, and RMS pKa errors. (A) Representative averages over a protein and active-site groups are shown for SAmax vs. VdWtol. Intersection of the SAmax = 0.75 and VdWtol = 1.4 lines gives FD/DH parameterization. (B) RMS pKa errors are shown for FD/DH calculations for each of the 110 and 117 group sets (with FD only and DH only at each extreme), and horizantal lines ("null") give RMS errors for model compound pKas.

Figure 5B shows that for the 110 groups, excluding most of the large active-site pK_as, FD performs poorly and DH well in comparison to the null hypothesis. For the 117 set, the large pK_as impair DH performance, but have relatively little impact overall on FD, which remains dominated by inaccuracy within the 110 set. A broad region of optimal fit to experiment is found as VdW_tol is varied for FD/DH. Increasing VdW_tol adds flexibility that reduces the erroneously large pK_as for surface groups in the 110 set (schematically at the top of Fig. 1), while too much flexibility (beyond VdW_tol = 1.4 Å) allows buried group relaxation (bottom of Fig. 1) to an extent that is inconsistent with experiment. From the optimal region, VdW_tol = 1.4 Å is taken for FD/DH pK_a calculations.

Comparison of pK_a calculation methods
Table 6 gives overall RMS pK_a errors for various methods. FD/_p = 20 (Antosiewicz et al. 1994) is similar to DH, consistent with a water-like environment (Warwicker 1999). Similarly, FD/_p = 4 with SA derived from VdW radii rather than a probe reentrant surface (Dong and Zhou 2002), has a large degree of internal water, and is comparable overall to DH. FD/DH performs moderately better with hydration entropy modification than without, with the largest difference for the active-site cysteine pK_as of papain and DsbA (Table 7). Overall results are improved somewhat with the simple polar hydrogen optimization algorithm, consistent with previous observations (Nielsen et al. 1999). Differences associated with application of either the hydration entropy or polar hydrogen placement schemes are generally smaller than those from FD/DH introduction relative to FD (Table 6). The MCCE method gives an RMS pK_a error for the 110 groups of 0.92 (Georgescu et al. 2002), compared with 0.58 for DH and 0.79 for FD/DH, while the single conformer FDPB method used by these authors gives 2.05, close to the current FD value of 2.10.

View larger version (26K): [in this window] [in a new window]   Figure 6. Scatter plots (calculation vs. experiment) for pKas and pKas with FD, DH, and FD/DH calculations (117 group set).

For 3rn3 [PDB] , groups with pK_a to experiment >1.5 are H48 and H119. There are alternate locations for H119 ("A" was used), and a sulphate ion that may give overestimation of stabilizing interactions (FD/DH pK_a = 8.3 and measured pK_a = 6.1). H48 of 3rn3 [PDB] has an FD/DH pK_a of 4.6 and a measured pK_a of 6.3. H48 is relatively buried and excluded from the DH scheme. FD interactions are dominated by an ion pair with D14 that varies between ribonuclease A structures (Georgescu et al. 2002).

Calculations for 2rn2 [PDB] were made without magnesium to match experiment. FD/DH discrepancies with experiment >1.5 pK_a units are E48 and H114. For E48 (1.7 calculated and 4.4 measured), DH is incorporated at VdW_tol = 1.4 Å but not at VdW_tol = 1.2 Å. However, the calculated pK_a remains the same, because interactions are dominated by hydrogen bonding to background charges in the FD scheme. It is possible that unfavorable interactions with the D10, D70, and D134 charge cluster surrounding E48 have been underestimated. These groups are DH accessible, so that relatively weak DH interactions dominate. For H114 an FD/DH pK_a at VdW_tol = 1.4 Å of 1.0 compares with a measured pK_a of 5.0. This residue is relatively buried, with an enforced FD-only scheme. The key (FD) interaction is a hydrogen bond to the peptide NH of C63, so that H114 protonation must involve some change in structure.

Of the seven additional proteins with individual groups, 1xnb [PDB] /E172 and 9pap [PDB] /C25 give the worst results for FD/DH (Table 7). For 1xnb [PDB] , this discrepancy (calculated pK_a = 4.4, measured = 6.7) is a failure of the FD/DH scheme, because FD alone gives a large positive pK_a for E172 with or without hydration entropy modification. Again, more detailed FD analysis of coupled clusters may improve the FD/DH model, a conclusion supported for xylanase by pK_a calculations using subsets of ionizable groups (Nielsen and McCammon 2003).

The FD/DH pK_a for C25 in papain (9pap [PDB] ) is 5.9, compared with 3.3 by experiment, and in contrast the FD pK_a is a good match. Both C25 and H159 are accessible to the DH scheme at VdW_tol = 1.4 Å. Because the FD result corresponds to significant stabilization of C25, it might be expected to figure more strongly in FD/DH calculation at VdW_tol = 1.4 Å. Of the four possible samplings of C25 and H159 net charged forms, only one (both FD) will give a large favorable interaction. Alteration from FD to FD/DH may therefore result from relative weighting toward weak (DH) interactions. Active-site residue Y9 of 1gsd [PDB] is allowed access to the DH scheme at VdW_tol = 1.4 Å, (FD/DH pK_a = 9.5 and measured 8.1), but not at VdW_tol = 1.3 Å (calculated pK_a = 8.7). It is tightly coupled to R20, which also has a DH accessibility transition over this VdW_tol range.

Ionizable group electrostatic energy contributions
A move away from neutral pH generally destabilizes the folded state (Fink et al. 1994), consistent with a stabilizing network of positive and negative charges (Wada and Nakamura 1981), and with destabilization by charge reversal of surface amino groups (Hollecker and Creighton 1982). Computational methods with FD/_p = 20 or DH interactions accommodate such pH dependence, albeit with adjustment for pH-dependent effects in unfolded states (Schaefer et al. 1997; Warwicker 1999).

However, use of single conformer FD/_p = 4 has given rise to questions of salt-bridge stability (Hendsch and Tidor 1994; Dong and Zhou 2002), as is shown with the mostly unfavorable contributions of ionizable group interactions to stability at pH 7 (FD column of Table 8). In contrast to the DH column, FD calculations suggest that most of these folded proteins would be more stable with global replacement of ionizable groups, at variance with the prior discussion. The FD/DH method recovers favorable contributions to neutral pH stability. In many cases FD/DH matches DH closely, with DH interactions dominating FD overall. The particularly favorable values in Table 8 for 3icb [PDB] /calbindin with FD and FD/DH calculation result from the strong binding of calcium ions by acidic side chains.

View larger version (59K): [in this window] [in a new window]   Figure 7. Active-site identification using FD/DH. Four proteins (green backbones) used in pKa calculations are shown, with ionizable groups assessed as disallowed from DH interactions in orange, and a subset with calculated pKas between 3 and 10 shown in purple and labeled.

The effectiveness of FD/DH calculations in combining FD/_p = 4 and DH computation is illustrated in Figure 5B. The simple DH model is relatively good for the mostly small pK_as of the 110 groups set (Table 6), with an RMS pK_a error of 0.58, matching the performance of more complex methods (Mehler and Guarnieri 1999; Georgescu et al. 2002). However, DH performance is eroded when larger pK_as are included (117 groups). The FD/_p = 4 single-conformer method fails to consistently predict pK_as in the 110 set. As the system is relaxed in the FD/DH model, with selective access to DH interactions mimicking side-chain and possibly main-chain readjustment on a limited scale, an optimal relaxation is evident before key groups are erroneously solvent-exposed. Large-scale conformational change cannot be reliably predicted, and is not the subject of this work. A VdW_tol parameter controls conformational relaxation through mean-field rotamer packing and estimation of SA_max for each ionizable group (Figs. 3, 5A). At VdW_tol = 1.4 Å, FD/DH gives an RMS pK_a error of 0.86 for the 117 groups, compared with null, DH and FD values of 1.24, 1.18, and 2.06, respectively. The tendencies of FD to overestimate small pK_as and DH to underestimate large pK_as are clear in Figure 6. Factors that could further improve the FD/DH scheme include modeling of larger scale flexibility and pH-dependent ion binding, and detailed water structure (Koumanov et al. 2002). Tightly coupled clusters, with groups that are borderline for access to the DH scheme at VdW_tol = 1.4 Å, may benefit from FD analysis of individual rotamer combinations. This refers to relatively few groups; the majority would be included in FD/DH sampling.

Generally, the FD/DH model is improved by hydration entropy modification and polar hydrogen optimization, but not to the extent of FD/DH improvement over FD (Tables 6,7). Two areas further demonstrate potential. First (Table 8), FD/DH recovers overall stabilizing contributions for ionizable group interactions at pH 7, in contrast to the de-stabilization of the FD/_p = 4 single-conformer model. Second, Figure 7 demonstrates active-site identification, which could be used in a structural genomics context or provide subsets of residues to seed active-site–centered electrostatics calculations (Nielsen and McCammon 2003). The computational speed of the FD/DH method, coupled to accuracy over large and small pK_a deviations, make it suitable for large-scale database analysis.

	Materials and methods

TOP Abstract Introduction Results and Discussion Materials and methods References

 
Coordinates and pK_as
Proteins and coordinate sets (Berman et al. 2000) are given in Table 1. A previous study of computational methods (Georgescu et al. 2002) has collated pK_a data for the first eight proteins in Table 1, covering 126 ionizable groups, although the authors report that some of these ionizations may overlap either a limit of the measured pH range and/or a protein unfolding transition. Further investigation of the literature cited by Georgescu et al. (2002), with regard to these criteria, lead to the exclusion of the following groups in the current study: D54, E73, D93, D101 of 1a2p [PDB] ; Y23, Y35, and K41 of 4pti [PDB] ; D66 of 1b0d [PDB] ; D7, D27, Y31, and C-t of 1ppf [PDB] ; K13 of 1pga [PDB] ; D14 of 3rn3 [PDB] ; D102 and D148 of 2rn2 [PDB] . These 16 exclusions reduce 126 groups for the eight proteins studied by Georgescu et al. (2002) to 110, which also excludes the C-t of 1pga [PDB] .

This set of 110 groups was supplemented by a further seven active-site groups, with large pK_as for FD/DH analysis, forming the "117" set. Groups and coordinates used for development of the empirical hydration model (Warwicker 1997), including these seven active-site residues, are listed in Table 1. In addition, coordinate set 1ebd [PDB] was used to study interactions within the pyruvate dehydrogenase complex.

Electrostatics calculations (FD and DH)
DH (Warwicker 1999) and FD (_p = 4, single conformer) (Warwicker 1998) calculations followed previous work, except for hydration entropy and polar hydrogen optimization modifications (next section). Model compound pK_as used and charge assignment for ionizable group atoms were Arg 12.0 0.5/NH1 0.5/NH2; Lys 10.4 1.0/NZ; His 6.3 0.5/ND1 0.5/NE2; Asp 4.0–0.5/OD1–0.5/ OD2; Glu 4.4–0.5/OE1–0.5/OE2; Cys 8.3–1.0/SG (where not disulfide bonded); Tyr 10.2–1.0/OH; N-t 7.5 1.0/N; C-t 3.8–0.5/O–0.5/OXT. Cysteine side chains were excluded from pK_a calculations, except for papain C25, DsbA C30 and C33, and thioredoxin C32 and C35. Partial charges were assigned from the GROMOS library (van Gunsteren and Berendsen 1987). Where specified, _p = 20 was used in place of _p = 4 for the FD model. A water relative dielectric of 78.4 was used throughout. An alternative to a solvent-probed (probe radius =1.4 Å) reentrant surface was tested, a VdW sphere derived surface giving greater water penetration into the protein. All DH and FD computations were made with a linear ion response at 0.15 Molar.

The energy of a protonation microstate relative to the fully unprotonated state for M titratable groups is (modified from Bashford and Karplus 1990):

where x_m,x_n are members of vectors (lengths M) taking the values 0 (unprotonated) and 1 (protonated), q_m^o,q_n^o are the charges of the unprotonated sites m,n, and pK_m,model is the model compound pK_a for site m. Interactions of group m with the protein dielectric and nonionizable (background) charge environments are differenced to calculations for each ionizable group extracted from the protein, but remaining on the same FD grid (G_Born, G_back), while G_mn gives the interaction between ionizable groups. The term G_(TS)QWAT arises from the hydration entropy model discussed in the next section. This definition of G_s applies to either FD or DH calculations, but in the DH case G_Born and G_(TS)QWAT terms are zero. In addition G_(TS)QWAT is set to zero for _p = 20 calculations. Average fractional protonation for site m is derived as:

where a full evaluation of the partition function (Z) is possible for small numbers of groups using the reduced sites method (Bashford and Karplus 1990) at pH extremes, or alternatively monte carlo sampling of lowest energy states is used (Beroza et al. 1991). The calculated pK_a of site m is that pH for which x_m=0.5.

Electrostatic energy contributions from ionizable groups were calculated by summing increments over pH according to ( G)/ pH =2.303RTQ, where G, Q differences are relative to pH titration of the same set of ionizable groups with model compound pK_as (Antosiewicz et al. 1994). These sums were extended to an extreme pH at which a full evaluation of ionization states was possible, G =-RTlnZ, again with differencing to a calculation for the same groups titrating with model compound pK_as.

Hydration entropy modification and polar hydrogen placement
An empirical estimate of the contribution that water ordering makes to the change in ionization energy upon transfer into protein is written as G_(TS)QWAT = V_f•E_s (Warwicker 1997), where V_f is the fractional change in first hydration shell volume, and E_s is a free energy associated with water ordering for the complete shell. V_f is calculated from FD grids for a group in the protein relative to that whole amino acid extracted from the protein. It is convenient to consider E_s in terms of a notional number of water molecules (Warwicker 1997) in the first hydration shell (N_s), using 25 J/K/ mole entropy cost of immobilization of a single water molecule or 7.5 kJ/mole at 300 K. In pK_a calculations, E_s and N_s correspond to differences between water structure in ionized and unionized forms. Averages of N_s are taken over groups within a class, and individual N_s values derived with charge environments fixed according to expected ionization at physiological pH (see Results and Discussion). Thus, N_s derivation was separated from global pK_a computation in the FD/DH method through the use of separate calculations with fixed ionization on all groups but one. The N_s parameter gives the empirically derived difference in hydration numbers between ionized and unionized forms, for a complete hydration shell. It is uniform for a class of ionizable group, but each individual modification is proportional to burial (V_f).

Optimization of polar hydrogen locations can improve pK_a calculations (Nielsen et al. 1999; Koumanov et al. 2003). A simple procedure was implemented which optimizes polar hydrogen positions for the OH groups of Ser and Thr side chains (and Tyr if present in its neutral form). Torsions are sampled, with 12° resolution, against coulombic interactions from charges other than those belonging to the set of OH groups being optimized. Ionizable groups are specified in their ionized forms, to focus pK_a calculations on the net difference upon titration.

Mean-field packing algorithm and access to the DH interaction scheme
A mean-field algorithm (Koehl and Delarue 1994), adjusted with hard sphere clashes replacing Lennard-Jones interactions (Cole and Warwicker 2002), establishes a set of allowed rotamers from a standard library (Tuffery et al. 1997) as well as those in the experimental structure. The conformational matrix (CM) or weight for rotamer k of side chain i (of N side chains), depends on clashes with fixed atoms and the rotamers of other side chains j:

Weights for rotamers l =1 to Kj of side chain j are summed, and the multiplication of these side-chain and fixed-atom weights requires that each be nonzero for a packing solution for rotamer k of residue i. CM values are iterated from an initial assignment of equal weights. These mean-field methods were used previously to assess the entropy associated with side-chain rotamer distributions, but are used here simply to estimate possible rotamers for a given structure and packing tolerance (VdW_tol).

Hydrogen atoms are absent in the packing calculations, which use united atom VdW radii, so that a VdW_tol of about 0.8 Å is typically required to repack rotamers of the experimental structure, due particularly to side-chain/main-chain clashes. Above this value clashes may partially account for conformational relaxation, such as variation from library side-chain rotamers or small main-chain movements.

An estimate of maximal SA (SA_max) is obtained for each ionizable group. For each neighboring side chain, the (allowed) rotamer that provides the most SA for the charged atom of the ionizable group under study is used to mark SA information on a spherical polar grid around that atom. This grid is scanned after complete neighbor analysis for the overall SA_max. Maximal SA estimates for an ionizable group with more than one charged atom in our model are averaged. The likelihood that small conformational adjustment could lead to substantial solvent exposure of each group is assessed with SA_max (Fig. 3). If SA_max/SA_fixed-atoms 0.75, where SA_fixed-atoms is calculated in the context of fixed atoms only (i.e., excluding atoms that move in the rotamer library), then that group is allowed to sample DH (water-dominated) interactions.

Combined FD/DH calculations
In the FD/DH combination, each site m in the equation for microstate energy G_s is sampled as protonated or unprotonated in each of FD and DH (if allowed) schemes, with the component energies of G_s assigned appropriately. Thus, a standard single conformer pK_a calculation with 2^M states for M groups, approaches 4^M states in FD/DH, because most sites are relatively flexible and have access to the DH scheme. Note that for any group m sampled as DH, all mn (and nm) interactions are DH whether group n is sampled as FD or DH, that is, the water-dominated scheme persists for any interaction involving a group sampled in a presumed water-rich environment. Ionizable group array energies and pK_as in FD/DH calculations are derived as for FD or DH alone. No attempt is made to assign a density of conformational states, other than DH access or not, so that electrostatic energy determines the relative weightings of FD and DH. A moderate energy (favorable or unfavorable) from DH interaction would outweigh a highly unfavorable FD interaction, while a large and favorable FD term will dominate DH.

	Acknowledgments

 
The European Union is thanked for funding during the course of this work.

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 Results and Discussion Materials and methods References

 
Alberty, R.A. 1983. Physical chemistry, 6th ed. John Wiley & Sons, New York.

Alexov, E. 2003. Role of the protein side-chain fluctuations on the strength of pair-wise electrostatic interactions: Comparing experimental with computed pK_as. Proteins 50: 94–103.[CrossRef][Medline]

Alexov, E. and Gunner, M.R. 1997. Incorporating protein conformational flexibility into the calculation of pH-dependent protein properties. Biophys. J. 74: 2075–2093.

Anderson, D.E., Becktel, W.J., and Dahlquist, F.W. 1990. pH-induced denaturation of proteins: A single salt bridge contributes 3–5 kcal/mol to the free energy of folding of T4 lysozyme. Biochemistry 29: 2403–2408.[CrossRef][Medline]

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

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

Barbas III, C.F., Heine, A., Zhong, G., Hoffmann, T., Gramatikova, S., Björnestedt, R., List, B., Anderson, J., Stura, E.A., Wilson, I.A., et al. 1997. Immune versus natural selection: Antibody aldolases with enzymic rates but broader scope. Science 278: 2085–2092.[Abstract/Free Full Text]

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

Berman, H.M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T.N., Weissig, H., Shindyalov, I.N., and Bourne, P.E. 2000. The Protein Data Bank. Nucleic Acids Res. 28: 235–242.[Abstract/Free Full Text]

Beroza, P., Fredkin, D.R., Okamura, M.Y., and Feher, G. 1991. Protonation of interacting residues in a protein by a Monte Carlo method: Application to lysozyme and the photosynthetic reaction center of Rhodobacter sphaeroides. Proc. Natl. Acad. Sci. 88: 5804–5808.[Abstract/Free Full Text]

Björnestedt, R., Stenberg, G., Widersten, M., Board, P.G., Sinning, I., Jones, T.A., and Mannervik, B. 1995. Functional significance of arginine 15 in the active site of human class glutathione transferase A1–1. J. Mol. Biol. 247: 765–773.[CrossRef][Medline]

Blom, N. and Sygusch, J. 1997. Product binding and role of the C-terminal region in class I D-fructose 1,6-bisphosphate aldolase. Nat. Struct. Biol. 4: 36–39.[CrossRef][Medline]

Cameron, A.D., Sinning, I., L’Hermite, G., Olin, B., Board, P.G., Mannervik, B., and Jones, T.A. 1995. Structural analysis of human -class glutathione transferase A1–1 in the apo-form and in complexes with ethacrynic acid and its glutathione conjugate. Structure 3: 717–727.[Medline]

Chivers, P.T., Prehoda, K.E., Volkman, B.F., Kim, B.M., Markley, J.L., and Raines, R.T. 1997. Microscopic pK_a values of Escherichia coli thioredoxin. Biochemistry 36: 14985–14991.[CrossRef][Medline]

Cole, C. and Warwicker, J. 2002. Side-chain conformational entropy at protein–protein interfaces. Protein Sci. 11: 2860–2870.[Abstract/Free Full Text]

Demchuk, E. and Wade, R.C. 1996. Improving the continuum dielectric approach to calculating pK_as of ionizable groups in proteins. J. Phys. Chem. 100: 17373–17387.[CrossRef]

Dong, F. and Zhou, H.X. 2002. Electrostatic contributions to T4 lysozyme stability: Solvent-exposed charges versus semi-buried salt bridges. Biophys. J. 83: 1341–1347.[Abstract/Free Full Text]

Edgcomb, S.P. and Murphy, K.P. 2002. Variability in the pKa of histidine side-chains correlates with burial within proteins. Proteins 49: 1–6.[CrossRef][Medline]

Elcock, A.H. 2001. Prediction of functionally important residues based solely on the computed energetics of protein structure. J. Mol. Biol. 312: 885–896.[CrossRef][Medline]

Fersht, A.R. and Renard, M. 1974. pH-dependence of chymotrypsin catalysis. Appendix: Substrate binding to dimeric chymotrypsin studied by x-ray diffraction and equilibrium method. Biochemistry 13: 1416–1426.[CrossRef][Medline]

Fink, A.L., Calciano, L.J., Goto, Y., Kurotsu, T., and Palleros, D.R. 1994. Classification of acid denaturation of proteins: Intermediates and unfolded states. Biochemistry 33: 12505–12511.

Fitch, C.A., Karp, D.A., Lee, K.K., Stites, W.E., Lattman, E.E., and Garcia- Moreno, E.B. 2002. Experimental pK_a values of buried residues: Analysis with continuum methods and role of water penetration. Biophys. J. 82: 3289–3304.[Abstract/Free Full Text]

Georgescu, R.E., Alexov, E., and Gunner, M.R. 2002. Combining conformational flexibility and continuum electrostatics for calculating pK_as in proteins. Biophys. J. 83: 1731–1748.[Abstract/Free Full Text]

Gilson, M.K.