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

Constant-pH Brownian Dynamics Simulations of a Protein near a Charged Surface

We have developed a rigid-body Brownian dynamics algorithm that allows for simulations of a globular protein suspended in an ionic solution confined by a charged planar boundary, with an explicit treatment of pH-dependent protein protonation equilibria and their couplings to the electrostatic potential of the plane. Electrostatic interactions are described within a framework of the continuum Poisson-Boltzmann model, whereas protein-plane hydrodynamic interactions are evaluated based on analytical expressions for the position- and orientation-dependent near-wall friction tensor of a spheroid. The algorithm was applied to simulate near-surface diffusion of lysozyme in solutions having pH in the range 4-10 and ionic strengths of 10 and 150 mM. As a reference, we performed Brownian dynamics simulations in which the protein is assigned a fixed, most probable protonation state, appropriate for given solution conditions and unaffected by the presence of the charged plane, and Brownian dynamics simulations in which the protein probes possible protonation states with the pH-dependent probability, but these variations are not coupled to the electric field generated by the boundary. We show that electrostatic interactions with the negatively charged plane substantially modify probabilities of different protonation states of lysozyme and shift protonation equilibria of both acidic and basic amino acid side chains toward higher pH values. Consequently, equilibrium energy distributions, equilibrium position-orientation distributions, and functions that characterize rotational dynamics, which for a protein with multiple ionization sites, such as lysozyme, in the presence of a charged obstacle are pH-dependent, are significantly affected by the approach taken to incorporate the solution pH into simulations.

A Hierarchical Approach to Predict Conformation-Dependent Histidine Protonation States in Stable and Flexible Proteins

Solution acidity measured by pH is an important environmental factor that affects protein structure. It influences the protonation state of protein residues, which in turn may be coupled to protein conformational changes, unfolding, and ligand binding. It remains difficult to compute and measure the p K_a of individual residues, as well as to relate them to pH-dependent protein transitions. This paper presents a hierarchical approach to compute the p K_a of individual protonatable residues, specifically histidines, coupled with underlying structural changes of a protein. A fast and efficient free energy perturbation (FEP) algorithm has also been developed utilizing a fast implementation of standard molecular dynamics (MD) algorithms. Specifically, a CUDA version of the AMBER MD engine is used in this paper. Eight histidine p K_a's are computed in a diverse set of pH stable proteins to demonstrate the proposed approach's utility and assess the predictive quality of the AMBER FF99SB force field. A reference molecule is carefully selected and tested for convergence. A hierarchical approach is used to model p K_a's of the six histidine residues of the diphtheria toxin translocation domain (DTT), which exhibits a diverse ensemble of individual conformations and pH-dependent unfolding. The hierarchical approach consists of first sampling equilibrium conformational ensembles of a protein with protonated and neutral histidine residues via long equilibrium MD simulations (Flores-Canales, J. C.; et al. bioRxiv, 2019, 572040). A clustering method is then used to identify sampled protein conformations, and p K_a's of histidines in each protein conformation are computed. Finally, an ensemble averaging formalism is developed to compute weighted average histidine p K_a's. These can be compared with an apparent experimentally measured p K_a of the DTT protein and thus allows us to propose a mechanism of pH-dependent unfolding of the DTT protein.

pK(A) in proteins solving the Poisson-Boltzmann equation with finite elements

Knowledge on pK(A) values is an eminent factor to understand the function of proteins in living systems. We present a novel approach demonstrating that the finite element (FE) method of solving the linearized Poisson-Boltzmann equation (lPBE) can successfully be used to compute pK(A) values in proteins with high accuracy as a possible replacement to finite difference (FD) method. For this purpose, we implemented the software molecular Finite Element Solver (mFES) in the framework of the Karlsberg+ program to compute pK(A) values. This work focuses on a comparison between pK(A) computations obtained with the well-established FD method and with the new developed FE method mFES, solving the lPBE using protein crystal structures without conformational changes. Accurate and coarse model systems are set up with mFES using a similar number of unknowns compared with the FD method. Our FE method delivers results for computations of pK(A) values and interaction energies of titratable groups, which are comparable in accuracy. We introduce different thermodynamic cycles to evaluate pK(A) values and we show for the FE method how different parameters influence the accuracy of computed pK(A) values.

