A fast and accurate computational approach to protein ionization

We report a very fast and accurate physics-based method to calculate pH-dependent electrostatic effects in protein molecules and to predict the pK values of individual sites of titration. In addition, a CHARMm-based algorithm is included to construct and refine the spatial coordinates of all hydrogen atoms at a given pH. The present method combines electrostatic energy calculations based on the Generalized Born approximation with an iterative mobile clustering approach to calculate the equilibria of proton binding to multiple titration sites in protein molecules. The use of the GBIM (Generalized Born with Implicit Membrane) CHARMm module makes it possible to model not only water-soluble proteins but membrane proteins as well. The method includes a novel algorithm for preliminary refinement of hydrogen coordinates. Another difference from existing approaches is that, instead of monopeptides, a set of relaxed pentapeptide structures are used as model compounds. Tests on a set of 24 proteins demonstrate the high accuracy of the method. On average, the RMSD between predicted and experimental pK values is close to 0.5 pK units on this data set, and the accuracy is achieved at very low computational cost. The pH-dependent assignment of hydrogen atoms also shows very good agreement with protonation states and hydrogen-bond network observed in neutron-diffraction structures. The method is implemented as a computational protocol in Accelrys Discovery Studio and provides a fast and easy way to study the effect of pH on many important mechanisms such as enzyme catalysis, ligand binding, protein-protein interactions, and protein stability.     

Constant pH molecular dynamics with proton tautomerism

The current article describes a new two-dimensional lambda-dynamics method to include proton tautomerism in continuous constant pH molecular dynamics (CPHMD) simulations. The two-dimensional lambda-dynamics framework is used to devise a tautomeric state titration model for the CPHMD simulations involving carboxyl and histidine residues. Combined with the GBSW implicit solvent model, the new method is tested on titration simulations of blocked histidine and aspartic acid as well as two benchmark proteins, turkey ovomucoid third domain (OMTKY3) and ribonuclease A (RNase A). A detailed analysis of the errors inherent to the CPHMD methodology as well as those due to the underlying solvation model is given. The average absolute error for the computed pKa values in OMTKY3 is 1.0 pK unit. In RNase A the average absolute errors for the carboxyl and histidine residues are 1.6 and 0.6 pK units, respectively. In contrast to the previous work, the new model predicts the correct sign for all the pKa shifts, but one, in the benchmark proteins. The predictions of the tautomeric states of His12 and His48 and the conformational states of His48 and His119 are in agreement with experiment. Based on the simulations of OMTKY3 and RNase A, the current work has demonstrated the capability of the CPHMD technique in revealing pH-coupled conformational dynamics of protein side chains.     

Multiple pH regime molecular dynamics simulation for pK calculations

Ionisation equilibria in proteins are influenced by conformational flexibility, which can in principle be accounted for by molecular dynamics simulation. One problem in this method is the bias arising from the fixed protonation state during the simulation. Its effect is mostly exhibited when the ionisation behaviour of the titratable groups is extrapolated to pH regions where the predetermined protonation state of the protein may not be statistically relevant, leading to conformational sampling that is not representative of the true state. In this work we consider a simple approach which can essentially reduce this problem. Three molecular dynamics structure sets are generated, each with a different protonation state of the protein molecule expected to be relevant at three pH regions, and pK calculations from the three sets are combined to predict pK over the entire pH range of interest. This multiple pH molecular dynamics approach was tested on the GCN4 leucine zipper, a protein for which a full data set of experimental data is available. The pK values were predicted with a mean deviation from the experimental data of 0.29 pH units, and with a precision of 0.13 pH units, evaluated on the basis of equivalent sites in the dimeric GCN4 leucine zipper.     

Effects of Triton X-100 on proton transfer and in the photocycle of archaerhodopsin 4

