A parameterized, continuum electrostatic model for predicting protein pKa values

Recognizing the limits of trying to achieve chemical accuracy for pK(a) calculations with a purely electrostatic model, we include empirical corrections into the Poisson-Boltzmann solver macroscopic electrostatics with atomic detail (Bashford, Biochemistry 1990;29:10219-10225), to improve the reliability and accuracy of the model. The total number of parameters is kept to a minimum to maximize the robustness of the model for compounds outside of the fitting dataset. The parameters are based on: (a) the electrostatic interaction between functional groups close to the titratable site, (b) the electrostatic work required to desolvate the residue, and (c) the site-to-site interactions. These interactions are straightforward to calculate once the electrostatic field has been solved for each residue using the linearized Poisson-Boltzmann equation and are assumed to be linearly related to the intrinsic pK(a). Two hundred and eighty-six residues from 30 proteins are used to determine the empirical parameters, which result in a root mean square error (RMSE) of 0.70 for the entire set. Eight proteins with 46 experimentally known values were excluded from the parameterization to test the model. This test set had a RMSE of 1.08. We show that the parameterized model improves the results over other models, although like other models the error is strongly correlated with the degree to which a residue is buried. The parameters themselves indicate that local effects are most important for determining the pK(a), whereas site-to-site interactions are found to be less significant.     

Predicting protein homocysteinylation targets based on dihedral strain energy and pKa of cysteines

A multitude of complex diseases have been linked to elevated homocysteine levels; however, till date there is no plausible explanation for a single amino acid's involvement in so many diseases. Since homocysteine is a reactive thiol amino acid and the majority of plasma homocysteine is protein thiol bound, we hypothesized that homocysteine might bind to accessible cysteine residues in target proteins, thereby modulating its structure or function or both. The parameters that dictate homocysteine-protein interaction are not well understood, and the few known homocysteine binding proteins were identified by a candidate protein approach. In this study, we identified potential homocysteine interacting proteins based on cysteine content, solvent accessibility of cysteine residues, and dihedral strain energies and pKa of these cysteines. Pathway mapping of the cysteine-rich proteins revealed that proteins in the coagulation cascade, notch receptor-mediated signaling, LDL endocytosis, programmed cell death, and extracellular matrix proteins were significantly over-represented with cysteine-rich proteins, and we believe that homocysteine has a high probability to bind to proteins in these pathways. In fact, several clinical studies have implicated high homocysteine levels to be associated with diseases like thrombosis, neural tube defects, and so forth, which result from dysfunction of one or more of the proteins identified in our study. Further, we successfully validated our prediction parameters on the proteins that have already been experimentally shown to bind homocysteine, and our structural analysis argues a plausible explanation for these prior reported protein interactions with homocysteine that could not be previously explained.     

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.     

Bayesian model aggregation for ensemble‐based estimates of protein pKa values

This article investigates an ensemble‐based technique called Bayesian Model Averaging (BMA) to improve the performance of protein amino acid pK_a predictions. Structure‐based pK_a calculations play an important role in the mechanistic interpretation of protein structure and are also used to determine a wide range of protein properties. A diverse set of methods currently exist for pK_a prediction, ranging from empirical statistical models to ab initio quantum mechanical approaches. However, each of these methods are based on a set of conceptual assumptions that can effect a model's accuracy and generalizability for pK_a prediction in complicated biomolecular systems. We use BMA to combine eleven diverse prediction methods that each estimate pK_a values of amino acids in staphylococcal nuclease. These methods are based on work conducted for the pK_a Cooperative and the pK_a measurements are based on experimental work conducted by the García‐Moreno lab. Our cross‐validation study demonstrates that the aggregated estimate obtained from BMA outperforms all individual prediction methods with improvements ranging from 45 to 73% over other method classes. This study also compares BMA's predictive performance to other ensemble‐based techniques and demonstrates that BMA can outperform these approaches with improvements ranging from 27 to 60%. This work illustrates a new possible mechanism for improving the accuracy of pK_a prediction and lays the foundation for future work on aggregate models that balance computational cost with prediction accuracy. Proteins 2014; 82:354–363. © 2013 Wiley Periodicals, Inc.     

