Electrostatic polarization is crucial for reproducing pKa shifts of carboxylic residues in Turkey ovomucoid third domain

AMOEBA Force Field Trajectories Improve Predictions of Accurate pKa Values of the GFP Fluorophore: The Importance of Polarizability and Water Interactions

Precisely quantifying the magnitude, direction, and biological functions of electric fields in proteins has long been an outstanding challenge in the field. The most widely implemented experimental method to measure such electric fields at a particular residue in a protein has been through changes in pKa of titratable residues. While many computational strategies exist to predict these values, it has been difficult to do this accurately or connect predicted results to key structural or mechanistic features of the molecule. Here, we used experimentally determined pKa values of the fluorophore in superfolder green fluorescent protein (GFP) with amino acid mutations made at position Thr 203 to evaluate the pKa prediction ability of molecular dynamics (MD) simulations using a polarizable force field, AMOEBA. Structure ensembles from AMOEBA were used to calculate pKa values of the GFP fluorophore. The calculated pKa values were then compared to trajectories using a conventional fixed charge force field (Amber03 ff). We found that the position of water molecules included in the pKa calculation had opposite effects on the pKa values between the trajectories from AMOEBA and Amber03 force fields. In AMOEBA trajectories, the inclusion of water molecules within 35 Å of the fluorophore decreased the difference between the predicted and experimental values, resulting in calculated pKa values that were within an average of 0.8 pKa unit from the experimental results. On the other hand, in Amber03 trajectories, including water molecules that were more than 5 Å from the fluorophore increased the differences between the calculated and experimental pKa values. The inaccuracy of pKa predictions determined from Amber03 trajectories was caused by a significant stabilization of the deprotonated chromophore's free energy compared to the result in AMOEBA. We rationalize the cutoffs for explicit water molecules when calculating pKa to better predict the electrostatic environment surrounding the fluorophore buried in GFP. We discuss how the results from this work will assist the prospective prediction of pKa values or other electrostatic effects in a wide variety of folded proteins.

A Polarizable Cationic Dummy Metal Ion Model

A novel locally polarizable multisite model based on the original cation dummy atom (CDA) model is described for molecular dynamics simulations of ions in condensed phases. Polarization effects are introduced by the electronegativity equalization model (EEM) method where charges on the metal ion and its dummy atoms can fluctuate to respond to the environment. This model includes explicit polarization and ion-induced interactions and can be coupled with nonpolarizable or polarizable water models, making it more transferable to simpler force fields. This approach allows us to enhance the original fixed charge CDA model where the charge distribution cannot adapt to the local solvent structure. To illustrate the new CDApol model, we examined properties of the Zn2+, Al3+, and Zr4+ ions in aqueous solution. The polarizable model and Lennard-Jones parameters were refined for octahedrally coordinated Zn2+, Al3+, and Zr4+ CDAs to reproduce thermodynamic and geometrical properties. Using this locally polarizable model, we were able to obtain the experimental hydration free energy, ion-oxygen distance, and coordination number coupled with the standard 12-6 Lennard-Jones model. This model can be used in myriad additional applications where local polarization and charge transfer is important.

Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Realistic modeling of biomolecular systems requires an accurate treatment of electrostatics, including electronic polarization. Due to recent advances in physical models, simulation algorithms, and computing hardware, biomolecular simulations with advanced force fields at biologically relevant timescales are becoming increasingly promising. These advancements have not only led to new biophysical insights but also afforded opportunities to advance our understanding of fundamental intermolecular forces. This article describes the recent advances and applications, as well as future directions, of polarizable force fields in biomolecular simulations.

Using polarizable POSSIM force field and fuzzy-border continuum solvent model to calculate pK(a) shifts of protein residues

