Predicting protein pKa by environment similarity

A statistical method to predict protein pK_a has been developed by using the 3D structure of a protein and a database of 434 experimental protein pK_a values. Each pK_a in the database is associated with a fingerprint that describes the chemical environment around an ionizable residue. A computational tool, MoKaBio, has been developed to identify automatically ionizable residues in a protein, generate fingerprints that describe the chemical environment around such residues, and predict pK_a from the experimental pK_a values in the database by using a similarity metric. The method, which retrieved the pK_a of 429 of the 434 ionizable sites in the database correctly, was crossvalidated by leave‐one‐out and yielded root mean square error (RMSE) = 0.95, a result that is superior to that obtained by using the Null Model (RMSE 1.07) and other well‐established protein pK_a prediction tools. This novel approach is suitable to rationalize protein pK_a by comparing the region around the ionizable site with similar regions whose ionizable site pK_a is known. The pK_a of residues that have a unique environment not represented in the training set cannot be predicted accurately, however, the method offers the advantage of being trainable to increase its predictive power. Proteins 2009. © 2009 Wiley‐Liss, Inc.     

pH dependence of conformational fluctuations of the protein backbone

Proton binding equilibria (pK_a values) of ionizable groups in proteins are exquisitely sensitive to their microenvironments. Apparent pK_a values measured for individual ionizable residues with NMR spectroscopy are actually population‐weighted averages of the pK_a in different conformational microstates. NMR spectroscopy experiments with staphylococcal nuclease were used to test the hypothesis that pK_a values of surface Glu and Asp residues are affected by pH‐sensitive fluctuations of the backbone between folded and locally unfolded conformations. ¹⁵N spin relaxation studies showed that as the pH decreases from the neutral into the acidic range the amplitudes of backbone fluctuations in the ps‐ns timescale increase near carboxylic residues. Hydrogen exchange experiments suggested that backbone conformational fluctuations promoted by decreasing pH also reflect slower local or sub‐global unfolding near carboxylic groups. This study has implications for structure‐based pK_a calculations: (1) The timescale of the backbone's response to ionization events in proteins can range from ps to ms, and even longer; (2) pH‐sensitive fluctuations of the backbone can be localized to both the segment the ionizable residue is attached to or the one that occludes the ionizable group; (3) Structural perturbations are not necessarily propagated through Coulomb interactions; instead, local fluctuations appear to be coupled through the co‐operativity inherent to elements of secondary structure and to networks of hydrogen bonds. These results are consistent with the idea that local conformational fluctuations and stabilities are important determinants of apparent pK_a values of ionizable residues in proteins. Proteins 2014; 82:3132–3143. © 2014 Wiley Periodicals, Inc.     

Carbon-13 NMR studies of the lysine side chains of calmodulin and its proteolytic fragments

ThepH-titration and dynamic behaviour of the seven lysine side chains in bovine calmodulin were studied by carbon-13 NMR. The amino groups of the calcium saturated protein and its proteolytic fragments TR₁C(1–75) and TR₂C (78–148) were dimethylated with carbon-13 labeled formaldehyde; this modification did not alter the protein's structure or its ability to activate the enzyme cyclic nucleotide phosphodiesterase. Tentative assignments for 5 out of the 7 dimethyl lysine resonances could be obtained by comparing spectra of the fully and partially modified protein, with those of the proteolytic fragments. ThepK_a values measured for calcium saturated calmodulin ranged between 9.5 (Lys 75) and 10.2 (Lys 13); two residues (Lys 94 and Lys 13) showed a biphasic titration curve suggesting their possible involvement in ion-pairs. The dynamic behavior of the lysine side chains was deduced from spin lattice relaxation measurements. All side chains were flexible and this was not influenced by the removal of calcium, or the addition of the calmodulin antagonist trifluoperazine. The latter data suggest that the lysine side chains are not directly involved in calmodulin's target binding sites.     

Reductive methylation andpKa determination of the lysine side chains in calbindin D9k