TSAR,a new graph–theoretical approach to computational modeling of protein side‐chain flexibility: Modeling of ionization properties of proteins

A new graph–theoretical approach called thermodynamic sampling of amino acid residues (TSAR) has been elaborated to explicitly account for the protein side chain flexibility in modeling conformation‐dependent protein properties. In TSAR, a protein is viewed as a graph whose nodes correspond to structurally independent groups and whose edges connect the interacting groups. Each node has its set of states describing conformation and ionization of the group, and each edge is assigned an array of pairwise interaction potentials between the adjacent groups. By treating the obtained graph as a belief‐network—a well‐established mathematical abstraction—the partition function of each node is found. In the current work we used TSAR to calculate partition functions of the ionized forms of protein residues. A simplified version of a semi‐empirical molecular mechanical scoring function, borrowed from our Lead Finder docking software, was used for energy calculations. The accuracy of the resulting model was validated on a set of 486 experimentally determined pK_a values of protein residues. The average correlation coefficient (R) between calculated and experimental pK_a values was 0.80, ranging from 0.95 (for Tyr) to 0.61 (for Lys). It appeared that the hydrogen bond interactions and the exhaustiveness of side chain sampling made the most significant contribution to the accuracy of pK_a calculations. Proteins 2011; © 2011 Wiley‐Liss, Inc.     

Calculating pKa values in the cAMP-dependent protein kinase: the effect of conformational change and ligand binding

The conformational change observed upon ligand binding and phosphorylation for the cAMP-dependent protein kinase (protein kinase A-PKA) is of high importance for the regulation of its activity. We calculate pKa values and net charges for 18 3D structures of PKA in various conformations and liganded states to examine the role of electrostatics in ligand binding and activation. We find that the conformational change of PKA takes place without any significant net proton uptake/release at all pH values, thus indicating that PKA has evolved to reduce any pH-dependent barriers to the conformational motion. We furthermore find that the binding of ligands induces large changes in the net charge of PKA at most pH values, but significantly, we find that the net charge difference at physiological pH is close to zero, thus indicating that the active-site pKa values have been preorganized for substrate binding. We are unable to unequivocally resolve the identity of the groups responsible for determining the pH-activity profile of PKA but speculate that the titration of Lys 168 or the titration of ATP itself could be responsible for the loss of activity at high pH values. Finally, we examine the effect of point mutations on the pKa values of the PKA catalytic residues and find these to be relatively insensitive to both noncharge-altering and charge-altering mutations.     

A Fast Method for pKa Determination by Capillary Electrophoresis

In this article, we describe a recently developed capillary‐electrophoresis method for the determination of acidity constants and compare it with other existing methods. The new method is based on the use of an internal standard (compound similar in nature and pK_a value to the analyte), and offers several benefits, since it has all the advantages of capillary electrophoresis. In addition, it is very fast, because the exact measure of the pH of the separation electrolytes is not needed, and only a few electrophoretic runs are required to perform a pK_a determination. The acidity constants of some monoprotic weak acids and bases were determined by this fast method, yielding a very good agreement with literature values.     

In silico pKa Prediction and ADME Profiling

Improving the ADME profile of drug candidates is a critical step in lead optimization, and because pK_a affects most ADME properties such as lipophilicity, solubility, and metabolism, it is extremely advantageous to predict pK_a in order to guide the design of new compounds. However, accurately (<0.5 log units) predicting pK_a by empirical methods can be challenging especially for novel series, because of lack of knowledge on determinants of pK_a (steric effects, ring effects, H‐bonding, etc.), and because of limited experimental data on the effects of specific chemical groups on the ionization of an atom. To address these issues, we recently developed the computational package MoKa, which integrates graphical and command line tools designed for computational and medicinal chemists to predict the pK_a values of organic compounds. Here, we present the major issues considered when we developed MoKa, such as the accurate selection of training data, the fundamentals of the methodology (which has also been extended to predict protein pK_a), the treatment of multiprotic compounds, and the selection of the dominant tautomer for the calculation. Last, we illustrate some specific applications of MoKa to predict solubility, lipophilicity, and metabolism.     