Our Fuzzy-Border (FB) continuum solvent model has been extended and modified to produce hydration parameters for small molecules using POlarizable Simulations Second-order Interaction Model (POSSIM) framework with an average error of 0.136 kcal/mol. It was then used to compute pKa shifts for carboxylic and basic residues of the turkey ovomucoid third domain (OMTKY3) protein. The average unsigned errors in the acid and base pKa values were 0.37 and 0.4 pH units, respectively, versus 0.58 and 0.7 pH units as calculated with a previous version of polarizable protein force field and Poisson Boltzmann continuum solvent. This POSSIM/FB result is produced with explicit refitting of the hydration parameters to the pKa values of the carboxylic and basic residues of the OMTKY3 protein; thus, the values of the acidity constants can be viewed as additional fitting target data. In addition to calculating pKa shifts for the OMTKY3 residues, we have studied aspartic acid residues of Rnase Sa. This was done without any further refitting of the parameters and agreement with the experimental pKa values is within an average unsigned error of 0.65 pH units. This result included the Asp79 residue that is buried and thus has a high experimental pKa value of 7.37 units. Thus, the presented model is capable or reproducing pKa results for residues in an environment that is significantly different from the solvated protein surface used in the fitting. Therefore, the POSSIM force field and the FB continuum solvent parameters have been demonstrated to be sufficiently robust and transferable. © 2016 Wiley Periodicals, Inc.

Developing multisite empirical force field models for Pt(II) and cisplatin

We have developed empirical force field parameters for Pt(II) and cisplatin. Two force field frameworks were used-modified OPLS-AA and our second-order polarizable POSSIM. A seven-site model was used for the Pt(II) ion. The goal was to create transferable parameter sets compatible with the force field models for proteins and general organic compounds. A number of properties of the Pt(II) ion and its coordination compounds have been considered, including geometries and energies of the complexes, hydration free energy, and radial distribution functions in water. Comparison has been made with experimental and quantum mechanical results. We have demonstrated that both versions are generally capable of reproducing key properties of the system, but the second-order polarizable POSSIM formalism permits more accurate quantitative results to be obtained. For example, the energy of formation of cisplatin as calculated with the modified OPLS-AA exhibited an error of 9.9%, while the POSSIM error for the same quantity was 6.2%. The produced parameter sets are transferable and suitable to be used in protein-metal binding simulations in which position or even coordination of the ion does not have to be constrained using preexisting knowledge. © 2016 Wiley Periodicals, Inc.

Evaluating Parametrization Protocols for Hydration Free Energy Calculations with the AMOEBA Polarizable Force Field

Hydration free energy (HFE) calculations are often used to assess the performance of biomolecular force fields and the quality of assigned parameters. The AMOEBA polarizable force field moves beyond traditional pairwise additive models of electrostatics and may be expected to improve upon predictions of thermodynamic quantities such as HFEs over and above fixed-point-charge models. The recent SAMPL4 challenge evaluated the AMOEBA polarizable force field in this regard but showed substantially worse results than those using the fixed-point-charge GAFF model. Starting with a set of automatically generated AMOEBA parameters for the SAMPL4 data set, we evaluate the cumulative effects of a series of incremental improvements in parametrization protocol, including both solute and solvent model changes. Ultimately, the optimized AMOEBA parameters give a set of results that are not statistically significantly different from those of GAFF in terms of signed and unsigned error metrics. This allows us to propose a number of guidelines for new molecule parameter derivation with AMOEBA, which we expect to have benefits for a range of biomolecular simulation applications such as protein-ligand binding studies.

Direct computation of parameters for accurate polarizable force fields

We present an improved electronic linear response model to incorporate polarization and charge-transfer effects in polarizable force fields. This model is a generalization of the Atom-Condensed Kohn-Sham Density Functional Theory (DFT), approximated to second order (ACKS2): it can now be defined with any underlying variational theory (next to KS-DFT) and it can include atomic multipoles and off-center basis functions. Parameters in this model are computed efficiently as expectation values of an electronic wavefunction, obviating the need for their calibration, regularization, and manual tuning. In the limit of a complete density and potential basis set in the ACKS2 model, the linear response properties of the underlying theory for a given molecular geometry are reproduced exactly. A numerical validation with a test set of 110 molecules shows that very accurate models can already be obtained with fluctuating charges and dipoles. These features greatly facilitate the development of polarizable force fields.

POSSIM: Parameterizing Complete Second-Order Polarizable Force Field for Proteins

