Benchmarking pK<Subscript>a</Subscript> prediction methods for Lys115 in acetoacetate decarboxylase

Three different pK_a prediction methods were used to calculate the pK_a of Lys115 in acetoacetate decarboxylase (AADase): the empirical method PROPKA, the multiconformation continuum electrostatics (MCCE) method, and the molecular dynamics/thermodynamic integration (MD/TI) method with implicit solvent. As expected, accurate pK_a prediction of Lys115 depends on the protonation patterns of other ionizable groups, especially the nearby Glu76. However, since the prediction methods do not explicitly sample the protonation patterns of nearby residues, this must be done manually. When Glu76 is deprotonated, all three methods give an incorrect pK_a value for Lys115. If protonated, Glu76 is used in an MD/TI calculation, the pK_a of Lys115 is predicted to be 5.3, which agrees well with the experimental value of 5.9. This result agrees with previous site-directed mutagenesis studies, where the mutation of Glu76 (negative charge when deprotonated) to Gln (neutral) causes no change in K_m, suggesting that Glu76 has no effect on the pK_a shift of Lys115. Thus, we postulate that the pK_a of Glu76 is also shifted so that Glu76 is protonated (neutral) in AADase.

Open image in new window

Graphical abstract Simulated abundances of protonated species as pH is varied

     

The effects of pK(a) tuning on the thermodynamics and kinetics of folding: design of a solvent-shielded carboxylate pair at the a-position of a coiled-coil

The tuning of the pK_a of ionizable residues plays a critical role in various protein functions, such as ligand-binding, catalysis, and allostery. Proteins harness the free energy of folding to position ionizable groups in highly specific environments that strongly affect their pK_a values. To investigate the interplay among protein folding kinetics, thermodynamics, and pK_a modulation, we introduced a pair of Asp residues at neighboring interior positions of a coiled-coil. A single Asp residue was replaced for an Asn side chain at the a-position of the coiled-coil from GCN4, which was also crosslinked at the C-terminus via a flexible disulfide bond. The thermodynamic and kinetic stability of the system was measured by circular dichroism and stopped-flow fluorescence as a function of pH and concentration of guanidine HCl. Both sets of data are consistent with a two-state equilibrium between fully folded and unfolded forms. Distinct pK_a values of 6.3 and 5.35 are assigned to the first and second protonation of the Asp pair; together they represent an energetic difference of 5 kcal/mol relative to the protonation of two Asp residues with unperturbed pK_a values. Analysis of the rate data as a function of pH and denaturant concentration allowed calculation of the kinetic constants for the conformational transitions of the peptide with the Asp residues in the doubly protonated, singly protonated, and unprotonated forms. The doubly and singly protonated forms fold rapidly, and a ϕ-value analysis shows that their contribution to folding occurs subsequent to the transition state ensemble for folding. By contrast, the doubly charged state shows a reduced rate of folding and a ϕ-value near 0.5 indicative of a repulsive interaction, and possibly also heterogeneity in the transition state ensemble.     

194 Energetics and dynamics of adenine protonation in HIV-1’s dimerization initiation site