pKa values in proteins determined by electrostatics applied to molecular dynamics trajectories

For a benchmark set of 194 measured pKa values in 13 proteins, electrostatic energy computations are performed in which pKa values are computed by solving the Poisson-Boltzmann equation. In contrast to the previous approach of Karlsberg(+) (KB(+)) that essentially used protein crystal structures with variations in their side chain conformations, the present approach (KB2(+)MD) uses protein conformations from four molecular dynamics (MD) simulations of 10 ns each. These MD simulations are performed with different specific but fixed protonation patterns, selected to sample the conformational space for the different protonation patterns faithfully. The root-mean-square deviation between computed and measured pKa values (pKa RMSD) is shown to be reduced from 1.17 pH units using KB(+) to 0.96 pH units using KB2(+)MD. The pKa RMSD can be further reduced to 0.79 pH units, if each conformation is energy-minimized with a dielectric constant of εmin = 4 prior to calculating the electrostatic energy. The electrostatic energy expressions upon which the computations are based have been reformulated such that they do not involve terms that mix protein and solvent environment contributions and no thermodynamic cycle is needed. As a consequence, conformations of the titratable residues can be treated independently in the protein and solvent environments. In addition, the energy terms used here avoid the so-called intrinsic pKa and can therefore be interpreted without reference to arbitrary protonation states and conformations.

Understanding the pK(a) of redox cysteines: the key role of hydrogen bonding

Many cellular functions involve cysteine chemistry via thiol-disulfide exchange pathways. The nucleophilic cysteines of the enzymes involved are activated as thiolate. A thiolate is much more reactive than a neutral thiol. Therefore, determining and understanding the pK(a)s of functional cysteines are important aspects of biochemistry and molecular biology with direct implications for redox signaling. Here, we describe the experimental and theoretical methods to determine cysteine pK(a) values, and we examine the factors that control these pK(a)s. Drawing largely on experience gained with the thioredoxin superfamily, we examine the roles of solvation, charge-charge, helix macrodipole, and hydrogen bonding interactions as pK(a)-modulating factors. The contributions of these factors in influencing cysteine pK(a)s and the associated chemistry, including the relevance for the reaction kinetics and thermodynamics, are discussed. This analysis highlights the critical role of direct hydrogen bonding to the cysteine sulfur as a key factor modulating the equilibrium between thiol S-H and thiolate S(-). This role is easily understood intuitively and provides a framework for biochemical functional insights.

Influence of the Length of the Alanine Spacer on the Acidic-Basic Properties of the Ac-Lys-(Ala)(n)-Lys-NH(2) Peptides (n = 0, 1, 2, …, 5)

By using the potentiometric titration method, we have determined the pK_a values of the two terminal lysine groups in six alanine-based peptides differing in the length of the alanine chain: Ac–Lys–Lys–NH₂ (KK), Ac–Lys–Ala–Lys–NH₂ (KAK), Ac–Lys–Ala–Ala–Lys–NH₂ (KAK2), Ac–Lys–Ala–Ala–Ala–Lys–NH₂ (KAK3), Ac–Lys–Ala–Ala–Ala–Ala–Lys–NH₂ (KAK4), and Ac–Lys–Ala–Ala–Ala–Ala–Ala–Lys–NH₂ (KAK5) in aqueous solution. For each compound, the model of two stepwise acid–base equilibria was fitted to the potentiometric-titration data. As expected, the pK_a values of the lysine groups increase with increasing length of the alanine spacer, which means that the influence of the electrostatic field between one charged group on the other decreases with increasing length of the alanine spacer. However, for KAK3, the pK_a1 value (8.20) is unusually small and pK_a2 (11.41) is remarkably greater than pK_a1, suggesting that the two groups are close to each other and, in turn, that a chain-reversal conformation is present for this peptide. Starting with KAK3, the differences between pK_a1 and pK_a2 decrease; however, for the longest peptide (KAK5), the values of pK_a1 and pK_a2 still differ by about 1 unit, i.e., by more than the value of log₁₀ (4) = 0.60 that is a limiting value for the pK_a difference of dicarboxylic acids with increasing methylene-spacer length. Consequently, some interactions between the two charged groups are present and, in turn, a bent shape occurs even for the longest of the peptides studied.

