Free energy landscape of protein folding in water: explicit vs. implicit solvent

The Generalized Born (GB) continuum solvent model is arguably the most widely used implicit solvent model in protein folding and protein structure prediction simulations; however, it still remains an open question on how well the model behaves in these large-scale simulations. The current study uses the beta-hairpin from C-terminus of protein G as an example to explore the folding free energy landscape with various GB models, and the results are compared to the explicit solvent simulations and experiments. All free energy landscapes are obtained from extensive conformation space sampling with a highly parallel replica exchange method. Because solvation model parameters are strongly coupled with force fields, five different force field/solvation model combinations are examined and compared in this study, namely the explicit solvent model: OPLSAA/SPC model, and the implicit solvent models: OPLSAA/SGB (Surface GB), AMBER94/GBSA (GB with Solvent Accessible Surface Area), AMBER96/GBSA, and AMBER99/GBSA. Surprisingly, we find that the free energy landscapes from implicit solvent models are quite different from that of the explicit solvent model. Except for AMBER96/GBSA, all other implicit solvent models find the lowest free energy state not the native state. All implicit solvent models show erroneous salt-bridge effects between charged residues, particularly in OPLSAA/SGB model, where the overly strong salt-bridge effect results in an overweighting of a non-native structure with one hydrophobic residue F52 expelled from the hydrophobic core in order to make better salt bridges. On the other hand, both AMBER94/GBSA and AMBER99/GBSA models turn the beta-hairpin in to an alpha-helix, and the alpha-helical content is much higher than the previously reported alpha-helices in an explicit solvent simulation with AMBER94 (AMBER94/TIP3P). Only AMBER96/GBSA shows a reasonable free energy landscape with the lowest free energy structure the native one despite an erroneous salt-bridge between D47 and K50. Detailed results on free energy contour maps, lowest free energy structures, distribution of native contacts, alpha-helical content during the folding process, NOE comparison with NMR, and temperature dependences are reported and discussed for all five models.     

Recent advances in the development and application of implicit solvent models in biomolecule simulations

Advances have recently been made in the development of implicit solvent methodologies and their application to the modeling of biomolecules, particularly with regard to generalized Born approaches, dielectric screening function formulations and models based on solvent-accessible surface areas. Interesting new developments include more refined non-polar solvation energy estimators, and implicit methods for modeling low-dielectric and heterogeneous environments such as membrane systems. These have been successfully applied to molecular dynamics simulations, the scoring of protein conformations, and the calculation of binding affinities and folding free energy landscapes.     

An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins

Exploiting recent developments in generalized Born (GB) electrostatics theory, we have reformulated the calculation of the self-electrostatic solvation energy to account for the influence of biological membranes. Consistent with continuum Poisson-Boltzmann (PB) electrostatics, the membrane is approximated as an solvent-inaccessible infinite planar low-dielectric slab. The present membrane GB model closely reproduces the PB electrostatic solvation energy profile across the membrane. The nonpolar contribution to the solvation energy is taken to be proportional to the solvent-exposed surface area (SA) with a phenomenological surface tension coefficient. The proposed membrane GB/SA model requires minor modifications of the pre-existing GB model and appears to be quite efficient. By combining this implicit model for the solvent/bilayer environment with advanced computational sampling methods, like replica-exchange molecular dynamics, we are able to fold and assemble helical membrane peptides. We examine the reliability of this model and approach by applications to three membrane peptides: melittin from bee venom, the transmembrane domain of the M2 protein from Influenza A (M2-TMP), and the transmembrane domain of glycophorin A (GpA). In the context of these proteins, we explore the role of biological membranes (represented as a low-dielectric medium) in affecting the conformational changes in melittin, the tilt of transmembrane peptides with respect to the membrane normal (M2-TMP), helix-to-helix interactions in membranes (GpA), and the prediction of the configuration of transmembrane helical bundles (GpA). The present method is found to perform well in each of these cases and is anticipated to be useful in the study of folding and assembly of membrane proteins as well as in structure refinement and modeling of membrane proteins where a limited number of experimental observables are available.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

An extended aqueous solvation model based on atom-weighted solvent accessible surface areas: SAWSA v2.0 model

