A self-consistent, microenvironment modulated screened coulomb potential approximation to calculate pH-dependent electrostatic effects in proteins.

An improved approach is presented for calculating pH-dependent electrostatic effects in proteins using sigmoidally screened Coulomb potentials (SCP). It is hypothesized that a key determinant of seemingly aberrant behavior in pKa shifts is due to the properties of the unique microenvironment around each residue. To help demonstrate this proposal, an approach is developed to characterize the microenvironments using the local hydrophobicity/hydrophilicity around each residue of the protein. The quantitative characterization of the microenvironments shows that the protein is a complex mosaic of differing dielectric regions that provides a physical basis for modifying the dielectric screening functions: in more hydrophobic microenvironments the screening decreases whereas the converse applies to more hydrophilic regions. The approach was applied to seven proteins providing more than 100 measured pKa values and yielded a root mean square deviation of 0.5 between calculated and experimental values. The incorporation of the local hydrophobicity characteristics into the algorithm allowed the resolution of some of the more intractable problems in the calculation of pKa. Thus, the divergent shifts of the pKa of Glu-35 and Asp-66 in hen egg white lysozyme, which are both about 90% buried, was correctly predicted. Mechanistically, the divergence occurs because Glu-35 is in a hydrophobic microenvironment, while Asp-66 is in a hydrophilic microenvironment. Furthermore, because the calculation of the microenvironmental effects takes very little CPU time, the computational speed of the SCP formulation is conserved. Finally, results from different crystal structures of a given protein were compared, and it is shown that the reliability of the calculated pKa values is sufficient to allow identification of conformations that may be more relevant for the solution structure.     

Rapid Calculation of Protein pK(a) Values Using Rosetta

We developed a Rosetta-based Monte Carlo method to calculate the pK(a) values of protein residues that commonly exhibit variable protonation states (Asp, Glu, Lys, His, and Tyr). We tested the technique by calculating pK(a) values for 264 residues from 34 proteins. The standard Rosetta score function, which is independent of any environmental conditions, failed to capture pK(a) shifts. After incorporating a Coulomb electrostatic potential and optimizing the solvation reference energies for pK(a) calculations, we employed a method that allowed side-chain flexibility and achieved a root mean-square deviation (RMSD) of 0.83 from experimental values (0.68 after discounting 11 predictions with an error over 2 pH units). Additional degrees of side-chain conformational freedom for the proximal residues facilitated the capture of charge-charge interactions in a few cases, resulting in an overall RMSD of 0.85 pH units. The addition of backbone flexibility increased the overall RMSD to 0.93 pH units but improved relative pK(a) predictions for proximal catalytic residues. The method also captures large pK(a) shifts of lysine and some glutamate point mutations in staphylococcal nuclease. Thus, a simple and fast method based on the Rosetta score function and limited conformational sampling produces pK(a) values that will be useful when rapid estimation is essential, such as in docking, design, and folding.     

A statistical approach to the prediction of pK(a) values in proteins

We propose a simple model for the calculation of pK(a) values of ionizable residues in proteins. It is based on the premise that the pK(a) shift of ionizable residues is linearly correlated to the interaction between a particular residue and the local environment created by the surrounding residues. Despite its simplicity, the model displays good prediction performance. Under the sixfold cross test prediction over a data set of 405 experimental pK(a) values in 73 protein chains with known structures, the root-mean-square deviation (RMSD) between the experimental and calculated pK(a) was found to be 0.77. The accuracy of this model increases with increasing size of the data set: the RMSD is 0.609 for glutamate (the largest data set with 141 sites) and approximately 1 pH unit for lysine, with a data set containing 45 sites.     

Electrostatic environment of hemes in proteins: pK(a)s of hydroxyl ligands

The pK(a)s of ferric aquo-heme and aquo-heme electrochemical midpoints (E(m)s) at pH 7 in sperm whale myoglobin, Aplysia myoblogin, hemoglobin I, heme oxygenase 1, horseradish peroxidase and cytochrome c oxidase were calculated with Multi-Conformation Continuum Electrostatics (MCCE). The pK(a)s span 3.3 pH units from 7.6 in heme oxygenase 1 to 10.9 in peroxidase, and the E(m)s range from -250 mV in peroxidase to 125 mV in Aplysia myoglobin. Proteins with higher in situ ferric aquo-heme pK(a)s tend to have lower E(m)s. Both changes arise from the protein stabilizing a positively charged heme. However, compared with values in solution, the protein shifts the aquo-heme E(m)s more than the pK(a)s. Thus, the protein has a larger effective dielectric constant for the protonation reaction, showing that electron and proton transfers are coupled to different conformational changes that are captured in the MCCE analysis. The calculations reveal a breakdown in the classical continuum electrostatic analysis of pairwise interactions. Comparisons with DFT calculations show that Coulomb's law overestimates the large unfavorable interactions between the ferric water-heme and positively charged groups facing the heme plane by as much as 60%. If interactions with Cu(B) in cytochrome c oxidase and Arg 38 in horseradish peroxidase are not corrected, the pK(a) calculations are in error by as much as 6 pH units. With DFT corrected interactions calculated pK(a)s and E(m)s differ from measured values by less than 1 pH unit or 35 mV, respectively. The in situ aquo-heme pK(a) is important for the function of cytochrome c oxidase since it helps to control the stoichiometry of proton uptake coupled to electron transfer [Song, Michonova-Alexova, and Gunner (2006) Biochemistry 45, 7959-7975].     

