Protein Structure Prediction with a Combined Solvation Free Energy-Molecular Mechanics Force Field

Abstract

Models of protein structure are frequently used to determine the physical characteristics of a protein when the crystal structure is not available. We developed a procedure to optimize such models, by use of a combined solvation free energy and molecular mechanics force field. Appropriately chosen atomic solvation parameters were defined using the criterion that the resulting protein model should deviate least from the crystal structure upon a forty picosecond molecular dynamics simulation carried out using the combined force field. Several tests were performed to refine the set of atomic solvation parameters which best complement the molecular mechanics forces. Four sets of parameters from the literature were tested and an empirically optimized set is proposed. The parameters are defined on a well characterized small molecule (alanyl dipeptide) and on the highly refined crystal structure of rat trypsin, and then tested on a second highly refined crystal structure of α-lytic protease. The new set of atomic solvation parameters provides a significant improvement over molecular mechanics alone in energy minimization of protein structures. This combined force field also has advantages over the use of explicit solvent as it is possible to take solvent effects into account during energetic conformational searching when modeling a homologous protein structure from a known crystal structure.     

An Effective Solvation Term Based on Atomic Occupancies for Use in Protein Simulations

Abstract

Several approaches to the treatment of solvent effects based on continuum models are reviewed and a new method based on occupied atomic volumes (occupancies) is proposed and tested. The new method describes protein-water interactions in terms of atomic solvation parameters, which represent the solvation free energy per unit of volume. These parameters were determined for six different atoms types, using experimental free energies of solvation. The method was implemented in the GROMOS and PRESTO molecular simulation program suites. Simulations with the solvation term require 20-50% more CPU time than the corresponding vacuum simulations and are approximately 20 times faster than explicit water simulations. The method and parameters were tested by carrying out 200 ps simulations of BPTI in water, in vacuo, and with the solvation term. The performance of the solvation term was assessed by comparing the structures and energies from the solvation simulations with the equivalent quantities derived from several BPTI crystal structures and from the explicit water and vacuum simulations. The model structures were evaluated in terms of exposed total surface, buried and exposed polar surfaces, secondary structure preservation, number of hydrogen bonds, energy contributions, and positional deviations from BPTI crystal structures. Vacuum simulations produced unrealistic structures with respect to all criteria applied. The structures resulting from the simulations with explicit water were closer to the 5PTI crystal structure, although part of the secondary structure dissolved. The simulations with the effective solvation term produce structures that are normal according to all evaluations and in most respects are remarkably similar to the 5PTI crystal structure despite considerable positional fluctuations during the simulations. The segments where the model and crystal structures differ are known to be flexible and the observed difference may be physically realistic. The effective solvation term based on occupancies is not only very efficient in terms of computer time but also results in meaningful structural properties for BPTI. It may therefore be generally useful in molecular dynamics of macromolecules.     

Parameterization and Application of an Implicit Solvent Model for Macromolecules

Abstract

As the field of theoretical biophysics begins to recognize systems of longer timescales and larger magnitude, rapid approaches for investigating these systems are required. One promising simplification of the typical system of a solute surrounded by water is the use of implicit solvation models. The generalized Born implicit solvent offers a rapid approach for computing the electrostatic effects of bulk solvent without the explicit representation of water molecules. This report describes the parameterization of a generalized Born (GB) model for protein and nucleic acid structures. As a demonstration of the usefulness of this approach, the GB model is applied toward the discrimination of misfolded and properly folded protein structures. This study attempts to illustrate the potential of the GB model for molecular dynamics simulations over longer timescales as well as for screening large structural databases.     

Role of Hydrogen Bonds in Protein-DNA Recognition: A Comparison of Generalized Born and Finite Difference Poisson-Boltzmann Solvation Treatments

Abstract

Hydrogen bonds have been accredited with a major role historically, in the formation and stabilization of biomolecular structures. The formation of hydrogen bonds at protein-DNA interfaces in aqueous medium involves not only favorable interactions of the donor and acceptor functional groups but also a loss of interactions between these groups with the solvent water. We have investigated the energetics of about 500 potential hydrogen bonds occuring at protein-DNA interfaces incorporating some recent improvements in biomolecular force fields and solvation treatments. We present here results of our assessment of hydrogen bond contributions to the overall standard free energy of formation of protein-DNA complexes obtained with the generalized Born model and finite difference Poisson- Boltzmann methodology for solvation in conjunction with AMBER force field. Our results support the emerging view on the role of electrostatics in general and that of hydrogen bonds in particular which is that hydrogen bonds do not drive protein-DNA complex formation by virtue of the unfavourable cost of the electrostatics of desolvation. They however, act to stabilize the complex once it is formed.     

Free energies of protein decoys provide insight into determinants of protein stability