HIV-1, like most retroviruses, packages two homologous copies of its RNA genome. The two RNA strands are non-covalently linked near their 5’ends. The proposed dimerisation initiation site is a 35-nucleotide (nt) stem loop capable of forming loop–loop interactions (a kissing dimer) via its highly conserved 6-nt loop palindrome. In a structural transformation affected by temperature, salt concentration, and by the HIV-1 nucleocapsid protein, the initial, kinetically-stable kissing dimer (KD) converts to a thermodynamically-stable extended dimer. It has been suggested that this in vitro observed rearrangement is associated with the in vivo viral genome maturation. Mihailescu et al. demonstrated enhanced rearrangement dynamics triggered by the protonation of a specific adenine residue at the genome sequence location 272 (A272) (Mihailescu, 2004). They suggested that the local environment of A272 caused its N1 atom’s pK_a to shift upwards (more basic) by approximately 2.5 pH units when compared with 5’-adenosine monophosphate (5’-AMP). In this work, we investigated the dynamics and energetics of protonating A272’s N1 atom with explicit solvent molecular dynamics (MD) simulations, Thermodynamic Integration (TI), and the Poisson–Boltzmann equation (PB). Two initial structures were used, an NMR solution structure (PDB ID: 2D19) and an X-ray crystal structure (PDB ID: 1XP7), where A272 was found inside (NMR) and outside (X-ray) the helical axis respectively. MD simulations showed when A272 started,it remained inside the KD’s axis, and when outside, it attempted to insert itself within the axis. Calculated pK_a shifts obtained from solving the PB equation were approximately?+?2.2 and?+?1.0 pH units for the NMR and X-ray structures respectively; while TI calculations performed with the NMR structure yielded a shift of approximately?+?2.5 pH units. Our simulations confirmed the strong influence of A272's local environment on its calculated pK_a when inserted in the helical axis. A272's N1 atom was approximately 200 times more likely protonated when compared with 5'-AMP. Also, protonated A272s were more energetically favoured in the kissing dimer versus an isolated monomer due to a diminished positive electrostatic potential near A272, while in the dimer. Overall, our computational investigations affirm the experimental suggestion that A272 in the kissing dimer is more likely to be protonated at physiological conditions and this protonation may trigger the structural rearrangement of the initial kissing dimer to form the extended dimer.     

Graphical analysis of pH-dependent properties of proteins predicted using PROPKA

Background  

Charge states of ionizable residues in proteins determine their pH-dependent properties through their pK_a values. Thus, various theoretical methods to determine ionization constants of residues in biological systems have been developed. One of the more widely used approaches for predicting pK_a values in proteins is the PROPKA program, which provides convenient structural rationalization of the predicted pK_a values without any additional calculations.     

Improving the accuracy of protein pKa calculations: Conformational averaging versus the average structure

Several methods for including the conformational flexibility of proteins in the calculation of titration curves are compared. The methods use the linearized Poisson-Boltzmann equation to calculate the electrostatic free energies of solvation and are applied to bovine pancreatic trypsin inhibitor (BPTI) and hen egg-white lysozyme (HEWL). An ensemble of conformations is generated by a molecular dynamics simulation of the proteins with explicit solvent. The average titration curve of the ensemble is calculated in three different ways: an average structure is used for the pK_a calculation; the electrostatic interaction free energies are averaged and used for the pK_a calculation; and the titration curve for each structure is calculated and the curves are averaged. The three averaging methods give very similar results and improve the pK_a values to approximately the same degree. This suggests, in contrast to implications from other work, that the observed improvement of pK_a values in the present studies is due not to averaging over an ensemble of structures, but rather to the generation of a single properly averaged structure for the pK_a calculation. Proteins 33:145–158, 1998. © 1998 Wiley-Liss, Inc.     

Prediction of Secondary Ionization of the Phosphate Group in Phosphotyrosine Peptides

A computational approach, based on a continuum molecular electrostatics model, for the calculation of the pK_a values of secondary ionization of the phosphate group in phenyl phosphate derivatives is described. The method uses the ESP atomic charges of the mono-anionic and di-anionic forms of the ionizable phosphate group, computed with the use of the density functional method, and applies a new concept of the model group, being the reference state for the pK_a calculations. Both conformational flexibility and tautomeric degrees of freedom are taken into account in the calculations. The method was parameterized using experimentally available pK_a values of four derivatives of phenyl phosphates, and phosphotyrosine. Subsequently this parameterization was used to predict pK_a of the phosphate group in a short peptide Gly-Gly-Tyr(P)-Ala, and in a longer peptide consisting of 12 residues, the latter in water, and in a complex with a protein—phospholipase. The agreement between the computed and the experimental pK_a values is better than ±0.3 pH units for the optimized solute dielectric constant of 11–13. This approach is promising and its extension to other phospho-amino acids is in progress.     