Previously, we reported development of a fast polarizable force field and software named POSSIM (POlarizable Simulations with Second order Interaction Model). The second-order approximation permits the speed up of the polarizable component of the calculations by ca. an order of magnitude. We have now expanded the POSSIM framework to include a complete polarizable force field for proteins. Most of the parameter fitting was done to high-level quantum mechanical data. Conformational geometries and energies for dipeptides have been reproduced within average errors of ca. 0.5 kcal/mol for energies of the conformers (for the electrostatically neutral residues) and 9.7° for key dihedral angles. We have also validated this force field by running Monte Carlo simulations of collagen-like proteins in water. The resulting geometries were within 0.94 Å root-mean-square deviation (RMSD) from the experimental data. We have performed additional validation by studying conformational properties of three oligopeptides relevant in the context of N-glycoprotein secondary structure. These systems have been previously studied with combined experimental and computational methods, and both POSSIM and benchmark OPLS-AA simulations that we carried out produced geometries within ca. 0.9 Å RMSD of the literature structures. Thus, the performance of POSSIM in reproducing the structures is comparable with that of the widely used OPLS-AA force field. Furthermore, our fitting of the force field parameters for peptides and proteins has been streamlined compared with the previous generation of the complete polarizable force field and relied more on transferability of parameters for nonbonded interactions (including the electrostatic component). The resulting deviations from the quantum mechanical data are similar to those achieved with the previous generation; thus, the technique is robust, and the parameters are transferable. At the same time, the number of parameters used in this work was noticeably smaller than that of the previous generation of our complete polarizable force field for proteins; thus, the transferability of this set can be expected to be greater, and the danger of force field fitting artifacts is lower. Therefore, we believe that this force field can be successfully applied in a wide variety of applications to proteins and protein–ligand complexes.

Polarizable simulations with second order interaction model (POSSIM) force field: developing parameters for protein side-chain analogues

A previously introduced polarizable simulations with second-order interaction model (POSSIM) force field has been extended to include parameters for small molecules serving as models for peptide and protein side-chains. Parameters have been fitted to permit reproducing many-body energies, gas-phase dimerization energies, and geometries and liquid-phase heats of vaporization and densities. Quantum mechanical and experimental data have been used as the target for the fitting. The POSSIM framework combines accuracy of a polarizable force field and computational efficiency of the second-order approximation of the full-scale induced point dipole polarization formalism. The resulting parameters can be used for simulations of the parameterized molecules themselves or their analogues. In addition to this, these force field parameters are currently being used in further development of the POSSIM fast polarizable force field for proteins.

Polarizable force fields

This chapter provides an overview of the most common methods for including an explicit description of electronic polarization in molecular mechanics force fields: the induced point dipole, shell, and fluctuating charge models. The importance of including polarization effects in biomolecular simulations is discussed, and some of the most important achievements in the development of polarizable biomolecular force fields to date are highlighted.

Calculating pKa values for substituted phenols and hydration energies for other compounds with the first-order Fuzzy-Border continuum solvation model

We have computed pK(a) values for 11 substituted phenol compounds using the continuum Fuzzy-Border (FB) solvation model. Hydration energies for 40 other compounds, including alkanes, alkenes, alkynes, ketones, amines, alcohols, ethers, aromatics, amides, heterocycles, thiols, sulfides, and acids have been calculated. The overall average unsigned error in the calculated acidity constant values was equal to 0.41 pH units and the average error in the solvation energies was 0.076 kcal/mol. We have also reproduced pK(a) values of propanoic and butanoic acids within about 0.1 pH units from the experimental values by fitting the solvation parameters for carboxylate ion carbon and oxygen atoms. The FB model combines two distinguishing features. First, it limits the amount of noise which is common in numerical treatment of continuum solvation models by using fixed-position grid points. Second, it uses either second- or first-order approximation for the solvent polarization, depending on a particular implementation. These approximations are similar to those used for solute and explicit solvent fast polarization treatment which we developed previously. This article describes results of using the first-order technique. This approximation places the presented methodology between the Generalized Born and Poisson-Boltzmann continuum solvation models with respect to their accuracy of reproducing the many-body effects in modeling a continuum solvent.

Effects of lysine substitution on stability of polyalanine alpha-helix