Prediction of titration properties of structures of a protein derived from molecular dynamics trajectories.

This paper explores the dependence of the molecular dynamics (MD) trajectory of a protein molecule on the titration state assigned to the molecule. Four 100-ps MD trajectories of bovine pancreatic trypsin inhibitor (BPTI) were generated, starting from two different structures, each of which was held in two different charge states. The two starting structures were the X-ray crystal structure and one of the solution structures determined by NMR, and the charge states differed only in the ionization state of N terminus. Although it is evident that the MD simulations were too short to sample fully the equilibrium distribution of structures in each case, standard Poisson-Boltzmann titration state analysis of the resulting configurations shows general agreement between the overall titration behavior of the protein and the charge state assumed during MD simulation: at pH 7, the total net charge of the protein resulting from the titration analysis is consistently lower for the protein with the N terminus assumed to be neutral than for the protein with the N terminus assumed to be charged. For most of the ionizable residues, the differences in the calculated pKaS among the four trajectories are statistically negligible and remain in good agreement with the data obtained by crystal structure titration and by experiment. The exceptions include the N terminus, which responds directly to the change of its imposed charge; the C terminus, which in the NMR structure interacts strongly with the former; and a few other residues (Arg 1, Glu 7, Tyr 35, and Arg 42) whose pKaS reflect the initial structure and the limited trajectory lengths. This study illustrates the importance of the careful assignment of protonation states at the start of MD simulations and points to the need for simulation methods that allow for the variation of the protonation state in the calculation of equilibrium properties.     

pKa predictions for proteins,RNAs, and DNAs with the Gaussian dielectric function using DelPhi pKa

We developed a Poisson‐Boltzmann based approach to calculate the values of protein ionizable residues (Glu, Asp, His, Lys and Arg), nucleotides of RNA and single stranded DNA. Two novel features were utilized: the dielectric properties of the macromolecules and water phase were modeled via the smooth Gaussian‐based dielectric function in DelPhi and the corresponding electrostatic energies were calculated without defining the molecular surface. We tested the algorithm by calculating values for more than 300 residues from 32 proteins from the PPD dataset and achieved an overall RMSD of 0.77. Particularly, the RMSD of 0.55 was achieved for surface residues, while the RMSD of 1.1 for buried residues. The approach was also found capable of capturing the large shifts of various single point mutations in staphylococcal nuclease (SNase) from ‐cooperative dataset, resulting in an overall RMSD of 1.6 for this set of pKa's. Investigations showed that predictions for most of buried mutant residues of SNase could be improved by using higher dielectric constant values. Furthermore, an option to generate different hydrogen positions also improves predictions for buried carboxyl residues. Finally, the calculations on two RNAs demonstrated the capability of this approach for other types of biomolecules. Proteins 2015; 83:2186–2197. © 2015 Wiley Periodicals, Inc.     

构建适用于蛋白质家族分类的相似性网络

图聚类用于蛋白质分类问题可以获得较好结果,其前提是将蛋白质之间复杂的相互关系转化为适当的相似性网络作为图聚类分类的输入数据。本文提出一种基于BLAST检索的相似性网络构建方法,从目标蛋白质序列出发,通过若干轮次的BLAST检索逐步从数据库中提取与目标蛋白质直接或间接相关的序列,构成关联集。关联集中序列之间的相似性关系即相似性网络,可作为图聚类算法的分类依据。对Pfam数据库中依直接相似关系难以正确分类的蛋白质的计算表明,按本文方法构建的相似性网络取得了比较满意的结果。     

Simplified methods for pKa and acid pH-dependent stability estimation in proteins: removing dielectric and counterion boundaries