Defining the role of salt bridges in protein stability

Although the energetic balance of forces stabilizing proteins has been established qualitatively over the last decades, quantification of the energetic contribution of particular interactions still poses serious problems. The reasons are the strong cooperativity and the interdependence ofnoncovalent interactions. Salt bridges are a typical example. One expects that ionizable side chains frequently form ion pairs in innumerable crystal structures. Since electrostatic attraction between opposite charges is strong per se, salt bridges can intuitively be regarded as an important factor stabilizing the native structure. Is that really so? In this chapter we critically reassess the available methods to delineate the role ofelectrostatic interactions and salt bridges to protein stability, and discuss the progress and the obstacles in this endeavor. The basic problem is that formation of salt bridges depends on the ionization properties of the participating groups, which is significantly influenced by the protein environment. Furthermore, salt bridges experience thermal fluctuations, continuously break and re-form, and their lifespan in solution is governed by the flexibility of the protein. Finally, electrostatic interactions are long-range and might be significant in the unfolded state, thus seriously influencing the energetic profile. Elimination of salt bridges by protonation/deprotonation at extreme pH or by mutation provides only rough energetic estimates, since there is no way to account for the nonadditive response of the protein moiety. From what we know so far, the strength of electrostatic interactions is strongly context-dependent, yet it is unlikely that salt bridges are dominant factors governing protein stability. Nevertheless, proteins from thermophiles and hyperthermophiles exhibit more, and frequently networked, salt bridges than proteins from the mesophilic counterparts. Increasing the thermal (not the thermodynamic) stability of proteins by optimization of charge-charge interactions is a good example for an evolutionary solution utilizing physical factors.

Acidic-basic properties of three alanine-based peptides containing acidic and basic side chains: comparison between theory and experiment

The purpose of this work was to evaluate the effect of the nature of the ionizable end groups, and the solvent, on their acid-base properties in alanine-based peptides. Hence, the acid-base properties of three alanine-based peptides: Ac-KK-(A)(7)-KK-NH(2) (KAK), Ac-OO-(A)(7)-DD-NH(2) (OAD), Ac-KK-(A)(7)-EE-NH(2) (KAE), where A, D, E, K, and O denote alanine, aspartic acid, glutamic acid, lysine, and ornithine, respectively, were determined in water and in methanol by potentiometry. With the availability of these data, the ability of two theoretical methods to simulate pH-metric titration of those peptides was assessed: (i) the electrostatically driven Monte Carlo method with the ECEPP/3 force field and the Poisson-Boltzmann approach to compute solvation energy (EDMC/PB/pH), and (ii) the molecular dynamics method with the AMBER force field and the Generalized Born model (MD/GB/pH). For OAD and KAE, pK(a1) and pK(a2) correspond to the acidic side chains. For all three compounds in both solvents, the pK(a1) value is remarkably lower than the pK(a) of a compound modeling the respective isolated side chain, which can be explained by the influence of the electrostatic field from positively charged ornithine or lysine side chains. The experimental titration curves are reproduced well by the MD/GB/pH approach, the agreement being better if restraints derived from NMR measurements are incorporated in the conformational search. Poorer agreement is achieved by the EDMC/PB/pH method.

Structure-based pKa calculations using continuum electrostatics methods