We have studied stability of polyalanine alpha-helices with lysine residues added at C-and N-termini in gas-phase and aqueous solution. Monte Carlo simulations with the fixed-charges OPLS-AA and our polarizable POSSIM force fields were carried out. The results of the simulations confirm previously observed phenomena of the helix being stable with the LYS residue on the C-terminus and losing its helical structure if the charged LYS residue is located at the N-terminus of the polypeptide in gas-hase. Both OPLS-AA and POSSIM force fields performed essentially similarly, thus validity of the both for reproducing and predicting structures of such polypeptides has been confirmed. We have also studied the effect of replacing the normal N- and C-termini with methyl capping (this approach is often used in computational studies). Our results have demonstrated that the structure and stability of the polypeptides do not depend significantly on such a substitution although details of the resulting structure may change. The liquid-state simulations produced stable alpha-helixes regardless of the position of the protonated lysine residue. Overall, we have validated our polarizable POSSIM force field and the techniques used in the simulations, since the change of the helix structure as a function of the position of the LYS residue depends on a fine balance of energy contributions, and our methodology reproduced this balance well.

Importance of electrostatic polarizability in calculating cysteine acidity constants and copper(I) binding energy of Bacillus subtilis CopZ

CopZ is a copper chaperone from Bacillus subtilis. It is an important part of Cu(I) trafficking. We have calculated pK(a) values for the CXXC motif of this protein, which is responsible for the Cu(I) binding, and the Cu(I) binding constants. Polarizable and fixed-charges formalisms were used, and solvation parameters for the both models have been refitted. We had to partially redevelop parameters for the protonated and deprotonated cysteine residues. We have discovered that the polarizable force field (PFF) is qualitatively superior and allows a uniformly better level of energetic results. The PFF pK(a) values for cysteine are within about 0.8-2.8 pH units of the experimental data, while the fixed-charges OPLS formalism yields errors of up to tens of units. The PFF magnitude of the copper binding energy is about 10 kcal/mol or 50% higher than the experimental value, while the using the refitted OPLS parameters leads to an overall positive binding energy, thus predicting no thermodynamically stable complex. At the same time, the agreement of the polarizable S···Cu(I) distances with the experimental results is within 0.08 Å range, and the nonpolarizable calculations lead to an error of about 0.4 Å. Moreover, the accuracy of the PFF has been achieved without any explicit fitting to either pK(a) or CopZ···Cu(I) binding energies. We believe that this makes our polarizable technique a choice method in reproducing protein-copper binding and further supports the notion that explicit treatment of electrostatic polarization is crucial in many biologically relevant studies, especially ion binding and transport.

Differential geometry based solvation model II: Lagrangian formulation

Solvation is an elementary process in nature and is of paramount importance to more sophisticated chemical, biological and biomolecular processes. The understanding of solvation is an essential prerequisite for the quantitative description and analysis of biomolecular systems. This work presents a Lagrangian formulation of our differential geometry based solvation models. The Lagrangian representation of biomolecular surfaces has a few utilities/advantages. First, it provides an essential basis for biomolecular visualization, surface electrostatic potential map and visual perception of biomolecules. Additionally, it is consistent with the conventional setting of implicit solvent theories and thus, many existing theoretical algorithms and computational software packages can be directly employed. Finally, the Lagrangian representation does not need to resort to artificially enlarged van der Waals radii as often required by the Eulerian representation in solvation analysis. The main goal of the present work is to analyze the connection, similarity and difference between the Eulerian and Lagrangian formalisms of the solvation model. Such analysis is important to the understanding of the differential geometry based solvation model. The present model extends the scaled particle theory of nonpolar solvation model with a solvent-solute interaction potential. The nonpolar solvation model is completed with a Poisson-Boltzmann (PB) theory based polar solvation model. The differential geometry theory of surfaces is employed to provide a natural description of solvent-solute interfaces. The optimization of the total free energy functional, which encompasses the polar and nonpolar contributions, leads to coupled potential driven geometric flow and PB equations. Due to the development of singularities and nonsmooth manifolds in the Lagrangian representation, the resulting potential-driven geometric flow equation is embedded into the Eulerian representation for the purpose of computation, thanks to the equivalence of the Laplace-Beltrami operator in the two representations. The coupled partial differential equations (PDEs) are solved with an iterative procedure to reach a steady state, which delivers desired solvent-solute interface and electrostatic potential for problems of interest. These quantities are utilized to evaluate the solvation free energies and protein-protein binding affinities. A number of computational methods and algorithms are described for the interconversion of Lagrangian and Eulerian representations, and for the solution of the coupled PDE system. The proposed approaches have been extensively validated. We also verify that the mean curvature flow indeed gives rise to the minimal molecular surface and the proposed variational procedure indeed offers minimal total free energy. Solvation analysis and applications are considered for a set of 17 small compounds and a set of 23 proteins. The salt effect on protein-protein binding affinity is investigated with two protein complexes by using the present model. Numerical results are compared to the experimental measurements and to those obtained by using other theoretical methods in the literature.

