The pKa Cooperative: a collaborative effort to advance structure-based calculations of pKa values and electrostatic effects in proteins

The pK(a) Cooperative (http://www.pkacoop.org) was organized to advance development of accurate and useful computational methods for structure-based calculation of pK(a) values and electrostatic energies in proteins. The Cooperative brings together laboratories with expertise and interest in theoretical, computational, and experimental studies of protein electrostatics. To improve structure-based energy calculations, it is necessary to better understand the physical character and molecular determinants of electrostatic effects. Thus, the Cooperative intends to foment experimental research into fundamental aspects of proteins that depend on electrostatic interactions. It will maintain a depository for experimental data useful for critical assessment of methods for structure-based electrostatics calculations. To help guide the development of computational methods, the Cooperative will organize blind prediction exercises. As a first step, computational laboratories were invited to reproduce an unpublished set of experimental pK(a) values of acidic and basic residues introduced in the interior of staphylococcal nuclease by site-directed mutagenesis. The pK(a) values of these groups are unique and challenging to simulate owing to the large magnitude of their shifts relative to normal pK(a) values in water. Many computational methods were tested in this first Blind Prediction Challenge and critical assessment exercise. A workshop was organized in the Telluride Science Research Center to objectively assess the performance of many computational methods tested on this one extensive data set. This volume of Proteins: Structure, Function, and Bioinformatics introduces the pK(a) Cooperative, presents reports submitted by participants in the Blind Prediction Challenge, and highlights some of the problems in structure-based calculations identified during this exercise.     

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

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

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.     

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.     

High apparent dielectric constant inside a protein reflects structural reorganization coupled to the ionization of an internal Asp

The dielectric properties of proteins are poorly understood and difficult to describe quantitatively. This limits the accuracy of methods for structure-based calculation of electrostatic energies and pK(a) values. The pK(a) values of many internal groups report apparent protein dielectric constants of 10 or higher. These values are substantially higher than the dielectric constants of 2-4 measured experimentally with dry proteins. The structural origins of these high apparent dielectric constants are not well understood. Here we report on structural and equilibrium thermodynamic studies of the effects of pH on the V66D variant of staphylococcal nuclease. In a crystal structure of this protein the neutral side chain of Asp-66 is buried in the hydrophobic core of the protein and hydrated by internal water molecules. Asp-66 titrates with a pK(a) value near 9. A decrease in the far UV-CD signal was observed, concomitant with ionization of this aspartic acid, and consistent with the loss of 1.5 turns of alpha-helix. These data suggest that the protein dielectric constant needed to reproduce the pK(a) value of Asp-66 with continuum electrostatics calculations is high because the dielectric constant has to capture, implicitly, the energetic consequences of the structural reorganization that are not treated explicitly in continuum calculations with static structures.     

High apparent dielectric constants in the interior of a protein reflect water penetration

A glutamic acid was buried in the hydrophobic core of staphylococcal nuclease by replacement of Val-66. Its pK(a) was measured with equilibrium thermodynamic methods. It was 4.3 units higher than the pK(a) of Glu in water. This increase was comparable to the DeltapK(a) of 4.9 units measured previously for a lysine buried at the same location. According to the Born formalism these DeltapK(a) are energetically equivalent to the transfer of a charged group from water to a medium of dielectric constant of 12. In contrast, the static dielectric constants of dry protein powders range from 2 to 4. In the crystallographic structure of the V66E mutant, a chain of water molecules was seen that hydrates the buried Glu-66 and links it with bulk solvent. The buried water molecules have never previously been detected in >20 structures of nuclease. The structure and the measured energetics constitute compelling and unprecedented experimental evidence that solvent penetration can contribute significantly to the high apparent polarizability inside proteins. To improve structure-based calculations of electrostatic effects with continuum methods, it will be necessary to learn to account quantitatively for the contributions by solvent penetration to dielectric effects in the protein interior.     

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.     

Assessing Predictors of Changes in Protein Stability upon Mutation Using Self-Consistency

The ability to predict the effect of mutations on protein stability is important for a wide range of tasks, from protein engineering to assessing the impact of SNPs to understanding basic protein biophysics. A number of methods have been developed that make these predictions, but assessing the accuracy of these tools is difficult given the limitations and inconsistencies of the experimental data. We evaluate four different methods based on the ability of these methods to generate consistent results for forward and back mutations, and examine how this ability varies with the nature and location of the mutation. We find that, while one method seems to outperform the others, the ability of these methods to make accurate predictions is limited.     

Who cares for the protons?

It is tempting to use standard protonation states for the analysis of protein-ligand interactions. Two different pK(a) calculation methods, PROPKA (protein pK(a)) and MCCE (multi conformation continuum electrostatics), were applied to challenge this convenient behavior. As data basis, we selected five recently approved drugs for which structural information of the protein-drug complex is available. We analyzed the pK(a) calculations in terms of a measure termed BIPS (binary protonation states) recently introduced by us. Both methods agree in detecting the majority of the sites with atypical BIPS values. However, when using only one method, some of the atypcial BIPS value would have been missed. Therefore, we recommend using both methods to set such an interpretation on a solid basis.     

Estimating protein-ligand binding free energy: atomic solvation parameters for partition coefficient and solvation free energy calculation

Solvation energy calculation is one of the main difficulties for the estimation of protein-ligand binding free energy and the correct scoring in docking studies. We have developed a new solvation energy estimation method for protein-ligand binding based on atomic solvation parameter (ASP), which has been shown to improve the power of protein-ligand binding free energy predictions. The ASP set, designed to handle both proteins and organic compounds and derived from experimental n-octanol/water partition coefficient (log P) data, contains 100 atom types (united model that treats hydrogen atoms implicitly) or 119 atom types (all-atom model that treats hydrogen atoms explicitly). By using this unified ASP set, an algorithm was developed for solvation energy calculation and was further integrated into a score function for predicting protein-ligand binding affinity. The score function reproduced the absolute binding free energies of a test set of 50 protein-ligand complexes with a standard error of 8.31 kJ/mol. As a byproduct, a conformation-dependent log P calculation algorithm named ASPLOGP was also implemented. The predictive results of ASPLOGP for a test set of 138 compounds were r = 0.968, s = 0.344 for the all-atom model and r = 0.962, s = 0.367 for the united model, which were better than previous conformation-dependent approaches and comparable to fragmental and atom-based methods. ASPLOGP also gave good predictive results for small peptides. The score function based on the ASP model can be applied widely in protein-ligand interaction studies and structure-based drug design.     

Can correct protein models be identified?

The ability to separate correct models of protein structures from less correct models is of the greatest importance for protein structure prediction methods. Several studies have examined the ability of different types of energy function to detect the native, or native-like, protein structure from a large set of decoys. In contrast to earlier studies, we examine here the ability to detect models that only show limited structural similarity to the native structure. These correct models are defined by the existence of a fragment that shows significant similarity between this model and the native structure. It has been shown that the existence of such fragments is useful for comparing the performance between different fold recognition methods and that this performance correlates well with performance in fold recognition. We have developed ProQ, a neural-network-based method to predict the quality of a protein model that extracts structural features, such as frequency of atom-atom contacts, and predicts the quality of a model, as measured either by LGscore or MaxSub. We show that ProQ performs at least as well as other measures when identifying the native structure and is better at the detection of correct models. This performance is maintained over several different test sets. ProQ can also be combined with the Pcons fold recognition predictor (Pmodeller) to increase its performance, with the main advantage being the elimination of a few high-scoring incorrect models. Pmodeller was successful in CASP5 and results from the latest LiveBench, LiveBench-6, indicating that Pmodeller has a higher specificity than Pcons alone.     

Redesigning protein pKa values

The ability to re-engineer enzymatic pH-activity profiles is of importance for industrial applications of enzymes. We theoretically explore the feasibility of re-engineering enzymatic pH-activity profiles by changing active site pK(a) values using point mutations. We calculate the maximum achievable DeltapK(a) values for 141 target titratable groups in seven enzymes by introducing conservative net-charge altering point mutations. We examine the importance of the number of mutations introduced, their distance from the target titratable group, and the characteristics of the target group itself. The results show that multiple mutations at 10A can change pK(a) values up to two units, but that the introduction of a requirement to keep other pK(a) values constant reduces the magnitude of the achievable DeltapK(a). The algorithm presented shows a good correlation with existing experimental data and is available for download and via a web server at http://enzyme.ucd.ie/pKD.     

Progress in the prediction of pKa values in proteins

The pK(a) -cooperative aims to provide a forum for experimental and theoretical researchers interested in protein pK(a) values and protein electrostatics in general. The first round of the pK(a) -cooperative, which challenged computational labs to carry out blind predictions against pK(a) s experimentally determined in the laboratory of Bertrand Garcia-Moreno, was completed and results discussed at the Telluride meeting (July 6-10, 2009). This article serves as an introduction to the reports submitted by the blind prediction participants that will be published in a special issue of PROTEINS: Structure, Function and Bioinformatics. Here, we briefly outline existing approaches for pK(a) calculations, emphasizing methods that were used by the participants in calculating the blind pK(a) values in the first round of the cooperative. We then point out some of the difficulties encountered by the participating groups in making their blind predictions, and finally try to provide some insights for future developments aimed at improving the accuracy of pK(a) calculations.     

Protein flexibility in ligand docking and virtual screening to protein kinases

The main complicating factor in structure-based drug design is receptor rearrangement upon ligand binding (induced fit). It is the induced fit that complicates cross-docking of ligands from different ligand-receptor complexes. Previous studies have shown the necessity to include protein flexibility in ligand docking and virtual screening. Very few docking methods have been developed to predict the induced fit reliably and, at the same time, to improve on discriminating between binders and non-binders in the virtual screening process.We present an algorithm called the ICM-flexible receptor docking algorithm (IFREDA) to account for protein flexibility in virtual screening. By docking flexible ligands to a flexible receptor, IFREDA generates a discrete set of receptor conformations, which are then used to perform flexible ligand-rigid receptor docking and scoring. This is followed by a merging and shrinking step, where the results of the multiple virtual screenings are condensed to improve the enrichment factor. In the IFREDA approach, both side-chain rearrangements and essential backbone movements are taken into consideration, thus sampling adequately the conformational space of the receptor, even in cases of large loop movements.As a preliminary step, to show the importance of incorporating protein flexibility in ligand docking and virtual screening, and to validate the merging and shrinking procedure, we compiled an extensive small-scale virtual screening benchmark of 33 crystal structures of four different protein kinases sub-families (cAPK, CDK-2, P38 and LCK), where we obtained an enrichment factor fold-increase of 1.85±0.65 using two or three multiple experimental conformations. IFREDA was used in eight protein kinase complexes and was able to find the correct ligand conformation and discriminate the correct conformations from the “misdocked” conformations solely on the basis of energy calculation. Five of the generated structures were used in the small-scale virtual screening stage and, by merging and shrinking the results with those of the original structure, we show an enrichment factor fold increase of 1.89±0.60, comparable to that obtained using multiple experimental conformations.Our cross-docking tests on the protein kinase benchmark underscore the necessity of incorporating protein flexibility in both ligand docking and virtual screening. The methodology presented here will be extremely useful in cases where few or no experimental structures of complexes are available, while some binders are known.     

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.     

低能Ti+注入棉花种子的深度-浓度分布的模拟计算

国内许多单位开展了低能离子注入植物种子的实验研究,但在生物诱变机理方面却尚未有定论,其中低能离子注入植物种子的深度-浓度分布成为研究的一个焦点问题.许多人直接用LSS理论得到的TRIM程序,从理论上计算了低能重离子注入植物种子的深度-浓度的分布,却发现计算结果与实验测量结果相差甚远,于是得出LSS理论及TRIM程序不能直接用于计算离子注入非均匀、非致密、有孔道的植物种子这种特殊的靶材料的结论.针对这一研究课题,在目前尚未有更加成熟完善的理论模型和计算方法的前提下,为了便于计算,本文根据植物种子的结构特点,对靶材料棉花种子进行了处理,并对LSS理论进行了修正,用蒙特卡罗方法计算了初始能量为20keV的Ti 注入棉花种子的深度-浓度分布,得到了与实验结果符合较好的分布曲线,初步的确定了低能离子注入植物种子深度-浓度分布的一种理论计算方法.     

Is the prediction of pKa values by constant-pH molecular dynamics being hindered by inherited problems?

In this study, we investigate two factors that can hinder the performance of constant-pH molecular dynamics methods in predicting protein pK(a) values, using hen egg white lysozyme as a test system. The first factor is related to the molecular definition and pK(a) value of model compounds in the Poisson-Boltzmann framework. We address this by defining the model compound as a molecular fragment with an associated pK(a) value that is calibrated against experimental data, which results in a decrease of 0.12 units in pK(a) errors. The second addressed factor is the possibility that detrimental structural distortions are being introduced in the simulations by the underlying molecular mechanics force field. This issue is investigated by analyzing how the gradual structural rearrangements affect the predicted pK(a) values. The two GROMOS force fields studied here (43A1 and 53A6) yield good pK(a) predictions, although a time-dependent performance is observed: 43A1 performs better after a few nanoseconds of structural reorganization (pK(a) errors of ~0.45), while 53A6 gives the best prediction right at the first nanosecond (pK(a) errors of 0.42). These results suggest that the good performance of constant-pH molecular dynamics methods could be further improved if these force field limitations were overcome.     

A buried lysine that titrates with a normal pKa: role of conformational flexibility at the protein-water interface as a determinant of pKa values

Previously we reported that Lys, Asp, and Glu residues at positions 66 and 92 in staphylococcal nuclease (SNase) titrate with pK(a) values shifted by up to 5 pK(a) units in the direction that promotes the neutral state. In contrast, the internal Lys-38 in SNase titrates with a normal pK(a). The crystal structure of the L38K variant shows that the side chain of Lys-38 is buried. The ionizable moiety is approximately 7 A from solvent and ion paired with Glu-122. This suggests that the pK(a) value of Lys-38 is normal because the energetic penalty for dehydration is offset by a favorable Coulomb interaction. However, the pK(a) of Lys-38 was also normal when Glu-122 was replaced with Gln or with Ala. Continuum electrostatics calculations were unable to reproduce the pK(a) of Lys-38 unless the protein was treated with an artificially high dielectric constant, consistent with structural reorganization being responsible for the normal pK(a) value of Lys-38. This reorganization must be local because circular dichroism and NMR spectroscopy indicate that the L38K protein is native-like under all conditions studied. In molecular dynamics simulations, the ion pair between Lys-38 and Glu-122 is unstable. The simulations show that a minor rearrangement of a loop is sufficient to allow penetration of water to the amino moiety of Lys-38. This illustrates both the important roles of local flexibility and water penetration as determinants of pK(a) values of ionizable groups buried near the protein-water interface, and the challenges faced by structure-based pK(a) calculations in reproducing these effects.     

Salt effects on ionization equilibria of histidines in myoglobin

The salt dependence of histidine pK(a) values in sperm whale and horse myoglobin and in histidine-containing peptides was measured by (1)H-NMR spectroscopy. Structure-based pK(a) calculations were performed with continuum methods to test their ability to capture the effects of solution conditions on pK(a) values. The measured pK(a) of most histidines, whether in the protein or in model compounds, increased by 0.3 pH units or more between 0.02 M and 1.5 M NaCl. In myoglobin two histidines (His(48) and His(36)) exhibited a shallower dependence than the average, and one (His(113)) showed a steeper dependence. The (1)H-NMR data suggested that the salt dependence of histidine pK(a) values in the protein was determined primarily by the preferential stabilization of the charged form of histidine with increasing salt concentrations rather than by screening of electrostatic interactions. The magnitude and salt dependence of interactions between ionizable groups were exaggerated in pK(a) calculations with the finite-difference Poisson-Boltzmann method applied to a static structure, even when the protein interior was treated with arbitrarily high dielectric constants. Improvements in continuum methods for calculating salt effects on pK(a) values will require explicit consideration of the salt dependence of model compound pK(a) values used for reference in the calculations.