Electrostatic free energy is useful for correlating structure with function in proteins in which ionizable groups play essential functional roles. To this end, the pK(a) values of ionizable groups must be known and their molecular determinants must be understood. Structure-based calculations of electrostatic energies and pK(a) values are necessary for this purpose. This unit describes protocols for pK(a) calculations with continuum electrostatics methods based on the numerical solution of the linearized Poisson-Boltzmann equation by the method of finite differences. Critical discussion of key parameters, approximations, and shortcomings of these methods is included. Two protocols are described for calculations with methods modified empirically to maximize agreement between measured and calculated pK(a) values. Applied judiciously, these methods can contribute useful and novel insight into properties of surface ionizable groups in proteins.

Molecular docking and 3D-QSAR studies of 2-substituted 1-indanone derivatives as acetylcholinesterase inhibitors

AIM

To explore the binding mode of 2-substituted 1-indanone derivatives with acetylcholinesterase (AChE) and provide hints for the future design of new derivatives with higher potency and specificity.

METHODS

The GOLD-docking conformations of the compounds in the active site of the enzyme were used in subsequent studies. The highly reliable and predictive three-dimensional quantitative structure-activity relationship (3D-QSAR) models were achieved by comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA) methods. The predictive capabilities of the models were validated by an external test set. Moreover, the stabilities of the 3D-QSAR models were verified by the leave-4-out cross-validation method.

RESULTS

The CoMFA and CoMSIA models were constructed successfully with a good cross-validated coefficient (q(2)) and a non-cross-validated coefficient (r(2)). The q(2)and r(2)obtained from the leave-1-out cross validation method were 0.784 and 0.974 in the CoMFA model and 0.736 and 0.947 in the CoMSIA model, respectively. The coefficient isocontour maps obtained from these models were compatible with the geometrical and physicochemical properties of AChE.

CONCLUSION

The contour map demonstrated that the binding affinity could be enhanced when the small protonated nitrogen moiety was replaced by a more hydrophobic and bulky group with a highly partial positive charge. The present study provides a better understanding of the interaction between the inhibitors and AChE, which is helpful for the discovery of new compounds with more potency and selective activity.

Calculations of pH-dependent binding of proteins to biological membranes

Binding of proteins to membranes is often accompanied by titration of ionizable residues and is, therefore, dependent on pH. We present a theoretical treatment and computational approach for predicting absolute, pH-dependent membrane binding free energies. The standard free energy of binding, DeltaG, is defined as -RTln(P(b)/P(f)), where P(b) and P(f) are the amounts of bound and free protein. The apparent pK(a) of binding is the pH value at which P(b) and P(f) are equal. Proteins bind to the membrane in the pH range where DeltaG is negative. The components of the binding free energy are (a) the free energy cost of ionization state changes (DeltaG(ion)), (b) the effective energy of transfer from solvent to the membrane surface, (c) the translational/rotational entropy cost of binding, and (d) an ideal entropy term that depends on the relative volume of the bound and free state and therefore depends on lipid concentration. Calculation of the first term requires determination of pK(a) values in solvent and on the membrane surface. All energies required by the method are obtained from molecular dynamics trajectories on an implicit membrane (IMM1-GC). The method is tested on pentalysine and the helical peptide VEEKS, derived from the membrane-binding domain of phosphocholine cytidylyltransferase. The agreement between the measured and the calculated free energies of binding of pentalysine is good. The extent of membrane binding of VEEKS is, however, underestimated compared to experiment. Calculations of the interaction energy between two VEEKS helices on the membrane suggest that the discrepancy is mainly due to the neglect of protein-protein interactions on the membrane surface.

Linear-scaling molecular orbital calculations for the pKa values of ionizable residues in proteins