Electrostatic polarization is crucial in reproducing Cu(I) interaction energies and hydration

We have explored the suitability of fixed-charges and polarizable force fields for modeling interactions of the monovalent Cu(I) ion. Parameters for this ion have been tested and refitted within the fixed-charges OPLS-AA and polarizable force field (PFF) frameworks. While this ion plays an important role in many protein interactions, the attention to it in developing empirical force fields is limited. Our PFF parameters for the copper ion worked very well for the Cu(I) interactions with water, while both the original OPLS2005 and our refitted OPLS versions moderately underestimated the copper-water interaction energy. However, the greatest problem in using the nonpolarizable fixed-charges OPLS force field was observed while calculating interaction energies and distances for Cu(I)-benzene complexes. The OPLS2005 model underestimates the interaction energy by a factor of 4. Refitting the OPLS parameters reduced this underestimation to a factor of 2.2-2.4, but only at a cost of distorting the complex geometry. At the same time, the polarizable calculations had an error of about 4%. Moreover, we then used the PFF and nonpolarizable refitted OPLS models for finding free energy of hydration for copper ion via molecular dynamics simulations. While the OPLS calculations lead to a 22% error in the solvation energy, the PFF result was off by only 1.8%. This was achieved with no refitting of the parameters but simply by employing the model developed for the Cu(I) interaction with a single water molecule. We believe that the presented results not only lead to a conclusion about a qualitatively greater suitability of polarizable force fields for simulating molecular interactions with ions but also attest to the excellent level of transferability of PFF parameters.

Polarizable Simulations with Second order Interaction Model (POSSIM) force field: Developing parameters for alanine peptides and protein backbone

A previously introduced POSSIM (POlarizable Simulations with Second order Interaction Model) force field has been extended to include parameters for alanine peptides and protein backbones. New features were introduced into the fitting protocol, as compared to the previous generation of the polarizable force field for proteins. A reduced amount of quantum mechanical data was employed in fitting the electrostatic parameters. Transferability of the electrostatics between our recently developed NMA model and the protein backbone was confirmed. Binding energy and geometry for complexes of alanine dipeptide with a water molecule were estimated and found in a good agreement with high-level quantum mechanical results (for example, the intermolecular distances agreeing within ca. 0.06Å). Following the previously devised procedure, we calculated average errors in alanine di- and tetra-peptide conformational energies and backbone angles and found the agreement to be adequate (for example, the alanine tetrapeptide extended-globular conformational energy gap was calculated to be 3.09 kcal/mol quantim mechanically and 3.14 kcal/mol with the POSSIM force field). However, we have now also included simulation of a simple alpha-helix in both gas-phase and water as the ultimate test of the backbone conformational behavior. The resulting alanine and protein backbone force field is currently being employed in further development of the POSSIM fast polarizable force field for proteins.

Differential geometry based solvation model I: Eulerian formulation