Combining conformational flexibility and continuum electrostatics for calculating pK(a)s in proteins

Protein stability and function relies on residues being in their appropriate ionization states at physiological pH. In situ residue pK(a)s also provides a sensitive measure of the local protein environment. Multiconformation continuum electrostatics (MCCE) combines continuum electrostatics and molecular mechanics force fields in Monte Carlo sampling to simultaneously calculate side chain ionization and conformation. The response of protein to charges is incorporated both in the protein dielectric constant (epsilon(prot)) of four and by explicit conformational changes. The pK(a) of 166 residues in 12 proteins was determined. The root mean square error is 0.83 pH units, and >90% have errors of <1 pH units whereas only 3% have errors >2 pH units. Similar results are found with crystal and solution structures, showing that the method's explicit conformational sampling reduces sensitivity to the initial structure. The outcome also changes little with protein dielectric constant (epsilon(prot) 4-20). Multiconformation continuum electrostatics titrations show coupling of conformational flexibility and changes in ionization state. Examples are provided where ionizable side chain position (protein G), Asn orientation (lysozyme), His tautomer distribution (RNase A), and phosphate ion binding (RNase A and H) change with pH. Disallowing these motions changes the calculated pK(a).     

pK values of the ionizable groups of proteins

We have used potentiometric titrations to measure the pK values of the ionizable groups of proteins in alanine pentapeptides with appropriately blocked termini. These pentapeptides provide an improved model for the pK values of the ionizable groups in proteins. Our pK values determined in 0.1 M KCl at 25 degrees C are: 3.67+/-0.03 (alpha-carboxyl), 3.67+/-0.04 (Asp), 4.25+/-0.05 (Glu), 6.54+/-0.04 (His), 8.00+/-0.03 (alpha-amino), 8.55+/-0.03 (Cys), 9.84+/-0.11 (Tyr), and 10.40+/-0.08 (Lys). The pK values of some groups differ from the Nozaki and Tanford (N & T) pK values often used in the literature: Asp (3.67 this work vs. 4.0 N & T); His (6.54 this work vs. 6.3 N & T); alpha-amino (8.00 this work vs. 7.5 N & T); Cys (8.55 this work vs. 9.5 N & T); and Tyr (9.84 this work vs. 9.6 N & T). Our pK values will be useful to those who study pK perturbations in folded and unfolded proteins, and to those who use theory to gain a better understanding of the factors that determine the pK values of the ionizable groups of proteins.     

Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pK(a) calculations

pK(a) calculation methods that are based on finite difference solutions to the Poisson-Boltzmann equation (FDPB) require that energy calculations be performed for a large number of different protonation states of the protein. Normally, the differences between these protonation states are modeled by changing the charges on a few atoms, sometimes the differences are modeled by adding or removing hydrogens, and in a few cases the positions of these hydrogens are optimized locally. We present an FDPB-based pK(a) calculation method in which the hydrogen-bond network is globally optimized for every single protonation state used. This global optimization gives a significant improvement in the accuracy of calculated pK(a) values, especially for buried residues. It is also shown that large errors in calculated pK(a) values are often due to structural artifacts induced by crystal packing. Optimization of the force fields and parameters used in pK(a) calculations should therefore be performed with X-ray structures that are corrected for crystal artifacts.     

Determination of pK(a) values of carboxyl groups in the N-terminal domain of rat CD2: anomalous pK(a) of a glutamate on the ligand-binding surface