In this report, we present a computational methodology for the pKa prediction of proteins, based on linear-scaling molecular orbital calculations for their solution-conformations obtained from NMR measurements. The method is used to predict the pKa values of five carboxylic acids (Asp7, Glu10, Glu19, Asp27, and Glu43) in turkey ovomucoid third domain (OMTKY3), and six aspartates residues (Asp 22, Asp 44, Asp 54, Asp 75, Asp 83, and Asp 93) in barnase. For OMTKY3, all the predicted pKa values are within 1 pH units from the available experimental ones, except for the case of Glu 43. For barnase, the root-mean-square deviation from experiment is 1.46 pH units. As a result, the proposed pKa calculation method correctly reproduces the relative order of the pKa values among the carboxylic acids located in different sites of the proteins. The calculated pKa values are decomposed into the contributions of short- and long-range structural difference effects. The results indicate that in both proteins the pKa value of the given carboxylic acid is partially influenced by long-range interactions with distant charged residues, which significantly contribute to determining the relative order of the pKa values. The current methodology based on LSMO provides us useful information about the titration behavior in a protein.

Effect of HPr phosphorylation on structure, dynamics, and interactions in the course of transcriptional control

The serine46-phosphorylated form of the bacterial protein HPr fulfils an essential function in carbon catabolite repression (CCR). Using molecular dynamics (MD) we studied the effect of Ser46 phosphorylation on the molecular properties of HPr and its capability to act as the co-repressor of carbon catabolite protein A (CcpA). The calculated pK (a) values for a representative set of HPr(Ser46P) structures indicate that the phosphate group of HPr(Ser46P) exists predominantly in the unprotonated form under neutral conditions. A hydrogen bond detected in HPr(Ser46P) between one phosphate-group oxygen and a side-chain hydrogen of Asn43-an amino acid conserved in all HPr proteins of Gram-positive bacteria that regulate their carbon consumption by CCR-might fulfil an important role in CcpA-HPr(Ser46P) complex formation. MD simulations show that the Ser46P-Asn43 hydrogen bond present in the unbound structure is replaced by intermolecular interactions upon complex formation. The degree to which amino acids in the CcpA-HPr(Ser46P) interface contribute to cofactor binding was analyzed by in silico alanine scanning. Lys307, Arg303, Asp296, Val300, and Tyr295 of CcpA were identified as important amino acids for the CcpA-HPr(Ser46P) interaction. Three of these residues are directly involved in sensing the correct phosphorylation state at His15(HPr) and Ser46(HPr). A substitution of interface residues Val319, Val314, Ser316, Leu321 and Gln320 by alanine showed that these amino acids, which contact helix alpha2 of HPr(Ser46P), play a less prominent role for complex formation.

Factors influencing the energetics of electron and proton transfers in proteins. What can be learned from calculations

A protein structure should provide the information needed to understand its observed properties. Significant progress has been made in developing accurate calculations of acid/base and oxidation/reduction reactions in proteins. Current methods and their strengths and weaknesses are discussed. The distribution and calculated ionization states in a survey of proteins is described, showing that a significant minority of acidic and basic residues are buried in the protein and that most of these remain ionized. The electrochemistry of heme and quinones are considered. Proton transfers in bacteriorhodopsin and coupled electron and proton transfers in photosynthetic reaction centers, 5-coordinate heme binding proteins and cytochrome c oxidase are highlighted as systems where calculations have provided insight into the reaction mechanism.

Analysing the pH-dependent properties of proteins using pKa calculations

The results of protein pKa calculations are routinely being analysed to understand the pH-dependence of protein characteristics such as stability and catalysis. Systems of functionally important titratable groups are identified from protein from pKa calculations, but the rationalisation of the behaviour of such systems is inherently problematic due to a lack of theoretical tools and methods. I present a number of novel methods for analysing the results of protein pKa calculations which have been embedded in a graphical user interface (pKaTool). In the present paper I present novel methods for assessing the reliability of protein pKa calculations and for analysing the roles of individual residues in determining active site pKa values and the pH-dependence of protein stability. The methods presented are freely available to academic researchers at http://enzyme.ucd.ie/Science/pKa/pKaTool .