Tuning of Thioredoxin Redox Properties by Intramolecular Hydrogen Bonds

Thioredoxin-like proteins contain a characteristic C-x-x-C active site motif and are involved in a large number of biological processes ranging from electron transfer, cellular redox level maintenance, and regulation of cellular processes. The mechanism for deprotonation of the buried C-terminal active site cysteine in thioredoxin, necessary for dissociation of the mixed-disulfide intermediate that occurs under thiol/disulfide mediated electron transfer, is not well understood for all thioredoxin superfamily members. Here we have characterized a 8.7 kD thioredoxin (BC3987) from Bacillus cereus that unlike the typical thioredoxin appears to use the conserved Thr8 side chain near the unusual C-P-P-C active site to increase enzymatic activity by forming a hydrogen bond to the buried cysteine. Our hypothesis is based on biochemical assays and thiolate pK_a titrations where the wild type and T8A mutant are compared, phylogenetic analysis of related thioredoxins, and QM/MM calculations with the BC3987 crystal structure as a precursor for modeling of reduced active sites. We suggest that our model applies to other thioredoxin subclasses with similar active site arrangements.     

Evaluating amber force fields using computed NMR chemical shifts

NMR chemical shifts can be computed from molecular dynamics (MD) simulations using a template matching approach and a library of conformers containing chemical shifts generated from ab initio quantum calculations. This approach has potential utility for evaluating the force fields that underlie these simulations. Imperfections in force fields generate flawed atomic coordinates. Chemical shifts obtained from flawed coordinates have errors that can be traced back to these imperfections. We use this approach to evaluate a series of AMBER force fields that have been refined over the course of two decades (ff94, ff96, ff99SB, ff14SB, ff14ipq, and ff15ipq). For each force field a series of MD simulations are carried out for eight model proteins. The calculated chemical shifts for the ¹H, ¹⁵N, and ¹³C^a atoms are compared with experimental values. Initial evaluations are based on root mean squared (RMS) errors at the protein level. These results are further refined based on secondary structure and the types of atoms involved in nonbonded interactions. The best chemical shift for identifying force field differences is the shift associated with peptide protons. Examination of the model proteins on a residue by residue basis reveals that force field performance is highly dependent on residue position. Examination of the time course of nonbonded interactions at these sites provides explanations for chemical shift differences at the atomic coordinate level. Results show that the newer ff14ipq and ff15ipq force fields developed with the implicitly polarized charge method perform better than the older force fields.     

Conformational sampling and polarization of Asp26 in pKa calculations of thioredoxin

Thioredoxin is a protein that has been used as model system by various computational methods to predict the pK_a of aspartate residue Asp26 which is 3.5 units higher than a solvent exposed one (eg, Asp20). Here, we use extensive atomistic molecular dynamics simulations of two different protonation states of Asp26 in combination with conformational analysis based on RMSD clustering and principle component analysis to identify representative conformations of the protein in solution. For each conformation, the Gibbs free energy of proton transfer between Asp26 and Asp20, which is fully solvated in a loop region of the protein, is calculated with the Amber99sb force field in alchemical transformations. The varying polarization of the two residues in different molecular environments and protonation states is described by Hirshfeld-I (HI) atomic charges obtained from the averaged polarized electron density. Our results show that the Gibbs free energy of proton transfer is dependent on the protein conformation, the proper sampling of the neighboring Lys57 residue orientations and on water molecules entering the hydrophobic cavity upon deprotonating Asp26. The inclusion of the polarization of both aspartate residues in the free energy cycle by HI atomic charges corrects the results from the non-polarizable force field and reproduces the experimental ΔpK_a value of Asp26.     

Catalytic subunit of cAMP-dependent protein kinase: Electrostatic features and peptide recognition