The ligand-binding surface of the T-lymphocyte glycoprotein CD2 has an unusually high proportion of charged residues, and ionic interactions are thought to play a significant role in defining the ligand specificity and binding affinity of CD2 with the structurally homologous ligands CD48 (in rodents) and CD58 (in humans). The determination of the electrostatic properties of these proteins can therefore contribute to our understanding of structure-activity relationships for these adhesion complexes that underpin T-cell adhesion to antigen-presenting cells. In this study, we investigated the pH titration behavior of the carboxyl groups of the N-terminal domain of rat CD2 (CD2d1) using the chemical shifts of backbone amide nitrogen-15 ((15)N) and proton NMR resonances, and carboxyl carbon-13 ((13)C) signals. The analysis revealed the presence of a glutamate (Glu41) on the binding surface of rat CD2 with an unusually elevated acidity constant (pK(a) = 6.73) for CD2d1 samples at 1.2 mM concentration. pH titration of CD2d1 at low protein concentration (0.1 mM) resulted in a slight decrease of the measured pK(a) of Glu41 to 6.36. The ionization of Glu41 exhibited reciprocal interactions with a second glutamate (Glu29) in a neighboring location, with both residues demonstrating characteristic biphasic titration behavior of the carboxyl (13)C resonances. Measurements at pH 5.5 of the two-bond deuterium isotope shift for the (13)C carboxyl resonances for Glu41 and Glu29 [(2)DeltaC(delta)(O(epsilon)D) = 0.2 and 0.1 ppm, respectively] were consistent with the assignment of the anomalous pK(a) to Glu41, under the strong influence of Glu29. The characterization of single site mutations of CD2d1 residues Glu41 and Glu29 to glutamine confirmed the anomalous pK(a) for Glu41, and indicated that electrostatic interaction with the Glu29 side chain is a significant contributing influence for this behavior in the wild-type protein. The implications of these observations are discussed with respect to recent structural and functional analyses of the interaction of rat CD2 with CD48. In particular, CD2 Glu41 must be a candidate residue to explain the previously reported strong pH dependence of binding of these two proteins in vitro.     

Cysteine pK(a) values for the bacterial peroxiredoxin AhpC

Salmonella typhimurium AhpC is a founding member of the peroxiredoxin family, a ubiquitous group of cysteine-based peroxidases with high reactivity toward hydrogen peroxide, organic hydroperoxides, and peroxynitrite. For all of the peroxiredoxins, the catalytic cysteine, referred to as the peroxidatic cysteine (C(P)), acts as a nucleophile in attacking the peroxide substrate, forming a cysteine sulfenic acid at the active site. Because thiolates are far stronger nucleophiles than thiol groups, it is generally accepted that cysteine-based peroxidases should exhibit pK(a) values lower than an unperturbed value of 8.3-8.5. In this investigation, several independent approaches were used to assess the pK(a) of the two cysteinyl residues of AhpC. Methods using two different iodoacetamide derivatives yielded unperturbed pK(a) values (7.9-8.7) for both cysteines, apparently due to reactivity with the wrong conformation of C(P) (i.e., locally unfolded and flipped out of the active site), as supported by X-ray crystallographic analyses. A functional pK(a) of 5.94 +/- 0.10 presumably reflecting the titration of C(P) within the fully folded active site was obtained by measuring AhpC competition with horseradish peroxidase for hydrogen peroxide; this value is quite similar to that obtained by analyzing the pH dependence of the epsilon(240) of wild-type AhpC (5.84 +/- 0.02) and similar to those obtained for two typical 2-cysteine peroxiredoxins from Saccharomyces cerevisiae (5.4 and 6.0). Thus, the pK(a) value of AhpC balances the need for a deprotonated thiol (at pH 7, approximately 90% of the C(P) would be deprotonated) with the fact that thiolates with higher pK(a) values are stronger nucleophiles.     

Calculation of protein ionization equilibria with conformational sampling: pK(a) of a model leucine zipper, GCN4 and barnase.

The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pK(a) calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pK(a) estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve-averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pK(a) calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pK(a) predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pK(a) predictions for the model leucine zipper was reduced to as low as approximately 0.25 pK(a) units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pK(a) calculations, are discussed.     

pK(a) calculations along a bacteriorhodopsin molecular dynamics trajectory

Electrostatic calculations of pK(a-values) are reported along a 400 ps molecular dynamics trajectory of bacteriorhodopsin. The sensitivity of calculated pK(a) values to a number of structural factors and factors related to the modelling of the electrostatics are also studied. The results are very sensitive to the choice of internal dielectric constant of the protein (in the interval 2-4). Moreover it is important to include internal water molecules and to average over a long enough portion ( approximately 100 ps) of an equilibrium molecular dynamics trajectory. The internal waters are necessary to get an ion-counter ion complex with the Schiff base and Arg 82 protonated and the aspartic groups (85 and 212) deprotonated. The fluctuations along the MD-trajectory do not change the protonation state of internal residues at neutral pH. However, at other pH values the averaging along a trajectory maybe crucial to get correct protonation states. A relationship is found between the arginine group 82, the aspartic group 85 and the glutamate group 204. Glu 204 is protonated in the ground state but the pK(a) value decreases towards deprotonation when the chromophore isomerizes into the cis state.     