A new method is proposed for calculating aqueous solvation free energy based on atom-weighted solvent accessible surface areas. The method, SAWSA v2.0, gives the aqueous solvation free energy by summing the contributions of component atoms and a correction factor. We applied two different sets of atom typing rules and fitting processes for small organic molecules and proteins, respectively. For small organic molecules, the model classified the atoms in organic molecules into 65 basic types and additionally. For small organic molecules we proposed a correction factor of hydrophobic carbon to account for the aggregation of hydrocarbons and compounds with long hydrophobic aliphatic chains. The contributions for each atom type and correction factor were derived by multivariate regression analysis of 379 neutral molecules and 39 ions with known experimental aqueous solvation free energies. Based on the new atom typing rules, the correlation coefficient (r) for fitting the whole neutral organic molecules is 0.984, and the absolute mean error is 0.40 kcal mol^–1, which is much better than those of the model proposed by Wang et al. and the SAWSA model previously proposed by us. Furthermore, the SAWSA v2.0 model was compared with the simple atom-additive model based on the number of atom types (NA). The calculated results show that for small organic molecules, the predictions from the SAWSA v2.0 model are slightly better than those from the atom-additive model based on NA. However, for macromolecules such as proteins, due to the connection between their molecular conformation and their molecular surface area, the atom-additive model based on the number of atom types has little predictive power. In order to investigate the predictive power of our model, a systematic comparison was performed on seven solvation models including SAWSA v2.0, GB/SA_1, GB/SA_2, PB/SA_1, PB/SA_2, AM1/SM5.2R and SM5.0R. The results showed that for organic molecules the SAWSA v2.0 model is better than the other six solvation models. For proteins, the model classified the atoms into 20 basic types and the predicted aqueous free energies of solvation by PB/SA were used for fitting. The solvation model based on the new parameters was employed to predict the solvation free energies of 38 proteins. The predicted values from our model were in good agreement with those from the PB/SA model and were much better than those given by the other four models developed for proteins.Figure The definition of hydrophobic carbons. Here CA, CB and CD are three carbon atoms; X represents a heteroatom. According to our definition, CB is a hydrophobic carbon, CA is not a hydrophobic carbon because a heteroatom is within four atoms and CD is not a hydrophobic carbon because CD is sp²- hydridized and in a six-member ring.Electronic Supplementary Material Supplementary material is available for this article at     

Parameterization of an effective potential for protein–ligand binding from host–guest affinity data

Force field accuracy is still one of the “stalemates” in biomolecular modeling. Model systems with high quality experimental data are valuable instruments for the validation and improvement of effective potentials. With respect to protein–ligand binding, organic host–guest complexes have long served as models for both experimental and computational studies because of the abundance of binding affinity data available for such systems. Binding affinity data collected for cyclodextrin (CD) inclusion complexes, a popular model for molecular recognition, is potentially a more reliable resource for tuning energy parameters than hydration free energy measurements. Convergence of binding free energy calculations on CD host–guest systems can also be obtained rapidly, thus offering the opportunity to assess the robustness of these parameters. In this work, we demonstrate how implicit solvent parameters can be developed using binding affinity experimental data and the binding energy distribution analysis method (BEDAM) and validated using the Grid Inhomogeneous Solvation Theory analysis. These new solvation parameters were used to study protein–ligand binding in two drug targets against the HIV‐1 virus and improved the agreement between the calculated and the experimental binding affinities. This work illustrates how benchmark sets of high quality experimental binding affinity data and physics‐based binding free energy models can be used to evaluate and optimize force fields for protein–ligand systems. Copyright © 2015 John Wiley & Sons, Ltd.     

Free energy surfaces of miniproteins with a betabetaalpha motif: replica exchange molecular dynamics simulation with an implicit solvation model

Designed miniproteins with a betabetaalpha motif, such as BBA5, 1FSD, and 1PSV can serve as a benchmark set to test the validity of all-atom force fields with computer simulation, because they contain all the basic structural elements in protein folding. Unfortunately, it was found that the standard all-atom force fields with the generalized Born (GB) implicit solvation model tend to produce distorted free energy surfaces for the betabetaalpha proteins, not only because energetically those proteins need to be described by more balanced weights of the alpha- and beta-strands, but also because the GB implicit solvation model suffers from overestimated salt bridge effects. In an attempt to resolve these problems, we have modified one of the standard all-atom force fields in conjunction with the GB model, such that each native state of the betabetaalpha proteins is in its free energy minimum state with reasonable energy barriers separating local minima. With this modified energy model, the free energy contour map in each protein was constructed from the replica exchange molecular dynamics REMD simulation. The resulting free energy surfaces are significantly improved in comparison with previous simulation results and consistent with general views on small protein folding behaviors with realistic topology and energetics of all three proteins.     

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Implicit solvation models provide, for many applications, a reasonably accurate and computationally effective way to describe the electrostatics of aqueous solvation. Here, a popular analytical Generalized Born (GB) solvation model is modified to improve its accuracy in calculating the solvent polarization part of free energy changes in large-scale conformational transitions, such as protein folding. In contrast to an earlier GB model (implemented in the AMBER-6 program), the improved version does not overstabilize the native structures relative to the finite-difference Poisson-Boltzmann continuum treatment. In addition to improving the energy balance between folded and unfolded conformers, the algorithm (available in the AMBER-7 and NAB molecular modeling packages) is shown to perform well in more than 50 ns of native-state molecular dynamics (MD) simulations of thioredoxin, protein-A, and ubiquitin, as well as in a simulation of Barnase/Barstar complex formation. For thioredoxin, various combinations of input parameters have been explored, such as the underlying gas-phase force fields and the atomic radii. The best performance is achieved with a previously proposed modification to the torsional potential in the Amber ff99 force field, which yields stable native trajectories for all of the tested proteins, with backbone root-mean-square deviations from the native structures being approximately 1.5 A after 6 ns of simulation time. The structure of Barnase/Barstar complex is regenerated, starting from an unbound state, to within 1.9 A relative to the crystal structure of the complex.     