Archaerhodopsin 4 (AR4), a retinal-containing membrane protein, exhibits a reversed order of proton release and uptake at neutral pH, as compared to the well-known bacteriorhodopsin (BR). In a preceding report, we stated that Triton X-100 solubilized the claret membrane containing AR4 (CM) into monomeric proteins and altered the time order in AR4 at neutral pH. The present study examined the mechanism underlying this phenomenon. We employed a photoelectrochemical cell suitable for observation of the proton pumping behaviors of both the membrane patch and detergent-solubilized proteins over a wide pH range. The pK(a) values of the proton release complex (PRC) in the initial state and the M state were determined with this device. The pK(a) of PRC of monomeric AR4 decreased to a value lower than 7.0 in the photocycle, allowing early proton release at neutral pH. The pK(a) of PRC in the initial state was also strongly affected by solubilization.     

Solvent isotope effects and the pH dependence of laccase activity under steady-state conditions

Solvent isotope effects and the pH dependence of laccase catalysis under steady-state conditions were examined with a rapid reductant to assess the potential roles of protein protic groups and the catalytic mechanism. The pH dependence of both reductant-dependent and reductant-independent steps showed bell-shaped profiles implicating at least two protic groups in each case. The apparent pKa values were: for the reductant-independent step(s), pK alpha 1 = 8.98 +/- 0.02 and pK alpha 2 = 5.91 +/- 0.03; for the reductant-dependent step(s), pK' alpha 1 = 7.55 +/- 0.12, pK' alpha 2 = 8.40 +/- 0.23. No solvent isotope effect on reductant-dependent steps was detected other than a standard shift effect. However, a significant solvent isotope effect on a reductant-independent step(s) was observed; kH/kD = 2.12 at the pH optimum of 7.5. The concentration dependence of the D2O effect indicated that a single proton was involved. Simulations of the p(H,D) data suggested that the solvent isotope effect was associated with the protein protic group required in its undissociated form (pK alpha 2). The pH effects on reductant-dependent steps are apparently associated with reductant-dependent steps that occur between O2 binding and water formation in the catalytic reaction sequence.     

Predicting extreme pKa shifts in staphylococcal nuclease mutants with constant pH molecular dynamics

Accurate computational methods of determining protein and nucleic acid pK(a) values are vital to understanding pH-dependent processes in biological systems. In this article, we use the recently developed method constant pH molecular dynamics (CPHMD) to explore the calculation of highly perturbed pK(a) values in variants of staphylococcal nuclease (SNase). Simulations were performed using the replica exchange (REX) protocol for improved conformational sampling with eight temperature windows, and yielded converged proton populations in a total sampling time of 4 ns. Our REX-CPHMD simulations resulted in calculated pK(a) values with an average unsigned error (AUE) of 0.75 pK units for the acidic residues in Δ + PHS, a hyperstable variant of SNase. For highly pK(a)-perturbed SNase mutants with known crystal structures, our calculations yielded an AUE of 1.5 pK units and for those mutants based on modeled structures an AUE of 1.4 pK units was found. Although a systematic underestimate of pK shifts was observed in most of the cases for the highly perturbed pK mutants, correlations between conformational rearrangement and plasticity associated with the mutation and error in pK(a) prediction was not evident in the data. This study further extends the scope of electrostatic environments explored using the REX-CPHMD methodology and suggests that it is a reliable tool for rapidly characterizing ionizable amino acids within proteins even when modeled structures are employed.     

Discrete charge calculations of potentiometric titrations for globular proteins: sperm whale myoglobin, hemoglobin alpha chain, cytochrome c.

The modified Tanford-Kirkwood theory of Shire et al. for intramolecular electrostatic interactions has been applied to hydrogen ion equilibria of sperm whale ferrimyoglobin, human hemoglobin α-chain and horse cytochrome c. The model employs two sets of parameters derived from the crystalline protein structures, first, the atomic coordinates of charged amino acid residues and, second, static accessibility factors to reflect their solvent exposure. In addition, a consistent set of intrinsic pK values (pK_int) for the individual groups is employed. The theoretical pK values at half-titration for individual groups in each protein correspond to the available observed pK values, and the theoretical titration curves compare closely with experimental potentiometric curves.     

pH replica-exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.     

Thermodynamic limitations of photosynthetic water oxidation at high proton concentrations