Conformational preferences and pK(a) value of selenocysteine residue

The conformational preferences of the L-selenocysteine (Sec) dipeptides with selenol and selenolate groups (Ac-Sec-NHMe and Ac-Sec(-) -NHMe, respectively) and the apparent (i.e., macroscopic) pK(a) value of the Sec residue have been studied using the dispersion-corrected density functionals M06-2X and B2PLYP-D with the implicit solvation method in the gas phase and in water. In the gas phase, the backbone-to-backbone and/or side chain-to-backbone hydrogen bonds are found to contribute in stabilizing the most preferred conformations for the Sec and Sec(-) residues, as seen for the Cys and Cys(-) residues. However, the polyproline II-like conformations prevail over the conformations with the backbone-to-backbone hydrogen bonds in water because of the weakened hydrogen bonds by the favorable direct interactions between the backbone C?O and H?N groups and water molecules. The Sec and Sec(-) residues are found to adopt more various conformations than the Cys and Cys(-) residues in water, although the most preferred conformations of the neutral and/or anionic forms of the two residues are similar each other in the gas phase and in water. Using the statistically weighted free energies of the Sec and Sec(-) dipeptides in the gas phase and their solvation free energies, the pK(a) value of the Sec residue is estimated to be 5.47 at 25°C, which is in good agreement with the experimental value of 5.43 ± 0.02. It is found that the lower pK(a) value of the selenol side chain for the Sec residue by ～3 units than the thiol side chain for the Cys residue is ascribed to the higher gas-phase acidity of the Sec residue.     

Artificial receptors involved in enolization and pK(a) shifts.

Abiotic receptors used to enolize carbonyl compounds or to shift substrate pK(a) values are reviewed. These systems exhibit disparate frameworks and several approaches to binding and anion stabilization. Detailed emphasis is placed on a bicyclic cyclophane that induces pK(a) shifts in active methylene compounds through NH-pi hydrogen bonding with the resultant enolates.     

Moving a phenol hydroxyl group from the surface to the interior of a protein: effects on the phenol potential and pK(A)

De novo protein design and electrochemistry were used to measure changes in the potential and pK(A) of a phenol when its OH group is moved from a solvent-exposed to a sequestered protein position. A "phenol rotation strategy" was adopted in which phenols, containing a SH in position 4, 3, or 2 relative to the OH group, were bound to a buried protein site. The alpha(3)C protein used here is a tryptophan to cysteine variant of the structurally defined alpha(3)W protein (Dai et al. (2002) J. Am. Chem. Soc. 124, 10952-10953). The protein characteristics of alpha(3)C and the three mercaptophenol-alpha(3)C (MP-alpha(3)C) proteins are shown to be close to those of alpha(3)W. Moreover, the phenol OH group is fully solvent exposed in 4MP-alpha(3)C and more sequestered in 3MP-alpha(3)C and 2MP-alpha(3)C. Here we compare the redox properties of the three mercaptophenols when bound to alpha(3)C and to cysteine free in water. The pK(A) and E(peak) values are essential identical when 4MP is ligated to alpha(3)C relative to when it is free in solution. In contrast, these values are increased in 3MP-alpha(3)C and 2MP-alpha(3)C relative to the solvated compounds. The E(peak) vs pH plots all display a approximately 59 mV/pH unit dependence. We conclude that interactions with the OH group dominate the phenol redox characteristics. In 3MP-alpha(3)C and 2MP-alpha(3)C, hydrogen bonds between the protein and the bound phenols appear to either stabilize the reduced phenol or destabilize the radical, relative to the aqueous buffer, raising the potential by 0.11 and 0.12 V, respectively.     

Remeasuring HEWL pK(a) values by NMR spectroscopy: methods, analysis, accuracy, and implications for theoretical pK(a) calculations