This paper presents a differential geometry based model for the analysis and computation of the equilibrium property of solvation. Differential geometry theory of surfaces is utilized to define and construct smooth interfaces with good stability and differentiability for use in characterizing the solvent-solute boundaries and in generating continuous dielectric functions across the computational domain. A total free energy functional is constructed to couple polar and nonpolar contributions to the salvation process. Geometric measure theory is employed to rigorously convert a Lagrangian formulation of the surface energy into an Eulerian formulation so as to bring all energy terms into an equal footing. By minimizing the total free energy functional, we derive coupled generalized Poisson-Boltzmann equation (GPBE) and generalized geometric flow equation (GGFE) for the electrostatic potential and the construction of realistic solvent-solute boundaries, respectively. By solving the coupled GPBE and GGFE, we obtain the electrostatic potential, the solvent-solute boundary profile, and the smooth dielectric function, and thereby improve the accuracy and stability of implicit solvation calculations. We also design efficient second order numerical schemes for the solution of the GPBE and GGFE. Matrix resulted from the discretization of the GPBE is accelerated with appropriate preconditioners. An alternative direct implicit (ADI) scheme is designed to improve the stability of solving the GGFE. Two iterative approaches are designed to solve the coupled system of nonlinear partial differential equations. Extensive numerical experiments are designed to validate the present theoretical model, test computational methods, and optimize numerical algorithms. Example solvation analysis of both small compounds and proteins are carried out to further demonstrate the accuracy, stability, efficiency and robustness of the present new model and numerical approaches. Comparison is given to both experimental and theoretical results in the literature.

Molecular acidity: A quantitative conceptual density functional theory description

Accurate predictions of molecular acidity using ab initio and density functional approaches are still a daunting task. Using electronic and reactivity properties, one can quantitatively estimate pKa values of acids. In a recent paper [S. B. Liu and L. G. Pedersen, J. Phys. Chem. A 113, 3648 (2009)], we employed the molecular electrostatic potential (MEP) on the nucleus and the sum of valence natural atomic orbital (NAO) energies for the purpose. In this work, we reformulate these relationships on the basis of conceptual density functional theory and compare the results with those from the thermodynamic cycle method. We show that MEP and NAO properties of the dissociating proton of an acid should satisfy the same relationships with experimental pKa data. We employ 27 main groups and first to third row transition metal-water complexes as illustrative examples to numerically verify the validity of these strong linear correlations. Results also show that the accuracy of our approach and that of the conventional method through the thermodynamic cycle are statistically similar.

Polarizable Simulations with Second order Interaction Model - force field and software for fast polarizable calculations: Parameters for small model systems and free energy calculations

We are presenting POSSIM (POlarizable Simulations with Second order Interaction Model) - a software package and a set of parameters designed for molecular simulations. The key feature of POSSIM is that the electrostatic polarization is taken into account using a previously introduced fast formalism. This permits cutting computational cost of using the explicit polarization by about an order of magnitude. In this article, parameters for water, methane, ethane, propane, butane, methanol and NMA are introduced. These molecules are viewed as model systems for protein simulations. We have achieved our goal of ca. 0.5 kcal/mol accuracy for gas-phase dimerization energies and no more than 2% deviations in liquid state heats of vaporization and densities. Moreover, free energies of hydration of the polarizable methane, ethane and methanol have been calculated using the statistical perturbation theory. These calculations are a model for calculating protein pKa shifts and ligand binding affinities. The free energies of hydration were found to be 2.12 kcal/mol, 1.80 kcal/mol and -4.95 kcal/mol for methane, ethane and methanol, respectively. The experimentally determined literature values are 1.91 kcal/mol, 1.83 kcal/mol and -5.11 kcal/mol. The POSSIM average error in these absolute free energies of hydration is only about 0.13 kcal/mol. Using the statistical perturbation theory with polarizable force fields is not widespread, and we believe that this work opens road to further development of the POSSIM force field and its applications for obtaining accurate energies in protein-related computer modeling.

Reproducing basic pKa values for turkey ovomucoid third domain using a polarizable force field

We have extended our previous studies of calculating acidity constants for the acidic residues found in the turkey ovomucoid third domain protein (OMTKY3) by determining the relative pKa values for the basic residues (Lys13, Arg21, Lys29, Lys34, His52, and Lys55). A polarizable force field (PFF) was employed. The values of the pKa were found by direct comparison of energies of solvated protonated and deprotonated forms of the protein. Poisson-Boltzmann (PBF) and surface generalized Born (SGB) continuum solvation models represent the hydration, and a nonpolarizable fixed-charge OPLS-AA force field was used for comparison. Our results indicate that (i) the pKa values of the basic residues can be found in close agreement with the experimental values when a PFF is used in conjunction with the PBF solvation model, (ii) it is sufficient to take into the account only the residues which are in close proximity (hydrogen bonded) to the residue in question, and (iii) the PBF solvation model is superior to the SGB solvation model for these pKa calculations. The average error with the PBF/PFF model is only 0.7 pH unit, compared with 2.2 and 6.1 units for the PBF/OPLS and SGB/OPLS, respectively. The maximum deviation of the PBF/PFF results from the experimental values is 1.7 pH units compared with 6.0 pH units for the PBF/OPLS. Moreover, the best results were obtained while using an advanced nonpolar energy calculation scheme. The overall conclusion is that this methodology and force field are suitable for the accurate assessment of pKa shifts for both acidic and basic protein residues.