Assessing the acid–base and conformational properties of histidine residues in human prion protein (125–228) by means of pK a calculations and molecular dynamics simulations

A thorough study of the acid-base behavior of the four histidines and the other titratable residues of the structured domain of human prion protein (125-228) is presented. By using multi-tautomer electrostatic calculations, average titration curves have been built for all titratable residues, using the whole bundles of NMR structures determined at pH 4.5 and 7.0. According to our results, (1) only histidine residues are likely to be involved in the first steps of the pH-driven conformational transition of prion protein; (2) the pK(a)'s of His140 and His177 are approximately 7.0, whereas those of His155 and His187 are < 5.5. 10-ns long molecular dynamics simulations have been performed on five different models, corresponding to the most significant combinations of histidine protonation states. A critical comparison between the available NMR structures and our computational results (1) confirms that His155 and His187 are the residues whose protonation is involved in the conformational rearrangement of huPrP in mildly acidic condition, and (2) shows how their protonation leads to the destructuration of the C-terminal part of HB and to the loss of the last turn of HA that represent the crucial microscopic steps of the rearrangement.

Assessment of two theoretical methods to estimate potentiometric titration curves of peptides: comparison with experiment

We compared the ability of two theoretical methods of pH-dependent conformational calculations to reproduce experimental potentiometric titration curves of two models of peptides: Ac-K5-NHMe in 95% methanol (MeOH)/5% water mixture and Ac-XX(A)7OO-NH2 (XAO) (where X is diaminobutyric acid, A is alanine, and O is ornithine) in water, methanol (MeOH), and dimethyl sulfoxide (DMSO), respectively. The titration curve of the former was taken from the literature, and the curve of the latter was determined in this work. The first theoretical method involves a conformational search using the electrostatically driven Monte Carlo (EDMC) method with a low-cost energy function (ECEPP/3 plus the SRFOPT surface-solvation model, assumming that all titratable groups are uncharged) and subsequent reevaluation of the free energy at a given pH with the Poisson-Boltzmann equation, considering variable protonation states. In the second procedure, molecular dynamics (MD) simulations are run with the AMBER force field and the generalized Born model of electrostatic solvation, and the protonation states are sampled during constant-pH MD runs. In all three solvents, the first pKa of XAO is strongly downshifted compared to the value for the reference compounds (ethylamine and propylamine, respectively); the water and methanol curves have one, and the DMSO curve has two jumps characteristic of remarkable differences in the dissociation constants of acidic groups. The predicted titration curves of Ac-K5-NHMe are in good agreement with the experimental ones; better agreement is achieved with the MD-based method. The titration curves of XAO in methanol and DMSO, calculated using the MD-based approach, trace the shape of the experimental curves, reproducing the pH jump, while those calculated with the EDMC-based approach and the titration curve in water calculated using the MD-based approach have smooth shapes characteristic of the titration of weak multifunctional acids with small differences between the dissociation constants. Nevertheless, quantitative agreement between theoretically predicted and experimental titration curves is not achieved in all three solvents even with the MD-based approach, which is manifested by a smaller pH range of the calculated titration curves with respect to the experimental curves. The poorer agreement obtained for water than for the nonaqueous solvents suggests a significant role of specific solvation in water, which cannot be accounted for by the mean-field solvation models.

Constant pH molecular dynamics in generalized Born implicit solvent

A new method is proposed for constant pH molecular dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). pK(a) predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to experimental values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quantitatively analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.

Constant-pH molecular dynamics study of protonation-structure relationship in a heptapeptide derived from ovomucoid third domain