Site-specific pK(a) values measured by NMR spectroscopy provide essential information on protein electrostatics, the pH-dependence of protein structure, dynamics and function, and constitute an important benchmark for protein pK(a) calculation algorithms. Titration curves can be measured by tracking the NMR chemical shifts of several reporter nuclei versus sample pH. However, careful analysis of these curves is needed to extract residue-specific pK(a) values since pH-dependent chemical shift changes can arise from many sources, including through-bond inductive effects, through-space electric field effects, and conformational changes. We have re-measured titration curves for all carboxylates and His 15 in Hen Egg White Lysozyme (HEWL) by recording the pH-dependent chemical shifts of all backbone amide nitrogens and protons, Asp/Glu side chain protons and carboxyl carbons, and imidazole protonated carbons and protons in this protein. We extracted pK(a) values from the resulting titration curves using standard fitting methods, and compared these values to each other, and with those measured previously by 1H NMR (Bartik et al., Biophys J 1994;66:1180–1184). This analysis gives insights into the true accuracy associated with experimentally measured pK(a) values. We find that apparent pK(a) values frequently differ by 0.5–1.0 units depending upon the nuclei monitored, and that larger differences occasionally can be observed. The variation in measured pK(a) values, which reflects the difficulty in fitting and assigning pH-dependent chemical shifts to specific ionization equilibria, has significant implications for the experimental procedures used for measuring protein pK(a) values, for the benchmarking of protein pK(a) calculation algorithms, and for the understanding of protein electrostatics in general.     

Hydrogen bonding markedly reduces the pK of buried carboxyl groups in proteins

The ionizable groups in proteins with the lowest pKs are the carboxyl groups of aspartic acid side-chains. One of the lowest, pK=0.6, is observed for Asp76 in ribonuclease T1. This low pK appeared to result from hydrogen bonds to a water molecule and to the side-chains of Asn9, Tyr11, and Thr91. The results here confirm this by showing that the pK of Asp76 increases to 1.7 in N9A, to 4.0 in Y11F, to 4.2 in T91V, to 4.4 in N9A+Y11F, to 4.9 in N9A+T91V, to 5.9 in Y11F+T91V, and to 6.4 in the triple mutant: N9A+Y11F+T91V. In ribonuclease Sa, the lowest pK=2.4 for Asp33. This pK increases to 3.9 in T56A, which removes the hydrogen bond to Asp33, and to 4.4 in T56V, which removes the hydrogen bond and replaces the -OH group with a -CH(3) group. It is clear that hydrogen bonds are able to markedly lower the pK values of carboxyl groups in proteins. These same hydrogen bonds make large contributions to the conformational stability of the proteins. At pH 7, the stability of D76A ribonuclease T1 is 3.8 kcal mol(-1) less than wild-type, and the stability of D33A ribonuclease Sa is 4.1 kcal mol(-1) less than wild-type. There is a good correlation between the changes in the pK values and the changes in stability. The results suggest that the pK values for these buried carboxyl groups would be greater than 8 in the absence of hydrogen bonds, and that the hydrogen bonds and other interactions of the carboxyl groups contribute over 8 kcal mol(-1) to the stability.     

Combined Density Functional Theory (DFT) and Continuum Calculations of pK(a) in Carbonic Anhydrase

Deprotonation of zinc-bound water in carbonic anhydrase II is the rate-limiting step in the catalysis of carbon dioxide between gas- and water-soluble forms. To understand the factors determining the extent of dissociation, or pK(a), of the zinc-bound water, we apply quantum chemistry calculations to the active site coupled with a continuum model of the surrounding environment. Experimentally determined changes in pK(a) associated with mutations of the active site are well reproduced by this approach. Analysis of the active site structure and charge/dipole values provides evidence that mutations cause changes in both conformation of the active site structure and local polarization, which accounts for the shifts in pK(a). More specifically, the shifts in pK(a) correlate with the dipole moments of the zinc-bound water upon deprotonation. The data further support the conclusion that the distinct pK(a) values found in mutations of the same type, but applied to different sites, result from asymmetric ligation and different electronic environments around the zinc ion.     

Empirical parametrization of pK values for carboxylic acids in proteins using a genetic algorithm

Considerable effort has been devoted to the development of theoretical electrostatic methods to predict the pK values of ionizable residues in proteins. However, predictions appear often to be still at the qualitative or semi-quantitative level. We believe that, with the increasing number experimentally available pK values for proteins of known structure, an alternative approach becomes feasible: the empirical parametrization of the experimental protein pK database. Of course, in the long term, this empirical approach is no substitute for rigorous electrostatic analysis but, in the short term, it may prove to have useful predictive power and it may help to pinpoint the main structural determinants of pK values in proteins. Here we demonstrate the feasibility of the parametrization approach by fitting (using a genetic algorithm as fitting tool) the database for carboxylic acid pK values in proteins on the basis of an empirical equation that takes into account the two following kinds of effects: (1) long-range charge-charge interactions; (2) interactions of the given carboxylic acid group with its environment in the protein, which are described in terms of contributions from the different kind of atoms present in the protein (atomic contributions).