We have calculated the stability of decoy structures of several proteins (from the CASP3 models and the Park and Levitt decoy set) relative to the native structures. The calculations were performed with the force field-consistent ES/IS method, in which an implicit solvent (IS) model is used to calculate the average solvation free energy for snapshots from explicit simulations (ESs). The conformational free energy is obtained by adding the internal energy of the solute from the ESs and an entropic term estimated from the covariance positional fluctuation matrix. The set of atomic Born radii and the cavity-surface free energy coefficient used in the implicit model has been optimized to be consistent with the all-atom force field used in the ESs (cedar/gromos with simple point charge (SPC) water model). The decoys are found to have a consistently higher free energy than that of the native structure; the gap between the native structure and the best decoy varies between 10 and 15 kcal/mole, on the order of the free energy difference that typically separates the native state of a protein from the unfolded state. The correlation between the free energy and the extent to which the decoy structures differ from the native (as root mean square deviation) is very weak; hence, the free energy is not an accurate measure for ranking the structurally most native-like structures from among a set of models. Analysis of the energy components shows that stability is attained as a result of three major driving forces: (1) minimum size of the protein-water surface interface; (2) minimum total electrostatic energy, which includes solvent polarization; and (3) minimum protein packing energy. The detailed fit required to optimize the last term may underlie difficulties encountered in recovering the native fold from an approximate decoy or model structure.     

A Practical Implicit Membrane Potential for NMR Structure Calculations of Membrane Proteins

The highly anisotropic environment of the lipid bilayer membrane imposes significant constraints on the structures and functions of membrane proteins. However, NMR structure calculations typically use a simple repulsive potential that neglects the effects of solvation and electrostatics, because explicit atomic representation of the solvent and lipid molecules is computationally expensive and impractical for routine NMR-restrained calculations that start from completely extended polypeptide templates. Here, we describe the extension of a previously described implicit solvation potential, eefxPot, to include a membrane model for NMR-restrained calculations of membrane protein structures in XPLOR-NIH. The key components of eefxPot are an energy term for solvation free energy that works together with other nonbonded energy functions, a dedicated force field for conformational and nonbonded protein interaction parameters, and a membrane function that modulates the solvation free energy and dielectric screening as a function of the atomic distance from the membrane center, relative to the membrane thickness. Initial results obtained for membrane proteins with structures determined experimentally in lipid bilayer membranes show that eefxPot affords significant improvements in structural quality, accuracy, and precision. Calculations with eefxPot are straightforward to implement and can be used to both fold and refine structures, as well as to run unrestrained molecular-dynamics simulations. The potential is entirely compatible with the full range of experimental restraints measured by various techniques. Overall, it provides a useful and practical way to calculate membrane protein structures in a physically realistic environment.     

Solvation free energies of nucleic acid bases in ionic liquids

The solvation free energies of five nucleic acid bases in [C_nbim]Br (where n = 2, 4, 6) ionic liquids (ILs) were computed using the Bennett acceptance ratio (BAR) method employing molecular dynamics simulations. The computed free energies using BAR were in agreement with other methods. The large and negative predicted free energies of the bases in ILs indicated that the bases were better solvated in the ILs rather than in water. Hydrogen bonding interactions between polar sites of the bases and ILs’ ions significantly contributed to the solvation mechanism.     

Conformational sampling and polarization of Asp26 in pKa calculations of thioredoxin

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

Parametric sensitivity analysis of avian pancreatic polypeptide (APP)

Computer simulations utilizing a classical force field have been widely used to study biomolecular properties. It is important to identify the key force field parameters or structural groups controlling the molecular properties. In the present paper the sensitivity analysis method is applied to study how various partial charges and solvation parameters affect the equilibrium structure and free energy of avian pancreatic polypeptide (APP). The general shape of APP is characterized by its three principal moments of inertia. A molecular dynamics simulation of APP was carried out with the OPLS/Amber force field and a continuum model of solvation energy. The analysis pinpoints the parameters which have the largest (or smallest) impact on the protein equilibrium structure (i.e., the moments of inertia) or free energy. A display of the protein with its atoms colored according to their sensitivities illustrates the patterns of the interactions responsible for the protein stability. The results suggest that the electrostatic interactions play a more dominant role in protein stability than the part of the solvation effect modeled by the atomic solvation parameters. © 1995 Wiley-Liss, Inc.     

Solvation Free Energy of the Peptide Group: Its Model Dependence and Implications for the Additive-Transfer Free-Energy Model of Protein Stability