The Lys residues in the 75-residue Ca²⁺-binding protein calbindin D_9k were reductively methylated with¹³C-enriched formaldehyde. The possible structural effects resulting from the chemical modification were critically investigated by comparing two-dimensional NMR spectra and the exchange rates of some of the amide protons of the native and the modified protein. Our results show that the protein retains its structure even though 10 Lys out of a total of 75 amino acid residues were modified. In the Ca²⁺- and apo-forms of the protein, the¹³C-methylated Lys residues can be detected with high sensitivity and resolution using two-dimensional (¹H,¹³C)-heteronuclear multiple quantum coherence (HMQC) NMR spectroscopy. ThepKa values of the individual Lys residues in Ca²⁺-calbindin D_9k and apo-calbindin D_9k were obtained by combiningpH titration experiments and (¹H,¹³C)-HMQC NMR spectroscopy. Each Lys residue in the Ca²⁺- and apo-forms of calbindin D_9k has a uniquepKa value. The LyspKa values in the calcium protein range from 9.3 to 10.9, while those in the apo-protein vary between 9.7 and 10.7. Although apo-calbindin D_9k has a very similar structure compared to Ca²⁺-calbindin D_9k, the removal of two Ca²⁺ ions from the protein leads to an increase of thepKa values of the Lys residues.     

Comparative studies of the interaction of human and bovine platelet factor 4 with heparin using histidine NMR resonances as spectroscopic probes

ThepK _a values of His-38 and His-50 of the heparin-binding protein, bovine platelet factor 4, are 5.6 and 6.5, respectively, as determined by¹H NMR spectroscopy. The¹H NMR resonance of His-38 of bovine platelet factor 4 which exhibits the lowerpK _a value is perturbed upon heparin binding to a greater degree than the resonance of His-50. Human platelet factor 4 contains the homologous residues His-23 and His-35. ThepK _a values of the two histidine residues of human platelet factor 4 are 5.3 and 6.4. The¹H NMR resonance of the histidine of human platelet factor 4 exhibiting the lowerpK _a value also is perturbed upon heparin binding to a greater degree than the histidine resonance exhibiting the higherpK _a , thereby suggesting comparable heparin-protein interactions in bovine and human platelet factor 4.     

Electrostatic contributions to the stability of the GCN4 leucine zipper structure

Ion pairs are ubiquitous in X-ray structures of coiled coils, and mutagenesis of charged residues can result in large stability losses. By contrast, pK_a values determined by NMR in solution often predict only small contributions to stability from charge interactions. To help reconcile these results we used triple-resonance NMR to determine pK_a values for all groups that ionize between pH 1 and 13 in the 33 residue leucine zipper fragment, GCN4p. In addition to the native state we also determined comprehensive pK_a values for two models of the GCN4p denatured state: the protein in 6 M urea, and unfolded peptide fragments of the protein in water. Only residues that form ion pairs in multiple X-ray structures of GCN4p gave large pK_a differences between the native and denatured states. Moreover, electrostatic contributions to stability were not equivalent for oppositely charged partners in ion pairs, suggesting that the interactions between a charge and its environment are as important as those within the ion pair. The pH dependence of protein stability calculated from NMR-derived pK_a values agreed with the stability profile measured from equilibrium urea-unfolding experiments as a function of pH. The stability profile was also reproduced with structure-based continuum electrostatic calculations, although contributions to stability were overestimated at the extremes of pH. We consider potential sources of errors in the calculations, and how pK_a predictions could be improved. Our results show that although hydrophobic packing and hydrogen bonding have dominant roles, electrostatic interactions also make significant contributions to the stability of the coiled coil.     

Intrinsically disordered inhibitor of glutamine synthetase is a functional protein with random‐coil‐like pKa values

The sequential action of glutamine synthetase (GS) and glutamate synthase (GOGAT) in cyanobacteria allows the incorporation of ammonium into carbon skeletons. In the cyanobacterium Synechocystis sp. PCC 6803, the activity of GS is modulated by the interaction with proteins, which include a 65‐residue‐long intrinsically disordered protein (IDP), the inactivating factor IF7. This interaction is regulated by the presence of charged residues in both IF7 and GS. To understand how charged amino acids can affect the binding of an IDP with its target and to provide clues on electrostatic interactions in disordered states of proteins, we measured the pK_a values of all IF7 acidic groups (Glu32, Glu36, Glu38, Asp40, Asp58, and Ser65, the backbone C‐terminus) at 100 mM NaCl concentration, by using NMR spectroscopy. We also obtained solution structures of IF7 through molecular dynamics simulation, validated them on the basis of previous experiments, and used them to obtain theoretical estimates of the pK_a values. Titration values for the two Asp and three Glu residues of IF7 were similar to those reported for random‐coil models, suggesting the lack of electrostatic interactions around these residues. Furthermore, our results suggest the presence of helical structure at the N‐terminus of the protein and of conformational changes at acidic pH values. The overall experimental and in silico findings suggest that local interactions and conformational equilibria do not play a role in determining the electrostatic features of the acidic residues of IF7.     

