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).     

Calculation of pK(a) in proteins with the microenvironment modulated-screened coulomb potential

The MM-SCP has been applied to predict pK(a) values of titratable residues in wild type and mutants of staphylococcal nuclease (SNase). The calculations were based on crystal structures made available by the Garcia-Moreno Laboratory. In the mutants, mostly deeply buried hydrophobic residues were replaced with ionizable residues, and thus their pK(a) values could be measured and calculated using various methods. The data set used here consisted of a set of WT SNase for which His pK(a) at several ionic strengths had been measured, a set of mutants for which measured pK(a) were available and a set of 11 mutants for which the measured pK(a) were not known at the time of calculation. For this latter set, blind predictions were submitted to the protein pK(a) cooperative, 2009 workshop at Telluride, where the results of the blind predictions were discussed (the RMSD of the submitted set was 1.10 pH units). The calculations on the structures with known pK(a) indicated that in addition to weaknesses of the method, structural issues were observed that led to larger errors (>1) in pK(a) predictions. For example, different crystallographic conditions or steric clashes can lead to differences in the local environment around the titratable residue, which can produce large differences in the calculated pK(a) . To gain further insight into the reliability of the MM-SCP, pK(a) of an extended set of 54 proteins belonging to several structural classes were carried out. Here some initial results from this study are reported to help place the SNase results in the appropriate context.     

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.     

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.     

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.     

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.     

Zinc-binding proteins (ligands) in brains of severely zinc-deficient rats

Previous studies from this laboratory reported the presence of a metallothionein-like protein in brain with an apparent estimated molecular weight of 13,000–15,000 daltons. The synthesis of this protein, which incorporates large quantity of cysteine, is stimulated following administration of zinc and copper and is blocked by actinomycin D. In this study, we report that the synthesis of this metallothionein-like protein is considerably lower in brains of severely zinc-deficient rats in comparison with pair-fed orad libitum fed groups. Furthermore, incubation of partially purified metallothionein-like protein with⁶⁵Zn and chromatography on DEAE A-25 Sephadex produced similar elution patterns in the three experimental groups. However, the extent of binding of⁶⁵Zn to the metallothionein-like protein from the zinc-deficient rats was significantly (p<0.05) lower than the control groups. On the other hand, the total concentration of zinc in brains of zinc deficient rats did not vary from control groups. Since the synthesis of this metallothionein-like protein is reduced by zinc deficiency and is stimulated following administration of zinc, we postulate that the free pool of zinc may regulate the synthesis of its binding protein in the brain.     

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.     

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.     

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.     

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.     

Role of the protein side-chain fluctuations on the strength of pair-wise electrostatic interactions: comparing experimental with computed pK(a)s

The effect of the protein side-chain fluctuations on the strength of electrostatic interactions was studied. The effect was modeled on 7 different crystal structures on the same enzyme as well as on 20 molecular dynamics snapshot structures. It was shown that the side-chain flexibility affects predominantly the magnitude of the strong pair-wise interactions, that is, the pair-wise interaction among ion pairs, and practically does not affect the interactions with the rest of the protein. This was used to suggest a correction function that should be applied to the original pair-wise electrostatic interaction to mimic the effects of the fluctuations. The procedure is applied on three ion pairs identified in lysozyme. It was shown that sampling different side-chain rotamers and modifying the strength of the pair-wise interaction energies makes calculated pK(a)s less sensitive to the fluctuations of the structure and improves the prediction accuracy.     

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.     

pK(a) calculations for class C beta-lactamases: the role of Tyr-150

The Poisson-Boltzmann method was used to compute the pK(a) values of titratable residues in a set of class C beta-lactamases. In these calculations, the pK(a) of the phenolic group of residue Tyr150 is the only one to stand out with an abnormally low value of 8.3, more than one pK(a) unit lower than the measured reference value for tyrosine in solution. Other important residues of the catalytic pocket, such as the conserved Lys67, Lys315, His314, and Glu272 (hydrogen-bonded to the ammonium group of Lys315), display normal protonation states at neutral pH. pK(a) values were also computed in catalytically impaired beta-lactamase mutants. Comparisons between the relative k(cat) values and the Tyr150 pK(a) value in these mutants revealed a striking correlation. In active enzymes, this pK(a) value is always lower than the solution reference value while it is close to normal in inactive enzymes. These results thus support the hypothesis that the phenolate form of Tyr150 is responsible for the activation of the nucleophilic serine. The possible roles of Lys67 and Lys315 during catalysis are also discussed.     

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.     

Post-translational chemical modification(s) of proteins.

1. The role played by the modification of protein in determining its fate is reported by us. 2. Post-translational modifications such as acetylation, phosphorylation, sulfation, methylation, hydroxylation, ADP-ribosylation, maturation, amidation, carboxylation, adenylylation, glycosylation, ubiquitination, and prenylation are extensively reviewed. 3. Each post-translational modification's significance and its role played in biological function(s) is summarized in the general discussion and the conclusion's remark is directed at the problems left to solve (e.g. post-translational modification reactions in recombinant protein in modern genetic engineering).