The group-additive decomposition of the unfolding free energy of a protein in an osmolyte solution relative to that in water poses a fundamental paradox: whereas the decomposition describes the experimental results rather well, theory suggests that a group-additive decomposition of free energies is, in general, not valid. In a step toward resolving this paradox, here we study the peptide-group transfer free energy. We calculate the vacuum-to-solvent (solvation) free energies of (Gly)_n and cyclic diglycine (cGG) and analyze the data according to experimental protocol. The solvation free energies of (Gly)_n are linear in n, suggesting group additivity. However, the slope interpreted as the free energy of a peptide unit differs from that for cGG scaled by a factor of half, emphasizing the context dependence of solvation. However, the water-to-osmolyte transfer free energies of the peptide unit are relatively independent of the peptide model, as observed experimentally. To understand these observations, a way to assess the contribution to the solvation free energy of solvent-mediated correlation between distinct groups is developed. We show that linearity of solvation free energy with n is a consequence of uniformity of the correlation contributions, with apparent group-additive behavior in the water-to-osmolyte transfer arising due to their cancellation. Implications for inferring molecular mechanisms of solvent effects on protein stability on the basis of the group-additive transfer model are suggested.     

Solvation Free Energy of the Peptide Group: Its Model Dependence and Implications for the Additive-Transfer Free-Energy Model of Protein Stability

The group-additive decomposition of the unfolding free energy of a protein in an osmolyte solution relative to that in water poses a fundamental paradox: whereas the decomposition describes the experimental results rather well, theory suggests that a group-additive decomposition of free energies is, in general, not valid. In a step toward resolving this paradox, here we study the peptide-group transfer free energy. We calculate the vacuum-to-solvent (solvation) free energies of (Gly)_n and cyclic diglycine (cGG) and analyze the data according to experimental protocol. The solvation free energies of (Gly)_n are linear in n, suggesting group additivity. However, the slope interpreted as the free energy of a peptide unit differs from that for cGG scaled by a factor of half, emphasizing the context dependence of solvation. However, the water-to-osmolyte transfer free energies of the peptide unit are relatively independent of the peptide model, as observed experimentally. To understand these observations, a way to assess the contribution to the solvation free energy of solvent-mediated correlation between distinct groups is developed. We show that linearity of solvation free energy with n is a consequence of uniformity of the correlation contributions, with apparent group-additive behavior in the water-to-osmolyte transfer arising due to their cancellation. Implications for inferring molecular mechanisms of solvent effects on protein stability on the basis of the group-additive transfer model are suggested.     

Molecular dynamics and quantum chemical studies of solvent effects on cyclo glycylglycine and glycylalanine dipeptides

Hydrogen bond (H-bond) interactions between the two cyclo dipeptides, cyclo(glycyl-glycine) (CGG) and cyclo(glycyl-alanine) (CGA), and water have been studied using molecular dynamics (MD) and quantum chemical methods. The MD studies have been carried out on CGG and CGA in water using fixed charge force field (AMBER ff03) for over 10 ns with a MD time step of 2 fs. The results of this study show that the solvation pattern influences the conformations of the cyclo dipeptides. Following molecular simulations, post Hartree–Fock and density functional theory methods have been used to explore the molecular properties of the cyclo dipeptides in gaseous and aqueous phase environments. The self-consistent reaction field theory has been used to optimise the cyclopeptides in diethyl ether (? = 4.3) and water (? = 78.5), and the solvent effects have been analysed. A cluster of eight water molecules leads to the formation of first solvation shell of CGG and CGA and the strong H-bonding mainly contributes to the interaction energies. The H-bond interactions have been analysed by the calculation of electron density ρ(r) and its Laplacian ▽²ρ(r) at bond critical points using atoms in molecules theory. The natural bond orbital analysis was carried out to reveal the nature of H-bond interactions. In the solvated complexes, the keto carbons registered the maximum NMR chemical shifts.     

Implicit solvent models

Implicit solvent models for biomolecular simulations are reviewed and their underlying statistical mechanical basis is discussed. The fundamental quantity that implicit models seek to approximate is the solute potential of mean force, which determines the statistical weight of solute conformations, and which is obtained by averaging over the solvent degrees of freedom. It is possible to express the total free energy as the reversible work performed in two successive steps. First, the solute is inserted in the solvent with zero atomic partial charges; second, the atomic partial charges of the solute are switched from zero to their full values. Consequently, the total solvation free energy corresponds to a sum of non-polar and electrostatic contributions. These two contributions are often approximated by simple geometrical models (such as solvent exposed area models) and by macroscopic continuum electrostatics, respectively. One powerful route is to approximate the average solvent density distribution around the solute, i.e. the solute-solvent density correlation functions, as in statistical mechanical integral equations. Recent progress with semi-analytical approximations makes continuum electrostatics treatments very efficient. Still more efficient are fully empirical, knowledge-based models, whose relation to explicit solvent treatments is not fully resolved, however. Continuum models that treat both solute and solvent as dielectric continua are also discussed, and the relation between the solute fluctuations and its macroscopic dielectric constant(s) clarified.     