The electrostatic field was calculated for the mammalian cAMP-dependent protein kinase (PKA) catalytic subunit (C-subunit) complexed with a 20-residue peptide from a heat stable protein kinase inhibitor (PKI: 5–24). The electrostatic field was also calculated for the C-subunit complexed with a modeled heptapeptide substrate that has been used extensively in structure/function studies for the C-subunit. Perturbations in the electrostatic free energy were calculated when single ionizable active site residues were mutated to alanine. These perturbations in electrostatic free energy were correlated to changes in the binding energy measured in a charge-to-alanine scan of the homologous yeast C-subunit by M. J. Zoller and C. S. Gibbs [(1991) Journal of Biological Chemistry, Vol. 266, pp. 8923–8931; C. S. Gibbs and M. J. Zoller (1991) Biochemistry, Vol. 30, p. 22]. This analysis indicated that the substrate binding parameters primarily depend on electrostatic interactions between a substrate or inhibitor and the C-subunit. Amino acid replacements that led to large perturbations in the electrostatic field are listed in the text. pK_a shifts were also calculated for the substrate's phosphate accepting atom, the serine hydroxyl oxygen, when the active site ionizable residues were changed to structurally similar uncharged amino acids. The theoretical mutation of three active site residues caused large shifts in this parameter: E91Q, D166N, and D184N. The calculated pK_a shifts for these mutants indicate that the rate of phosphotransfer should be markedly reduced in these cases. This prediction has been experimentally confirmed for the D166N mutant. The correlation between calculated electrostatic free energy changes and measured binding energy, and pK_a shifts with phosphotransfer for C-subunit mutants were within experimental error of the measurements. The calculations of electrostatic energy and ΔpK_a have identified previously unconsidered active site residues in the mammalian C-subunit that contribute to binding energy and phosphotransfer. © 1996 John Wiley & Sons, Inc.     

Microscopic mechanisms that govern the titration response and pKa values of buried residues in staphylococcal nuclease mutants

To probe the microscopic mechanisms that govern the titration behavior of buried ionizable groups, microsecond explicit solvent molecular dynamics simulations are carried out for several mutants of Staphylococcal nuclease using both fixed charge and polarizable force fields. While the ionization of Asp 66, Glu 66, and Lys 125 lead to enhanced structural fluctuations and partial unfolding of adjacent α‐helical regions, the ionization of Lys 25 causes local unfolding of adjacent β sheets. Using the sampled conformational ensembles, good agreement with experimental pK_a values is obtained with Poisson–Boltzmann calculations using a protein dielectric constant of 2–4 for V66D/E; slightly larger dielectric constants are needed for Lys mutants especially L25K, suggesting that structural responses beyond microseconds are involved in ionization of Lys 25. Overall, the set of unbiased simulations provides insights into the spatial and temporal scales of protein and solvent motions that dictate the diverse titration behaviors of buried protein residues. Proteins 2017; 85:268–281. © 2016 Wiley Periodicals, Inc.     

pK(a) values for the unfolded state under native conditions explain the pH-dependent stability of PGB1

Understanding the role of electrostatics in protein stability requires knowledge of these interactions in both the folded and unfolded states. Electrostatic interactions can be probed experimentally by characterizing ionization equilibria of titrating groups, parameterized as pK_a values. However, pK_a values of the unfolded state are rarely accessible under native conditions, where the unfolded state has a very low population. Here, we report pK_a values under nondenaturing conditions for two unfolded fragments of the protein G B1 domain that mimic the unfolded state of the intact protein. pK_a values were determined for carboxyl groups by monitoring their pH-dependent ¹³C chemical shifts. Monte Carlo simulations using a Gaussian chain model provide corrections for changes in electrostatic interactions that arise from fragmentation of the protein. Most pK_a values for the unfolded state agree well with model values, but some residues show significant perturbations that can be rationalized by local electrostatic interactions. The pH-dependent stability was calculated from the experimental pK_a values of the folded and unfolded states and compared to experimental stability data. The use of experimental pK_a values for the unfolded state results in significantly improved agreement with experimental data, as compared to calculations based on model data alone.     