Molecular dynamics (MD) simulations with implicit solvent and variable protonation states for titratable residues at constant pH are performed for a short peptide derived from ovomucoid third domain (OMTKY3), acetyl-Ser-Asp-Asn-Lys-Thr-Tyr-Gly-amide (residues 26-32 of OMTKY3). Nuclear magnetic resonance (NMR) measurements indicate that the pK(a) for Asp is 3.6. However, if the charge on Lys is neutralized by acetylation, then the pK(a) for Asp is 4.0. These pK(a)'s, and therefore the Asp-Lys interaction, are insensitive to changes in ionic strength. The constant-pH MD simulations for both variants of the heptapeptide yield Asp pK(a) values that are 0.6-0.9 pH units greater than experimental values, but the difference between the variants that is observed in the NMR experiments is reproduced much better. Moreover, the simulations suggest that Asp-Lys interactions do not dominate the behavior of this heptapeptide, even for normal Lys residue where there is a possibility of forming a salt bridge between negatively charged Asp and positively charged Lys. This is consistent with the experimentally observed independence of Asp pK(a) values with respect to ionic strength. Another important result of the simulations with variable protonation states is that they lead to ensembles of the heptapeptide structures that are different from those derived from simulations with fixed protonation states. It should be stressed that these results are for structures generated entirely by computer simulations without any restrictions imposed by experimental data.

Calculating pKa values in enzyme active sites

The ionization properties of the active-site residues in enzymes are of considerable interest in the study of the catalytic mechanisms of enzymes. Knowledge of these ionization constants (pKa values) often allows the researcher to identify the proton donor and the catalytic nucleophile in the reaction mechanism of the enzyme. Estimates of protein residue pKa values can be obtained by applying pKa calculation algorithms to protein X-ray structures. We show that pKa values accurate enough for identifying the proton donor in an enzyme active site can be calculated by considering in detail only the active-site residues and their immediate electrostatic interaction partners, thus allowing for a large decrease in calculation time. More specifically we omit the calculation of site-site interaction energies, and the calculation of desolvation and background interaction energies for a large number of pairs of titratable groups. The method presented here is well suited to be applied on a genomic scale, and can be implemented in most pKa calculation algorithms to give significant reductions in calculation time with little or no impact on the accuracy of the results. The work presented here has implications for the understanding of enzymes in general and for the design of novel biocatalysts.

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.

Langevin dynamics of proteins at constant pH

An application of the Langevin dynamics algorithm for simulation of protein conformational equilibria at constant pH is presented. The algorithm is used to compute average protonation of titratable groups in ovomucoid third domain, as functions of pH, resulting in data, basically equivalent to the pH dependencies of chemical shifts obtained from multidimensional nuclear magnetic resonance (NMR) spectroscopy, for the protein titratable residues. The pK(a) values obtained from the simulation are in reasonable agreement with experimental data. Possible improvements of this methodology, using achievements from other fields of mesoscopic biomolecular simulations, are also discussed.

Calculation of protein ionization equilibria with conformational sampling: pK(a) of a model leucine zipper, GCN4 and barnase

The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pK(a) calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pK(a) estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve-averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pK(a) calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pK(a) predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pK(a) predictions for the model leucine zipper was reduced to as low as approximately 0.25 pK(a) units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pK(a) calculations, are discussed.

Structure, dynamics and electrostatics of the active site of glutaredoxin 3 from Escherichia coli: comparison with functionally related proteins