Estimation of molecular acidity via electrostatic potential at the nucleus and valence natural atomic orbitals

An effective approach of estimating molecular pK(a) values from simple density functional calculations is proposed in this work. Both the molecular electrostatic potential (MEP) at the nucleus of the acidic atom and the sum of valence natural atomic orbitals are employed for three categories of compounds, amines and anilines, carbonyl acids and alcohols, and sulfonic acids and thiols. A strong correlation between experimental pK(a) values and each of these two quantities for each of the three categories has been discovered. Moreover, if the MEP is subtracted by the isolated atomic MEP for each category of compounds, we observe a single unique linear relationship between the resultant MEP difference and experimental pK(a) data of amines, anilines, carbonyl acids, alcohols, sulfonic acids, thiols, and their substituents. These results can generally be utilized to simultaneously estimate pK(a) values at multiple sites with a single calculation for either relatively small molecules in drug design or amino acids in proteins and macromolecules.

Assessing the performance of implicit solvation models at a nucleic acid surface

Implicit solvation models are popular alternatives to explicit solvent methods due to their ability to "pre-average" solvent behavior and thus reduce the need for computationally-expensive sampling. Previously, we have demonstrated that Poisson-Boltzmann models for polar solvation and integral-based models for nonpolar solvation can reproduce explicit solvation forces in a low-charge density protein system. In the present work, we examine the ability of these continuum models to describe solvation forces at the surface of a RNA hairpin. While these models do not completely describe all of the details of solvent behavior at this highly-charged biomolecular interface, they do provide a reasonable description of average solvation forces and therefore show significant promise for developing more robust implicit descriptions of solvent around nucleic acid systems for use in biomolecular simulation and modeling. Additionally, we observe fairly good transferability in the nonpolar model parameters optimized for protein systems, suggesting its robustness for modeling general nonpolar solvation phenomena in biomolecular systems.

pKa of residue 66 in Staphylococal nuclease. I. Insights from QM/MM simulations with conventional sampling

A combined quantum mechanical/molecular mechanical (QM/MM) potential function is used in a thermodynamic integration approach to calculate the pK(a) of residue 66 in two mutants (V66E, V66D) of Staphylococal nuclease relative to solution. Despite the similarity in chemical nature and experimentally measured pK(a) of the two buried titritable residues, the behaviors of the two mutants and the computed pK(a) values vary greatly in the simulations. For Glu66, the side chain is consistently observed to spontaneously flip out from the protein interior during titration, and the overall protein structure remains stable throughout the simulations. The computed pK(a) shifts using conventional sampling techniques with multiple nanoseconds per lambda window (Set A and B) are generally close to the experimental value, therefore indicating that large-scale conformational rearrangements are not as important for V66E as suggested by the recent study of Warshel and co-worker. For Asp66, by contrast, flipping of the shorter side chain is not sufficient for getting adequate solvent stabilization of the ionized state. As a result, more complex behaviors such as partial unfolding of a nearby beta-sheet region is observed, and the computed pK(a) shift is substantially higher than the experimental value unless Asp66 is biased to adopt the similar configurations as Glu66 in the V66E simulations. Collectively, these studies suggest that the lack of electronic polarization is not expected to be the dominant source of error in microscopic pK(a) shift calculations, while the need of enhanced sampling is more compelling for predicting the pK(a) of buried residues. Furthermore, the comparison between V66E and V66D also highlights that the microscopic interpretation of similar apparent pK(a) values and effective "dielectric constants" of proteins can vary greatly in terms of the residues that make key contributions and the scale of structural/hydration response to titration, the latter of which is difficult to predict a priori. Perturbative analyses of interactions that contribute to the titration free energy point to mutants that can be used to verify the microscopic mechanisms of titration in V66E/D SNase proteins.