In oxygenic photosynthesis, solar energy drives the oxidation of water catalyzed by a Mn(4)Ca complex bound to the proteins of Photosystem II. Four protons are released during one turnover of the water oxidation cycle (S-state cycle), implying thermodynamic limitations at low pH. For proton concentrations ranging from 1 nm (pH 9) to 1 mm (pH 3), we have characterized the low-pH limitations using a new experimental approach: a specific pH-jump protocol combined with time-resolved measurement of the delayed chlorophyll fluorescence after nanosecond flash excitation. Effective pK values were determined for low-pH inhibition of the light-induced S-state transitions: pK(1)=3.3 ± 0.3, pK(2)=3.5 ± 0.2, and pK(3)≈pK(4)=4.6 ± 0.2. Alkaline inhibition was not observed. An extension of the classical Kok model facilitated assignment of these four pK values to specific deprotonation steps in the reaction cycle. Our results provide important support to the extended S-state cycle model and criteria needed for assessment of quantum chemical calculations of the mechanism of water oxidation. They also imply that, in intact organisms, the pH in the lumen compartment can hardly drop below 5, thereby limiting the ΔpH contribution to the driving force of ATP synthesis.     

pK(a) values for side-chain carboxyl groups of a PGB1 variant explain salt and pH-dependent stability

Determination of pK(a) values of titrating residues in proteins provides a direct means of studying electrostatic coupling as well as pH-dependent stability. The B1 domain of protein G provides an excellent model system for such investigations. In this work, we analyze the observed pK(a) values of all carboxyl groups in a variant of PGB1 (T2Q, N8D, N37D) at low and high ionic strength as determined using (1)H-(13)C heteronuclear NMR in a structural context. The pK(a) values are used to calculate the pH-dependent stability in low and high salt and to investigate electrostatic coupling in the system. The observed pK(a) values can explain the pH dependence of protein stability but require pK(a) shifts relative to model values in the unfolded state, consistent with persistent residual structure in the denatured state. In particular, we find that most of the deviations from the expected random coil values can be explained by a significantly upshifted pK(a) value. We show also that (13)C backbone carbonyl data can be used to study electrostatic coupling in proteins and provide specific information on hydrogen bonding and electrostatic potential at nontitrating sites.     

pH corrections and protein ionization in water/guanidinium chloride.

More than 30 years ago, Nozaki and Tanford reported that the pK values for several amino acids and simple substances in 6 M guanidinium chloride differed little from the corresponding values in low salt (Nozaki, Y., and C. Tanford. 1967. J. Am. Chem. Soc. 89:736-742). This puzzling and counter-intuitive result hinders attempts to understand and predict the proton uptake/release behavior of proteins in guanidinium chloride solutions, behavior which may determine whether the DeltaG(N-D) values obtained from guanidinium chloride-induced denaturation data can actually be interpreted as the Gibbs energy difference between the native and denatured states (Bolen, D. W., and M. Yang. 2000. Biochemistry. 39:15208-15216). We show in this work that the Nozaki-Tanford result can be traced back to the fact that glass-electrode pH meter readings in water/guanidinium chloride do not equal true pH values. We determine the correction factors required to convert pH meter readings in water/guanidinium chloride into true pH values and show that, when these corrections are applied, the effect of guanidinium chloride on the pK values of simple substances is found to be significant and similar to that of NaCl. The results reported here allow us to propose plausible guanidinium chloride concentration dependencies for the pK values of carboxylic acids in proteins and, on their basis, to reproduce qualitatively the proton uptake/release behavior for the native and denatured states of several proteins (ribonuclease A, alpha-chymotrypsin, staphylococcal nuclease) in guanidinium chloride solutions. Finally, the implications of the pH correction for the experimental characterization of protein folding energetics are briefly discussed.     

Optimizing pKa computation in proteins with pH adapted conformations