A Simple Solvation Model Along with A Multibody Dynamics Strategy (MBO(N)D) Produces Stable DNA Simulations that are Faster than Traditional Atomistic Methods

Abstract

We are developing solvation strategies that complement the speed advantage of MBO(N)D (a multibody simulation approach developed by Moldyn) for simulating biomolecular systems. In this report we propose to approximate the effect of bulk waters on DNA by using only a thin layer of waters proximate to the surface of DNA (which we will call the ‘thin shell approach’ or TSA). We will show that the TSA combined with substructuring (the grouping of atoms into rigid or flexible bodies) of the Dickerson dodecamer produces good comparisons with standard atomistic methods (over a nanosecond trajectory) as judged by a variety of DNA specific geometric (e.g., CURVES output) and dynamics (power spectra) properties. The MBO(N)D method, however, was faster than atomistic by a factor of six using the same solvation strategy and factor of 70 when compared to fully solvated atomistic system. The key to the speed of MBO(N)D is in its ability to use large time steps during dynamics. By keeping only a shell of molecules of water proximate to the dodecamer, we limit artifacts due to surface tension at the water-vacuum interface. These proximate waters are fairly immobile as compared to those in bulk and therefore do not severely limit the time step in the simulation. The strengths and limitations of this solvation approach, and future directions, will also be discussed.     

Description of hydration free energy density as a function of molecular physical properties

A method to calculate the solvation free energy density (SFED) at any point in the cavity surface or solvent volume surrounding a solute is proposed. In the special case in which the solvent is water, the SFED is referred to as the hydration free energy density (HFED). The HFED is described as a function of some physical properties of the molecules. These properties are represented by simple basis functions. The hydration free energy of a solute was obtained by integrating the HFED over the solvent volume surrounding the solute, using a grid model. Of 34 basis functions that were introduced to describe the HFED, only six contribute significantly to the HFED. These functions are representations of the surface area and volume of the solute, of the polarization and dispersion of the solute, and of two types of electrostatic interactions between the solute and its environment. The HFED is described as a linear combination of these basis functions, evaluated by summing the interaction energy between each atom of the solute with a grid point in the solvent, where each grid point is a representation of a finite volume of the solvent. The linear combination coefficients were determined by minimizing the error between the calculated and experimental hydration free energies of 81 neutral organic molecules that have a variety of functional groups. The calculated hydration free energies agree well with the experimental results. The hydration free energy of any other solute molecule can then be calculated by summing the product of the linear combination coefficients and the basis functions for the solute.     

Salting-in of neopentane in the aqueous solutions of urea and glycine-betaine

The effects of urea and glycine-betaine (GB) osmolytes on the hydrophobic interactions of neopentane in water have been studied using molecular dynamics simulations. From the study of the potentials of mean force, it is observed that both urea and GB decrease the association and solvation of neopentane. The calculated equilibrium constants show that urea and GB decrease the population of solvent-separated minima of neopentane. The hydrophobic association as well as solvation of neopentane molecules are stabilised by entropy and enthalpy in the mixtures. The radial distribution functions (RDFs) and running coordination numbers of water, urea and GB molecules show that neopentane shows salting-in behaviour in aqueous-GB, aqueous-urea and aqueous-urea-GB mixtures. Neopentane is preferentially solvated by GB in aqueous-GB and preferentially solvated by urea in aqueous-urea-GB solutions. The preferential solvation of neopentane by GB suggests that GB decreases the interaction between neopentane molecules i.e. salting-in of neopentane. The calculated solvation free energies and radial density profiles of neopentane also support the salting-in behaviour of neopentane in the mixtures of these osmolytes.     

Intrinsic pKas of ionizable residues in proteins: An explicit solvent calculation for lysozyme

Molecular dynamics simulations of triclinic hen egg white lysozyme in aqueous solution were performed to calculate the intrinsic pK_as of 14 ionizable residues. An all-atom model was used for both solvent and solute, and a single 180 ps simulation in conjunction with a Gaussian fluctuation analysis method was used. An advantage of the Gaussian fluctuation method is that it only requires a single simulation of the system in a reference state to calculate all the pK_as in the protein, in contrast to multiple simulations for the free energy perturbation method. pK_int shifts with respect to reference titratable residues were evaluated and compared to results obtained using a finite difference Poisson-Boltzmann (FDPB) method with a continuum solvent model; overall agreement with the direction of the shifts was generally observed, though the magnitude of the shifts was typically larger with the explicit solvent model. The contribution of the first solvation shell to the total charging free energies of the titratable groups was explicitly evaluated and found to be significant. Dielectric shielding between pairs of titratable groups was examined and found to be smaller than expected. The effect of the approximations used to treat the long-range interactions on the pK_int shifts is discussed. © 1994 Wiley-Liss, Inc.