Empirical solvation models can be used to differentiate native from near-native conformations of bovine pancreatic trypsin inhibitor.

Several hydration models for peptides and proteins based on solvent accessible surface area have been proposed previously. We have evaluated some of these models as well as four new ones in the context of near-native conformations of a protein. In addition, we propose an empirical site-site distance-dependent correction that can be used in conjunction with any of these models. The set of near-native structures consisted of 39 conformations of bovine pancreatic trypsin inhibitor (BPTI) each of which was a local minimum of an empirical energy function (ECEPP) in the absence of solvent. Root-mean-square (rms) deviations from the crystallographically determined structure were in the following ranges: 1.06-1.94 A for all heavy atoms, 0.77-1.36 A for all backbone heavy atoms, 0.68-1.33 A for all alpha-carbon atoms, and 1.41-2.72 A for all side-chain heavy atoms. We have found that there is considerable variation among the solvent models when evaluated in terms of concordance between the solvation free energy and the rms deviations from the crystallographically determined conformation. The solvation model for which the best concordance (0.939) with the rms deviations of the C alpha atoms was found was derived from NMR coupling constants of peptides in water combined with an exponential site-site distance dependence of the potential of mean force. Our results indicate that solvation free energy parameters derived from nonpeptide free energies of hydration may not be transferrable to peptides. Parameters derived from peptide and protein data may be more applicable to conformational analysis of proteins. A general approach to derive parameters for free energy of hydration from ensemble-averaged properties of peptides in solution is described.     

Recent advances in implicit solvent-based methods for biomolecular simulations

Implicit solvent-based methods play an increasingly important role in molecular modeling of biomolecular structure and dynamics. Recent methodological developments have mainly focused on the extension of the generalized Born (GB) formalism for variable dielectric environments and accurate treatment of nonpolar solvation. Extensive efforts in parameterization of GB models and implicit solvent force fields have enabled ab initio simulation of protein folding to native or near-native structures. Another exciting area that has benefited from the advances in implicit solvent models is the development of constant pH molecular dynamics methods, which have recently been applied to the calculations of protein pK(a) values and the studies of pH-dependent peptide and protein folding.     

DOCK 6: Combining techniques to model RNA–small molecule complexes

With an increasing interest in RNA therapeutics and for targeting RNA to treat disease, there is a need for the tools used in protein-based drug design, particularly DOCKing algorithms, to be extended or adapted for nucleic acids. Here, we have compiled a test set of RNA–ligand complexes to validate the ability of the DOCK suite of programs to successfully recreate experimentally determined binding poses. With the optimized parameters and a minimal scoring function, 70% of the test set with less than seven rotatable ligand bonds and 26% of the test set with less than 13 rotatable bonds can be successfully recreated within 2 Å heavy-atom RMSD. When DOCKed conformations are rescored with the implicit solvent models AMBER generalized Born with solvent-accessible surface area (GB/SA) and Poisson–Boltzmann with solvent-accessible surface area (PB/SA) in combination with explicit water molecules and sodium counterions, the success rate increases to 80% with PB/SA for less than seven rotatable bonds and 58% with AMBER GB/SA and 47% with PB/SA for less than 13 rotatable bonds. These results indicate that DOCK can indeed be useful for structure-based drug design aimed at RNA. Our studies also suggest that RNA-directed ligands often differ from typical protein–ligand complexes in their electrostatic properties, but these differences can be accommodated through the choice of potential function. In addition, in the course of the study, we explore a variety of newly added DOCK functions, demonstrating the ease with which new functions can be added to address new scientific questions.     

Optimizing the affinity and specificity of ligand binding with the inclusion of solvation effect