pK(A) in proteins are determined by electrostatic energy computations using a small number of optimized protein conformations derived from crystal structures. In these protein conformations hydrogen positions and geometries of salt bridges on the protein surface were determined self-consistently with the protonation pattern at three pHs (low, ambient, and high). Considering salt bridges at protein surfaces is most relevant, since they open at low and high pH. In the absence of these conformational changes, computed pK(A)(comp) of acidic (basic) groups in salt bridges underestimate (overestimate) experimental pK(A)(exp), dramatically. The pK(A)(comp) for 15 different proteins with 185 known pK(A)(exp) yield an RMSD of 1.12, comparable with two other methods. One of these methods is fully empirical with many adjustable parameters. The other is also based on electrostatic energy computations using many non-optimized side chain conformers but employs larger dielectric constants at short distances of charge pairs that diminish their electrostatic interactions. These empirical corrections that account implicitly for additional conformational flexibility were needed to describe the energetics of salt bridges appropriately. This is not needed in the present approach. The RMSD of the present approach improves if one considers only strongly shifted pK(A)(exp) in contrast to the other methods under these conditions. Our method allows interpreting pK(A)(comp) in terms of pH dependent hydrogen bonding pattern and salt bridge geometries. A web service is provided to perform pK(A) computations.     

Negatively charged residues and hydrogen bonds tune the ligand histidine pKa values of Rieske iron-sulfur proteins

Rieske proteins carry a redox-active iron-sulfur cluster, which is bound by two histidine and two cysteine side chains. The reduction potential of Rieske proteins depends on pH. This pH dependence can be described by two pK(a) values, which have been assigned to the two iron-coordinating histidines. Rieske proteins are commonly grouped into two major classes: Rieske proteins from quinol-oxidizing cytochrome bc complexes, in which the ligand histidines titrate in the physiological pH range, and bacterial ferredoxin Rieske proteins, in which the ligand histidines are protonated at physiological pH. In the study presented here, we have calculated pK(a) values of the cluster ligand histidines using a combined density functional theory/continuum electrostatics approach. Experimental pK(a) values for a bc-type and a ferredoxin Rieske protein could be reproduced. We could identify functionally important differences between the two proteins: hydrogen bonds toward the cluster, which are present in bc-type Rieske proteins, and negatively charged residues, which are present in ferredoxin Rieske proteins. We removed these differences by mutating the proteins in our calculations. The Rieske centers in the mutated proteins have very similar pK(a) values. We thus conclude that the studied structural differences are the main reason for the different pH-titration behavior of the proteins. Interestingly, the shift caused by neutralizing the negative charges in ferredoxin Rieske proteins is larger than the shift caused by removing the hydrogen bonds toward the cluster in bc-type Rieske proteins.     

Assessing the quality of the homology-modeled 3D structures from electrostatic standpoint: test on bacterial nucleoside monophosphate kinase families

In this study, we address the issue of performing meaningful pK(a) calculations using homology modeled three-dimensional (3D) structures and analyze the possibility of using the calculated pK(a) values to detect structural defects in the models. For this purpose, the 3D structure of each member of five large protein families of a bacterial nucleoside monophosphate kinases (NMPK) have been modeled by means of homology-based approach. Further, we performed pK(a) calculations for the each model and for the template X-ray structures. Each bacterial NMPK family used in the study comprised on average 100 members providing a pool of sequences and 3D models large enough for reliable statistical analysis. It was shown that pK(a) values of titratable groups, which are highly conserved within a family, tend to be conserved among the models too. We demonstrated that homology modeled structures with sequence identity larger than 35% and gap percentile smaller than 10% can be used for meaningful pK(a) calculations. In addition, it was found that some highly conserved titratable groups either exhibit large pK(a) fluctuations among the models or have pK(a) values shifted by several pH units with respect to the pK(a) calculated for the X-ray structure. We demonstrated that such case usually indicates structural errors associated with the model. Thus, we argue that pK(a) calculations can be used for assessing the quality of the 3D models by monitoring fluctuations of the pK(a) values for highly conserved titratable residues within large sets of homologous proteins.     

Acid-base equilibria in rhodopsin: dependence of the protonation state of glu134 on its environment