Structural origins of high apparent dielectric constants experienced by ionizable groups in the hydrophobic core of a protein

The side chains of Lys66, Asp66, and Glu66 in staphylococcal nuclease are fully buried and surrounded mainly by hydrophobic matter, except for internal water molecules associated with carboxylic oxygen atoms. These ionizable side chains titrate with pK_a values of 5.7, 8.8, and 8.9, respectively. To reproduce these pK_a values with continuum electrostatics calculations, we treated the protein with high dielectric constants. We have examined the structural origins of these high apparent dielectric constants by using NMR spectroscopy to characterize the structural response to the ionization of these internal side chains. Substitution of Val66 with Lys66 and Asp66 led to increased conformational fluctuations of the microenvironments surrounding these groups, even under pH conditions where Lys66 and Asp66 are neutral. When Lys66, Asp66, and Glu66 are charged, the proteins remain almost fully folded, but resonances for a few backbone amides adjacent to the internal ionizable residues are broadened. This suggests that the ionization of the internal groups promotes a local increase in dynamics on the intermediate timescale, consistent with either partial unfolding or increased backbone fluctuations of helix 1 near residue 66, or, less likely, with increased fluctuations of the charged side chains at position 66. These experiments confirm that the high apparent dielectric constants reported by internal Lys66, Asp66, and Glu66 reflect localized changes in conformational fluctuations without incurring detectable global structural reorganization. To improve structure-based pK_a calculations in proteins, we will need to learn how to treat this coupling between ionization of internal groups and local changes in conformational fluctuations explicitly.     

Carboxyl pKa Values and Acid Denaturation of BBL

The protein BBL undergoes structural transitions and acid denaturation between pH 1.2 and 8.0. Using NMR spectroscopy, we measured the pK_a values of all the carboxylic residues in this pH range. We employed ¹³C direct-detection two-dimensional IPAP (in-phase antiphase) CACO NMR spectroscopy to monitor the ionization state of different carboxylic groups and demonstrated its advantages over other NMR techniques in measuring pK_a values of carboxylic residues. The two residues Glu161 and Asp162 had significantly lowered pK_a values, showing that these residues are involved in a network of stabilizing electrostatic interactions, as is His166. The other carboxylates had unperturbed values. The pH dependence of the free energy of denaturation was described quantitatively by the ionizations of those three residues of perturbed pK_a, and, using thermodynamic cycles, we could calculate their pK_as in the native and denatured states as well as the equilibrium constants for denaturation of the different protonation states. We also measured ¹³C^α chemical shifts of individual residues as a function of pH. These shifts sense structural transitions rather than ionizations, and they titrated with pH consistent with the change in equilibrium constant for denaturation. Kinetic measurements of the folding of BBL E161Q indicated that, at pH 7, the stabilizing interactions with Glu161 are formed mainly in the transition state. We also found that local interactions still exist in the acid-denatured state of BBL, which attenuate somewhat the flexibility of the acid-denatured state.     

Physicochemical properties of hydroxylated polychlorinated biphenyls aid in predicting their interactions with rat sulfotransferase 1A1 (rSULT1A1)

Hydroxylated metabolites of polychlorinated biphenyls (OHPCBs) interact with rat sulfotransferase 1A1 (rSULT1A1) as substrates and inhibitors. Previous studies have shown that there are complex and incompletely understood structure–activity relationships governing the interaction of rSULT1A1 with these molecules. Furthermore, modification of the enzyme with glutathione disulfide (GSSG) results in the conversion of some OHPCBs from inhibitors to substrates. We have now examined estimated values for the acid-dissociation constant (K_a) and the octanol–water distribution coefficient (D), as well as experimentally determined dissociation constants for enzyme complexes, to assist in the prediction of interactions of OHPCBs with rSULT1A1. Under reducing conditions, initial velocities for rSULT1A1-catalyzed sulfation exhibited a positive correlation with pK_a and a negative correlation with log D of the OHPCBs. IC₅₀ values of inhibitory OHPCBs decreased with decreasing pK_a values for both the glutathione (GSH)-pretreated and GSSG-pretreated forms of rSULT1A1. Comparison of GSH- and GSSG-pretreated forms of rSULT1A1 with respect to binding of OHPCB in the presence and absence of adenosine 3′,5′-diphosphate (PAP) revealed that the dissociation constants with the two redox states of the enzyme were similar for each OHPCB. Thus, pK_a and log D values are useful in predicting the binding of OHPCBs to the two redox forms of rSULT1A1 as well as the rates of sulfation of those OHPCBs that are substrates. However, the differences in substrate specificity for OHPCBs that are seen with changes in redox status of the enzyme are not directly related to specific structural effects of individual OHPCBs within inhibitory enzyme–PAP–OHPCB complexes.     