Solvation effect is an important factor for protein–ligand binding in aqueous water. Previous scoring function of protein–ligand interactions rarely incorporates the solvation model into the quantification of protein–ligand interactions, mainly due to the immense computational cost, especially in the structure‐based virtual screening, and nontransferable application of independently optimized atomic solvation parameters. In order to overcome these barriers, we effectively combine knowledge‐based atom–pair potentials and the atomic solvation energy of charge‐independent implicit solvent model in the optimization of binding affinity and specificity. The resulting scoring functions with optimized atomic solvation parameters is named as specificity and affinity with solvation effect (SPA‐SE). The performance of SPA‐SE is evaluated and compared to 20 other scoring functions, as well as SPA. The comparative results show that SPA‐SE outperforms all other scoring functions in binding affinity prediction and “native” pose identification. Our optimization validates that solvation effect is an important regulator to the stability and specificity of protein–ligand binding. The development strategy of SPA‐SE sets an example for other scoring function to account for the solvation effect in biomolecular recognitions. Proteins 2015; 83:1632–1642. © 2015 Wiley Periodicals, Inc.     

Does Flory–Huggins theory help in interpreting solute partitioning experiments?

The development of Flory–Huggins (FH) theory is reviewed, particularly with regard to the molecular significance of the interaction parameter that scales the contact interaction of solute and solvent. The chemical potential given by FH theory for an “idealute” solute is then compared with that provided by a more general, statistical thermodynamic approach. It is found that the FH contact term does not directly correspond to the solvation free energy. The significance of this result for the interpretation of free energies of transfer of a solute from one solvent to another is examined. It is shown that neither the earlier recommended standard free energy change for the process (using the infinitely dilute reference state, mole fraction units) nor the recently recommended FH-corrected standard free energy change provides the solvation energy desired. Instead, the standard free energy using the infinitely dilute reference state and molarity units, as long advocated by Ben-Naim, provides the desired solvation free energy. Correction of extant values, based on mole fraction units, is easily made. However, application of such results to problems of protein folding is not straightforward. © 1994 John Wiley & Sons, Inc.     

Systematic search for the rate constants that control the exocytotic process from chromaffin cells by a genetic algorithm

We have recently created a kinetic model that reproduces the dynamics of exocytosis with high accuracy. The reconstruction necessitated a search, in a multi-dimensional parameter space, for 37 parameters that described the system, with no assurance that the parameters, which reconstructed the observations, are a unique set. In the present study, a Genetic Algorithm (GA) was used for a thorough search in the unknown parameter space, using a strategy of gradual increase of the complexity of the analyzed input data. Upon systematic incorporation of one to four measurable parameters, used as input signals for the analysis, the constraint set on the GA search imposed the convergence of the free parameters into a single narrow range. The mean values for each adjustable parameter represent a minimum for the fitness function in the multi-dimensional parameter space. The GA search demonstrates that the parameters that control the kinetics of exocytosis are the rate constants of the steps downstream to synaptotagmin binding, and that the equilibrium constant of the binding of calcium to Munc13 controls the calcium-dependent priming process. Thus, the systematic use of the GA creates a link between specific reactions in the process of exocytosis and experimental phenotypes.     

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.     

Linear response theory: An alternative to PB and GB methods for the analysis of molecular dynamics trajectories?

We explore the use of classical Linear Response Theory (LRT) as an alternative strategy to the use of Molecular Mechanics/Poisson-Boltzmann strategies to compute the solvation free energy of macromolecules from molecular dynamics simulations using an explicit representation of solvent. The method reproduces well the free energy of solvation of standard amino acid side chains, small peptides, and proteins. The use of a fully discrete representation of solvent avoids the possible problems of continuum models to represent the solvation of systems containing tightly bound water molecules.     

Blind Test of Physics-Based Prediction of Protein Structures

We report here a multiprotein blind test of a computer method to predict native protein structures based solely on an all-atom physics-based force field. We use the AMBER 96 potential function with an implicit (GB/SA) model of solvation, combined with replica-exchange molecular-dynamics simulations. Coarse conformational sampling is performed using the zipping and assembly method (ZAM), an approach that is designed to mimic the putative physical routes of protein folding. ZAM was applied to the folding of six proteins, from 76 to 112 monomers in length, in CASP7, a community-wide blind test of protein structure prediction. Because these predictions have about the same level of accuracy as typical bioinformatics methods, and do not utilize information from databases of known native structures, this work opens up the possibility of predicting the structures of membrane proteins, synthetic peptides, or other foldable polymers, for which there is little prior knowledge of native structures. This approach may also be useful for predicting physical protein folding routes, non-native conformations, and other physical properties from amino acid sequences.