Glutamic acid E134 in rhodopsin is part of a highly conserved triad, D(E)RY, located near the cytoplasmic lipid/water interface in transmembrane helix 3 of G protein-coupled receptors (GPCRs). A large body of experimental evidence suggests that the protonation of E134 plays a role in the mechanism of activation of rhodopsin and other GPCRs as well. For E134 to change its protonation state, its pK(a) value must shift from values below physiological pH to higher values. Because of the proximity of the triad to the lipid/water interface, it was hypothesized that a change in solvent around E134 from water to lipid could induce such a shift in pK(a). To test this hypothesis, the pK(a) values of the titratable amino acid residues in rhodopsin have been calculated and the change in solvent around E134 was modeled by shifting the position of the lipid/water interface. The approach used to carry out the pK(a) calculations takes into account the partial immersion of transmembrane proteins in lipid. Qualitative experimental evidence is available for several residues regarding their likely protonation state in rhodopsin at or near physiological pH. Comparison of the calculated pK(a) values with these experimental findings shows good agreement between the two. Notably, glutamic acids E122 and E181 were found to be protonated. The pK(a) values were then calculated for a range of lipid/water interface positions. Although the surrounding solvent of several titratable residues changed from water to lipid in this range, leading to pK(a) shifts in most cases, only for E134 would the shift lead to a change in protonation state at physiological pH. Thus, our results show that the protonation state of E134 is particularly sensitive to its environment. This sensitivity together with the location of E134 near the actual position of the lipid/water interface could be a strategic element in the mechanism of activation of rhodopsin.     

The uptake of protons by heme-linked ionizable groups on azide binding to methemoglobin

When azide ion reacts with methemoglobin in unbuffered solution the pH of the solution increases. This phenomenon is associated with increases in the pK values of heme-linked ionizable groups on the protein which give rise to an uptake of protons from solution. We have determined as a functional of pH the proton uptake, delta h+, on azide binding to methemoglobin at 20 degrees C. Data for methemoglobins A (human), guinea pig and pigeon are fitted to a theoretical expression based on the electrostatic effect of these sets of heme-linked ionizable groups on the binding of the ligand. From these fits the pK values of heme-linked ionizable groups are obtained for liganded and unliganded methemoglobins. In unliganded methemoglobin pK1, which is associated with carboxylic acid groups, ranges between 4.0 and 5.5 for the three methemoglobins; pK2, which is associated with histidines and terminal amino groups, ranges from 6.2 to 6.7. In liganded methemoglobin pK1 lies between 5.8 and 6.3 and pK2 varies from 8.1 to 8.5. The pH dependences of the apparent equilibrium constants for azide binding to the three methemoglobins at 20 degrees C are well accounted for with the pK values calculated from the variation of delta h+ with pH.     

Evidence for the rate of the final step in the bacteriorhodopsin photocycle being controlled by the proton release group: R134H mutant

Light absorbed by bacteriorhodopsin (bR) leads to a proton being released at the extracellular surface of the purple membrane. Structural studies as well as studies of mutants of bR indicate that several groups form a pathway for proton transfer from the Schiff base to the extracellular surface. These groups include D85, R82, E204, E194, and water molecules. Other residues may be important in tuning the initial state pK(a) values of these groups and in mediating light-induced changes of the pK(a) values. A potentially important residue is R134: it is located close to E194 and might interact electrostatically to affect the pK(a) of E194 and light-induced proton release. In this study we investigated effects of the substitution of R134 with a histidine on light-induced proton release and on the photocycle transitions associated with proton transfer. By measuring the light-induced absorption changes versus pH, we found that the R134H mutation results in an increase in the pK(a) of the proton release group in both the M (0.6 pK unit) and O (0.7 pK unit) intermediate states. This indicates the importance of R134 in tuning the pK(a) of the group that, at neutral and high pH, releases the proton upon M formation (fast proton release) and that, at low pH, releases the proton simultaneously with O decay (slow proton release). The higher pK(a) of the proton release group found in R134H correlates with the slowing of the rate of the O --> bR transition at low pH and probably is the cause of this slowing. The pH dependence of the fraction of the O intermediate is altered in R134H compared to the WT but is similar to that in the E194D mutant: a very small amount of O is present at neutral pH, but the fraction of O increases greatly upon decreasing the pH. These results provide further support for the hypothesis that the O --> bR transition is controlled by the rate of deprotonation of the proton release group. These data also provide further evidence for the importance of the R134-E194 interaction in modulating proton release from D85 after light has led to its being protonated.