Mutation of non-conserved amino acids surrounding catalytic site to shift pH optimum of Bacillus circulans xylanase

Improvement of enzyme function by engineering pH dependence of enzymatic activity is of importance for industrial application of Bacillus circulans xylanases. Target mutation sites were selected by structural alignment between B. circulans xylanase and other xylanases having different pH optima. We selected non-conserved mutant sites within 8 Å from the catalytic residues, to see whether these residues have some role in modulating pK_as of the catalytic residues. We hypothesized that the non-conserved residues which may not have any role in enzyme catalysis might perturb pK_as of the catalytic residues. Change in pK_a of a titratable group due to change in electrostatic potential of a mutation was calculated and the change in pH optimum was predicted from the change in pK_a of the catalytic residues. Our strategy is proved to be useful in selection of promising mutants to shift the pH optimum of the xylanases towards desired side.     

Mitochondrial protein acetylation as a cell-intrinsic,evolutionary driver of fat storage: Chemical and metabolic logic of acetyl-lysine modifications

Abstract

Hormone systems evolved over 500 million years of animal natural history to motivate feeding behavior and convert excess calories to fat. These systems produced vertebrates, including humans, who are famine-resistant but sensitive to obesity in environments of persistent overnutrition. We looked for cell-intrinsic metabolic features, which might have been subject to an evolutionary drive favoring lipogenesis. Mitochondrial protein acetylation appears to be such a system. Because mitochondrial acetyl-coA is the central mediator of fuel oxidation and is saturable, this metabolite is postulated to be the fundamental indicator of energy excess, which imprints a memory of nutritional imbalances by covalent modification. Fungal and invertebrate mitochondria have highly acetylated mitochondrial proteomes without an apparent mitochondrially targeted protein lysine acetyltransferase. Thus, mitochondrial acetylation is hypothesized to have evolved as a nonenzymatic phenomenon. Because the pK_a of a nonperturbed Lys is 10.4 and linkage of a carbonyl carbon to an ? amino group cannot be formed with a protonated Lys, we hypothesize that acetylation occurs on residues with depressed pK_a values, accounting for the propensity of acetylation to hit active sites and suggesting that regulatory Lys residues may have been under selective pressure to avoid or attract acetylation throughout animal evolution. In addition, a shortage of mitochondrial oxaloacetate under ketotic conditions can explain why macronutrient insufficiency also produces mitochondrial hyperacetylation. Reduced mitochondrial activity during times of overnutrition and undernutrition would improve fitness by virtue of resource conservation. Micronutrient insufficiency is predicted to exacerbate mitochondrial hyperacetylation. Nicotinamide riboside and Sirt3 activity are predicted to relieve mitochondrial inhibition.     

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.     

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.     

Electrostatics of cysteine residues in proteins: Parameterization and validation of a simple model

One of the most popular and simple models for the calculation of pK_as from a protein structure is the semi‐macroscopic electrostatic model MEAD. This model requires empirical parameters for each residue to calculate pK_as. Analysis of current, widely used empirical parameters for cysteine residues showed that they did not reproduce expected cysteine pK_as; thus, we set out to identify parameters consistent with the CHARMM27 force field that capture both the behavior of typical cysteines in proteins and the behavior of cysteines which have perturbed pK_as. The new parameters were validated in three ways: (1) calculation across a large set of typical cysteines in proteins (where the calculations are expected to reproduce expected ensemble behavior); (2) calculation across a set of perturbed cysteines in proteins (where the calculations are expected to reproduce the shifted ensemble behavior); and (3) comparison to experimentally determined pK_a values (where the calculation should reproduce the pK_a within experimental error). Both the general behavior of cysteines in proteins and the perturbed pK_a in some proteins can be predicted reasonably well using the newly determined empirical parameters within the MEAD model for protein electrostatics. This study provides the first general analysis of the electrostatics of cysteines in proteins, with specific attention paid to capturing both the behavior of typical cysteines in a protein and the behavior of cysteines whose pK_a should be shifted, and validation of force field parameters for cysteine residues. Proteins 2012. © 2012 Wiley Periodicals, Inc.     