Highly perturbed pKa values in the unfolded state of hen egg white lysozyme

The majority of pK_a values in protein unfolded states are close to the amino acid model pK_a values, thus reflecting the weak intramolecular interactions present in the unfolded ensemble of most proteins. We have carried out thermal denaturation measurements on the WT and eight mutants of HEWL from pH 1.5 to pH 11.0 to examine the unfolded state pK_a values and the pH dependence of protein stability for this enzyme. The availability of accurate pK_a values for the folded state of HEWL and separate measurements of mutant-induced effects on the folded state pK_a values, allows us to estimate the pK_a values of seven acidic residues in the unfolded state of HEWL. Asp-48 and Asp-66 display pK_a values of 2.9 and 3.1 in our analysis, thus representing the most depressed unfolded state pK_a values observed to date. We observe a strong correlation between the folded state pK_a values and the unfolded state pK_a values of HEWL, thus suggesting that the unfolded state of HEWL possesses a large degree of native state characteristics.     

pH-Dependent Conformational Changes in Proteins and Their Effect on Experimental pKas: The Case of Nitrophorin 4

The acid-base behavior of amino acids is an important subject of study due to their prominent role in enzyme catalysis, substrate binding and protein structure. Due to interactions with the protein environment, their pK_as can be shifted from their solution values and, if a protein has two stable conformations, it is possible for a residue to have different “microscopic”, conformation-dependent pK_a values. In those cases, interpretation of experimental measurements of the pK_a is complicated by the coupling between pH, protonation state and protein conformation. We explored these issues using Nitrophorin 4 (NP4), a protein that releases NO in a pH sensitive manner. At pH 5.5 NP4 is in a closed conformation where NO is tightly bound, while at pH 7.5 Asp30 becomes deprotonated, causing the conformation to change to an open state from which NO can easily escape. Using constant pH molecular dynamics we found two distinct microscopic Asp30 pK_as: 8.5 in the closed structure and 4.3 in the open structure. Using a four-state model, we then related the obtained microscopic values to the experimentally observed “apparent” pK_a, obtaining a value of 6.5, in excellent agreement with experimental data. This value must be interpreted as the pH at which the closed to open population transition takes place. More generally, our results show that it is possible to relate microscopic structure dependent pKa values to experimentally observed ensemble dependent apparent pK_as and that the insight gained in the relatively simple case of NP4 can be useful in several more complex cases involving a pH dependent transition, of great biochemical interest.     

Inductive effects in amino acids and peptides: Ionization constants and tryptophan fluorescence

Although inductive effects in organic compounds are known to influence chemical properties such as ionization constants, their specific contribution to the properties/behavior of amino acids and functional groups in peptides remains largely unexplored. In this study we developed a computationally economical algorithm for ab initio calculation of the magnitude of inductive effects for non-aromatic molecules. The value obtained by the algorithm is called the Inductive Index and we observed a high correlation (R² = 0.9427) between our calculations and the pK_a values of the alpha-amino groups of amino acids with non-aromatic side-chains. Using a series of modified amino acids, we also found similarly high correlations (R² > 0.9600) between Inductive Indexes and two wholly independent chemical properties: i) the pK_a values of ionizable side-chains and, ii) the fluorescence response of the indole group of tryptophan. After assessing the applicability of the method of calculation at the amino acid level, we extended our study to tryptophan-containing peptides and established that inductive contributions of neighboring side-chains are transmitted through peptide bonds. We discuss possible contributions to the study of proteins.     

Protonation state of Asp30 exerts crucial influence over surface loop rearrangements responsible for NO release in nitrophorin 4