Much computational research aimed at understanding ionizable group interactions in proteins has focused on numerical solutions of the Poisson-Boltzmann (PB) equation, incorporating protein exclusion zones for solvent and counterions in a continuum model. Poor agreement with measured pKas and pH-dependent stabilities for a (protein, solvent) relative dielectric boundary of (4,80) has lead to the adoption of an intermediate (20,80) boundary. It is now shown that a simple Debye-Huckel (DH) calculation, removing both the low dielectric and counterion exclusion regions associated with protein, is equally effective in general pKa calculations. However, a broad-based discrepancy to measured pH-dependent stabilities is maintained in the absence of ionizable group interactions in the unfolded state. A simple model is introduced for these interactions, with a significantly improved match to experiment that suggests a potential utility in predicting and analyzing the acid pH-dependence of protein stability. The methods are applied to the relative pH-dependent stabilities of the pore-forming domains of colicins A and N. The results relate generally to the well-known preponderance of surface ionizable groups with solvent-mediated interactions. Although numerical PB solutions do not currently have a significant advantage for overall pKa estimations, development based on consideration of microscopic solvation energetics in tandem with the continuum model could combine the large deltapKas of a subset of ionizable groups with the overall robustness of the DH model.     

Surface‐based protein binding pocket similarity

Protein similarity comparisons may be made on a local or global basis and may consider sequence information or differing levels of structural information. We present a local three‐dimensional method that compares protein binding site surfaces in full atomic detail. The approach is based on the morphological similarity method which has been widely applied for global comparison of small molecules. We apply the method to all‐by‐all comparisons two sets of human protein kinases, a very diverse set of ATP‐bound proteins from multiple species, and three heterogeneous benchmark protein binding site data sets. Cases of disagreement between sequence‐based similarity and binding site similarity yield informative examples. Where sequence similarity is very low, high pocket similarity can reliably identify important binding motifs. Where sequence similarity is very high, significant differences in pocket similarity are related to ligand binding specificity and similarity. Local protein binding pocket similarity provides qualitatively complementary information to other approaches, and it can yield quantitative information in support of functional annotation. Proteins 2011; © 2011 Wiley‐Liss, Inc.     

Empirical relationships between protein structure and carboxyl pKa values in proteins

Relationships between protein structure and ionization of carboxyl groups were investigated in 24 proteins of known structure and for which 115 aspartate and 97 glutamate pK(a) values are known. Mean pK(a) values for aspartates and glutamates are < or = 3.4 (+/-1.0) and 4.1 (+/-0.8), respectively. For aspartates, mean pK(a) values are 3.9 (+/-1.0) and 3.1 (+/-0.9) in acidic (pI < 5) and basic (pI > 8) proteins, respectively, while mean pK(a) values for glutamates are approximately 4.2 for acidic and basic proteins. Burial of carboxyl groups leads to dispersion in pK(a) values: pK(a) values for solvent-exposed groups show narrow distributions while values for buried groups range from < 2 to 6.7. Calculated electrostatic potentials at the carboxyl groups show modest correlations with experimental pK(a) values and these correlations are not improved by including simple surface-area-based terms to account for the effects of desolvation. Mean aspartate pK(a) values decrease with increasing numbers of hydrogen bonds but this is not observed at glutamates. Only 10 pK(a) values are > 5.5 and most are found in active sites or ligand-binding sites. These carboxyl groups are buried and usually accept no more than one hydrogen bond. Aspartates and glutamates at the N-termini of helices have mean pK(a) values of 2.8 (+/-0.5) and 3.4 (+/-0.6), respectively, about 0.6 units less than the overall mean values.     

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.     

Arginine: Its pKa value revisited

Using complementary approaches of potentiometry and NMR spectroscopy, we have determined that the equilibrium acid dissociation constant (pK_a value) of the arginine guanidinium group is 13.8 ± 0.1. This is substantially higher than that of ∼12 often used in structure-based electrostatics calculations and cited in biochemistry textbooks. The revised intrinsic pK_a value helps explains why arginine side chains in proteins are always predominantly charged, even at pH values as great as 10. The high pK_a value also reinforces the observation that arginine side chains are invariably protonated under physiological conditions of near neutral pH. This occurs even when the guanidinium moiety is buried in a hydrophobic micro-environment, such as that inside a protein or a lipid membrane, thought to be incompatible with the presence of a charged group.