The 40s Ω-loop plays a critical role in the stability and the alkaline conformational transition of cytochrome<Emphasis Type="Italic"> c</Emphasis>

The structural and redox properties of a non-covalent complex reconstituted upon mixing two non-contiguous fragments of horse cytochrome c, the residues 1–38 heme-containing N-fragment with the residues 57–104 C-fragment, have been investigated. With respect to native cyt c, the complex lacks a segment of 18 residues, corresponding, in the native protein, to an omega ()-loop region. The fragment complex shows compact structure, native-like -helix content but a less rigid atomic packing and reduced stability with respect to the native protein. Structural heterogeneity is observed at pH 7.0, involving formation of an axially misligated low-spin species and consequent partial displacement of Met80 from the sixth coordination position of the heme-iron. Spectroscopic data suggest that a lysine (located in the Met80-containing loop, namely Lys72, Lys73, or Lys79) replaces the methionine residue. The residues 1–38/57–104 fragment complex shows an unusual biphasic alkaline titration characterized by a low (pK_a1=6.72) and a high pK_a-associated state transition (pK_a2=8.56); this behavior differs from that of native cyt c, which shows a monophasic alkaline transition (pK_a=8.9). The data indicate that the 40s -loop plays an important role in the stability of cyt c and in ensuring a correct alkaline conformational transition of the protein.     

Self‐guided Langevin dynamics study of regulatory interactions in NtrC

Multiple self‐guided Langevin dynamics (SGLD) simulations were performed to examine structural and dynamical properties of the receiver domain of nitrogen regulatory protein C (NtrC^r). SGLD and MD simulations of the phosphorylated active form structure suggest a mostly stable but broad structural ensemble of this protein. The finite difference Poisson–Boltzmann calculations of the pK_a values of the active site residues suggest an increase in the pK_a of His‐84 on phosphorylation of Asp‐54. In SGLD simulations of the phosphorylated active form with charged His‐84, the average position of the regulatory helix α4 is found closer to the starting structure than in simulations with the neutral His‐84. To model the transition pathway, the phosphate group was removed from the simulations. After 7 ns of simulations, the regulatory helix α4 was found approximately halfway between positions in the NMR structures of the active and inactive forms. Removal of the phosphate group stimulated loss of helix α4, suggesting that the pathway of conformational transition may involve partial unfolding mechanism. The study illustrates the potential utility of the SGLD method in studies of the coupling between ligand binding and conformational transitions. Proteins 2009. © 2009 Wiley‐Liss, Inc.     

Error estimation and global fitting in transverse-relaxation dispersion experiments to determine chemical-exchange parameters

Off-resonance effects can introduce significant systematic errors in R₂ measurements in constant-time Carr-Purcell-Meiboom-Gill (CPMG) transverse relaxation dispersion experiments. For an off-resonance chemical shift of 500 Hz, ¹⁵N relaxation dispersion profiles obtained from experiment and computer simulation indicated a systematic error of ca. 3%. This error is three- to five-fold larger than the random error in R₂ caused by noise. Good estimates of total R₂ uncertainty are critical in order to obtain accurate estimates in optimized chemical exchange parameters and their uncertainties derived from χ² minimization of a target function. Here, we present a simple empirical approach that provides a good estimate of the total error (systematic + random) in ¹⁵N R₂ values measured for the HIV protease. The advantage of this empirical error estimate is that it is applicable even when some of the factors that contribute to the off-resonance error are not known. These errors are incorporated into a χ² minimization protocol, in which the Carver–Richards equation is used fit the observed R₂ dispersion profiles, that yields optimized chemical exchange parameters and their confidence limits. Optimized parameters are also derived, using the same protein sample and data-fitting protocol, from ¹H R₂ measurements in which systematic errors are negligible. Although ¹H and ¹⁵N relaxation profiles of individual residues were well fit, the optimized exchange parameters had large uncertainties (confidence limits). In contrast, when a single pair of exchange parameters (the exchange lifetime, τ_ex, and the fractional population, p_a), were constrained to globally fit all R₂ profiles for residues in the dimer interface of the protein, confidence limits were less than 8% for all optimized exchange parameters. In addition, F-tests showed that quality of the fits obtained using τ_ex, p_a as global parameters were not improved when these parameters were free to fit the R₂ profiles of individual residues. Finally, nearly the same optimized global τ_ex, p_a values were obtained, when the ¹H and ¹⁵N data sets for residues in the dimer interface, were fit independently; the difference in optimized global parameters, ca. 10%, was of marginal significance according to the F-test.