pK(a) values of ionizable residues were calculated for the crystal structures describing the pH and NO binding dependant conformations of nitrophorin 4, a pH sensitive NO carrier heme protein. Comparison of resultant H-bonding patterns allowed the identification of the amino acids that take part in signaling pH change. We carried out MD simulations to show that the protonation state of Asp30, buried in the closed conformation, is crucial for maintaining the tight packed conformation of the closed form of the complex - presenting a model for the functional decrease of NO binding affinity of nitrophorins at physiological pH.     

Backbone relaxation coupled to the ionization of internal groups in proteins: a self-guided Langevin dynamics study

Pathways of structural relaxation triggered by ionization of internal groups in staphylococcal nuclease (SNase) were studied through multiple self-guided Langevin dynamics (SGLD) simulations. Circular dichroism, steady-state Trp fluorescence, and nuclear magnetic resonance spectroscopy have shown previously that variants of SNase with internal Glu, Asp, and Lys at positions 66 or 92, and Arg at position 66, exhibit local reorganization or global unfolding when the internal ionizable group is charged. Except for Arg-66, these internal ionizable groups have unusual pK_a values and are neutral at physiological pH. The structural trends observed in the simulations are in general agreement with experimental observations. The I92D variant, which unfolds globally upon ionization of Asp-92, in simulations often exhibits extensive hydration of the protein core, and sometimes also significant perturbations of the β-barrel. In the crystal structure of the V66R variant, the β1 strand from the β-barrel is domain-swapped; in the simulations, the β1 strand is sometimes partially released. The V66K variant, which in solutions shows reorganization of six residues at the C-terminus of helix α1 and perturbations in the β-barrel structure, exhibits fraying of three residues of helix α1 in one simulation, and perturbations and partial unfolding of three β-strands in a few other simulations. In sharp contrast, very small structural changes were observed in simulations of the wild-type protein. The simulations indicate that charging of internal groups frequently triggers penetration of water into the protein interior. The pK_a values of Asp-92 and Arg-66 calculated with continuum methods on SGLD-relaxed structures reached the normal values in most simulations. Detailed analysis of accuracy and performance of SGLD demonstrates that SGLD outperforms LD in sampling of alternative protein conformations without loss of the accuracy and level of detail characteristic of regular LD.     

Importance of polarization effect in the study of metalloproteins: Application of polarized protein specific charge scheme in predicting the reduction potential of azurin

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein‐specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild‐type azurin (ΔE_cal) were calculated and compared with experimental values. The ΔE_cal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein‐specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations. Proteins 2014; 82:2209–2219. © 2014 Wiley Periodicals, Inc.     

Electrostatics of folded and unfolded bovine β-lactoglobulin

We report on electrophoretic, spectroscopic, and computational studies aimed at clarifying, at atomic resolution, the electrostatics of folded and unfolded bovine β-lactoglobulin (BLG) with a detailed characterization of the specific aminoacids involved. The procedures we used involved denaturant gradient gel electrophoresis, isoelectric focusing, electrophoretic titration curves, circular dichroism and fluorescence spectra in the presence of increasing concentrations of urea (up to 8 M), electrostatics computations and low-mode molecular dynamics. Discrepancy between electrophoretic and spectroscopic evidence suggests that changes in mobility induced by urea are not just the result of changes in gyration radius upon unfolding. Electrophoretic titration curves run across a pH range of 3.5–9 in the presence of urea suggest that more than one aminoacid residue may have anomalous pK _a value in native BLG. Detailed computational studies indicate a shift in pKa of Glu44, Glu89, and Glu114, mainly due to changes in global and local desolvation. For His161, the formation of hydrogen bond(s) could add up to desolvation contributions. However, since His161 is at the C terminus, the end-effect associated to the solvated form strongly influences its pK _a value with extreme variation between crystal structures on one side and NMR or low-mode molecular dynamics structures on the other. The urea concentration effective in BLG unfolding depends on pH, with higher stability of the protein at lower pH.