The chemistry of active-site cysteine residues is central to the activity of thiol-disulfide oxidoreductases of the thioredoxin superfamily. In these reactions, a nucleophilic thiolate is required, but the associated pK(a) values differ vastly in the superfamily, from less than 4 in DsbA to greater than 7 in Trx. The factors that stabilize this thiolate are, however, not clearly established. The glutaredoxins (Grxs), which are members of this superfamily, contain a Cys-Pro-Tyr-Cys motif in their active site. In reduced Grxs, the pK(a) of the N-terminal active-site nucleophilic cysteine residue is lowered significantly, and the stabilization of the corresponding thiolate is expected to influence the redox potential of these enzymes. Here, we use a combination of long molecular dynamics (MD) simulations, pK(a) calculations, and experimental investigations to derive the structure and dynamics of the reduced active site from Escherichia coli Grx3, and investigate the factors that stabilize the thiolate. Several different MD simulations converged toward a consensus conformation for the active-site cysteine residues (Cys11 and Cys14), after a number of local conformational changes. Key features of the model were tested experimentally by measurement of NMR scalar coupling constants, and determination of pK(a) values of selected residues. The pK(a) values of the Grx3 active-site residues were calculated during the MD simulations, and support the underlying structural model. The structure of Grx3, in combination with the pK(a) calculations, indicate that the pK(a) of the N-terminal active-site cysteine residue in Grx3 is intermediate between that of its counterpart in DsbA and Trx. The pK(a) values in best agreement with experiment are obtained with a low (<4) protein dielectric constant. The calculated pK(a) values fluctuate significantly in response to protein dynamics, which underscores the importance of the details of the underlying structures when calculating pK(a) values. The thiolate of Cys11 is stabilized primarily by direct hydrogen bonding with the amide protons of Tyr13 and Cys14 and the thiol proton of Cys14, rather than by long-range interactions from charged groups or from a helix macrodipole. From the comparison of reduced Grx3 with other members of the thioredoxin superfamily, a unifying theme for the structural basis of thiol pK(a) differences in this superfamily begins to emerge.

Side-chain ionization states in a potassium channel

KcsA is a bacterial K+ channel that is gated by pH. Continuum dielectric calculations on the crystal structure of the channel protein embedded in a low dielectric slab suggest that side chains E71 and D80 of each subunit, which lie adjacent to the selectivity filter region of the channel, form a proton-sharing pair in which E71 is neutral (protonated) and D80 is negatively charged at pH 7. When K+ ions are introduced into the system at their crystallographic positions the pattern of proton sharing is altered. The largest perturbation is for a K+ ion at site S3, i.e., interacting with the carbonyls of T75 and V76. The presence of multiple K+ ions in the filter increases the probability of E71 being ionized and of D80 remaining neutral (i.e., protonated). The ionization states of the protein side chains influence the potential energy profile experienced by a K+ ion as it is translated along the pore axis. In particular, the ionization state of the E71-D80 proton-sharing pair modulates the shape of the potential profile in the vicinity of the selectivity filter. Such reciprocal effects of ion occupancy on side-chain ionization states, and of side-chain ionization states on ion potential energy profiles will complicate molecular dynamics simulations and related studies designed to calculate ion permeation energetics.

Molecular dynamics simulations of protein-tyrosine phosphatase 1B. II. substrate-enzyme interactions and dynamics

Molecular dynamics simulations of protein tyrosine phosphatase 1B (PTP1B) complexed with the phosphorylated peptide substrate DADEpYL and the free substrate have been conducted to investigate 1) the physical forces involved in substrate-protein interactions, 2) the importance of enzyme and substrate flexibility for binding, 3) the electrostatic properties of the enzyme, and 4) the contribution from solvation. The simulations were performed for 1 ns, using explicit water molecules. The last 700 ps of the trajectories was used for analysis determining enthalpic and entropic contributions to substrate binding. Based on essential dynamics analysis of the PTP1B/DADEpYL trajectory, it is shown that internal motions in the binding pocket occur in a subspace of only a few degrees of freedom. In particular, relatively large flexibilities are observed along several eigenvectors in the segments: Arg(24)-Ser(28), Pro(38)-Arg(47), and Glu(115)-Gly(117). These motions are correlated to the C- and N-terminal motions of the substrate. Relatively small fluctuations are observed in the region of the consensus active site motif (H/V)CX(5)R(S/T) and in the region of the WPD loop, which contains the general acid for catalysis. Analysis of the individual enzyme-substrate interaction energies revealed that mainly electrostatic forces contribute to binding. Indeed, calculation of the electrostatic field of the enzyme reveals that only the field surrounding the binding pocket is positive, while the remaining protein surface is characterized by a predominantly negative electrostatic field. This positive electrostatic field attracts negatively charged substrates and could explain the experimentally observed preference of PTP1B for negatively charged substrates like the DADEpYL peptide.