Comparative Study of the Folding Free Energy Landscape of a Three-Stranded β-Sheet Protein with Explicit and Implicit Solvent Models

Solvation Thermodynamics of Solutes in Water and Ionic Liquids Using the Multiscale Solvation-Layer Interface Condition Continuum Model

Molecular assembly processes are generally driven by thermodynamic properties in solutions. Atomistic modeling can be very helpful in designing and understanding complex systems, except that bulk solvent is very inefficient to treat explicitly as discrete molecules. In this work, we develop and assess two multiscale solvation models for computing solvation thermodynamic properties. The new SLIC/CDC model combines continuum solvent electrostatics based on the solvent layer interface condition (SLIC) with new statistical thermodynamic models for hydrogen bonding and nonpolar modes: cavity formation, dispersion interactions, combinatorial mixing (CDC). Given the structures of 500 solutes, the SLIC/CDC model predicts Gibbs energies of solvation in water with an average accuracy better than 1 kcal/mol, when compared to experimental measurements, and better than 0.8 kcal/mol, when compared to explicit-solvent molecular dynamics simulations. The individual SLIC/CDC energy mode values agree quantitatively with those computed from explicit-solvent molecular dynamics. The previously published SLIC/SASA multiscale model combines the SLIC continuum electrostatic model with the solvent-accessible surface area (SASA) nonpolar energy mode. With our new, improved parametrization method, the SLIC/SASA model now predicts Gibbs energies of solvation with better than 1.4 kcal/mol average accuracy in aqueous systems, compared to experimental and explicit-solvent molecular dynamics, and better than 1.6 kcal/mol average accuracy in ionic liquids, compared to explicit-solvent molecular dynamics. Both models predict solvation entropies, and are the first implicit-solvation models capable of predicting solvation heat capacities.

Sampling of the conformational landscape of small proteins with Monte Carlo methods

Computer simulation provides an increasingly realistic picture of large-scale conformational change of proteins, but investigations remain fundamentally constrained by the femtosecond timestep of molecular dynamics simulations. For this reason, many biologically interesting questions cannot be addressed using accessible state-of-the-art computational resources. Here, we report the development of an all-atom Monte Carlo approach that permits the modelling of the large-scale conformational change of proteins using standard off-the-shelf computational hardware and standard all-atom force fields. We demonstrate extensive thermodynamic characterization of the folding process of the α-helical Trp-cage, the Villin headpiece and the β-sheet WW-domain. We fully characterize the free energy landscape, transition states, energy barriers between different states, and the per-residue stability of individual amino acids over a wide temperature range. We demonstrate that a state-of-the-art intramolecular force field can be combined with an implicit solvent model to obtain a high quality of the folded structures and also discuss limitations that still remain.

Assessing AMBER force fields for protein folding in an implicit solvent

Molecular dynamics (MD) simulation implemented with a state-of-the-art protein force field and implicit solvent model is an attractive approach to investigate protein folding, one of the most perplexing problems in molecular biology. But how well can force fields developed independently of implicit solvent models work together in reproducing diverse protein native structures and measuring the corresponding folding thermodynamics is not always clear. In this work, we performed enhanced sampling MD simulations to assess the ability of six AMBER force fields (FF99SBildn, FF99SBnmr, FF12SB, FF14ipq, FF14SB, and FF14SBonlysc) as coupled with a recently improved pair-wise GB-Neck2 model in modeling the folding of two helical and two β-sheet peptides. Whilst most of the tested force fields can yield roughly similar features for equilibrium conformational ensembles and detailed folding free-energy profiles for short α-helical TC10b in an implicit solvent, the measured counterparts are significantly discrepant in the cases of larger or β-structured peptides (HP35, 1E0Q, and GTT). Additionally, the calculated folding/unfolding thermodynamic quantities can only partially match the experimental data. Although a combination of the force fields and GB-Neck2 implicit model able to describe all aspects of the folding transitions towards the native structures of all the considered peptides was not identified, we found that FF14SBonlysc coupled with the GB-Neck2 model seems to be a reasonably balanced combination to predict peptide folding preferences.

Kinetics from Implicit Solvent Simulations of Biomolecules as a Function of Viscosity

Kinetic properties of alanine dipeptide, the B1 domain of streptococcal protein G, and ubiquitin are compared between explicit solvent and implicit solvent simulations with the generalized Born molecular volume (GBMV) method. The results indicate that kinetics from explicit solvent simulations and experiments can be matched closely when the implicit solvent simulations are combined with Langevin dynamics and a friction coefficient near 10 ps(-1). Smaller and larger friction coefficients accelerate and slow down conformational sampling. It is found that local conformational exploration without the crossing of significant barriers can be accelerated by a factor of 4-5; however, the acceleration of barrier crossings is limited to about a factor of 2. The use of a Nosé-Hoover thermostat instead of Langevin dynamics greatly enhances local conformational sampling but slows down the crossing of barriers by at least an order of magnitude because of the lack of solute-solvent stochastic collisions.

Computer Simulations of Peptides from the p53 DNA Binding Domain

We have studied the dynamics and thermodynamics of two of the four evolutionarily conserved segments from the p53 DNA binding domain, using molecular dynamics and replica exchange simulations. These two regions contain well-defined elements of secondary structure (a β hairpin for region II and an α helix for region V) and bind to DNA in the intact protein. They are also mutational hot spots. The goal of our study was to determine the stability and folding propensity of these peptides in isolation. We used three force fields and solvent models (CHARMM19 with EEF1, CHARMM27 with GBMV, GROMOS96 with SPC). The predicted stability, folding propensity, and secondary structures depend upon the potential. Secondary structure predictors identify helical propensity for region II, in agreement with one of the force fields (CHARMM/GBMV). However, the other two potentials favor β structure for this peptide, although the conformations may differ from the crystal. For region V secondary structure predictions are unclear. Only one force field (CHARMM/GBMV) produces low-lying free energy minima that retain some of the α helical structure from the crystal structure. The other two potentials appear to favor β structure for this peptide.

Prediction of the structure of a silk-like protein in oligomeric states using explicit and implicit solvent models

We perform Replica Exchange Molecular Dynamics (REMD) simulations on a silk-like protein design with amino-acid sequence [(Gly-Ala)3-Gly-Glu]5 to investigate the stability of a single protein, a dimer, a trimer and a tetramer made up of these proteins starting from β-roll and β-sheet structures in both explicit (TIP3P) and implicit (GBSA) solvent models. Our simulation results for the implicit solvent model agree with those for the explicit solvent model for simulation times up to the longest tested, being 30 ns per replica. From this we infer that the implicit solvent model that we use is reliable, allowing us to reach much longer time scales (up to 200 ns per replica). We find that the self-assembly of fibers of these proteins in solution must be a nucleated process, involving nuclei made up of at least three monomers. We also find that the conformation of the protein changes upon assembly, i.e., there is a transition from a disordered globular state to an ordered β-sheet structure in the self-assembled state of aggregates containing more than two monomers. This indicates that autosteric effects must be important in the polymerization of this protein, reminiscent of what is observed for β-amyloids. Our findings are consistent with recent experimental results on a protein with an amino acid sequence similar to that of the protein we study.

Effect of sequence and stereochemistry reversal on p53 peptide mimicry

Peptidomimetics effective in modulating protein-protein interactions and resistant to proteolysis have potential in therapeutic applications. An appealing yet underperforming peptidomimetic strategy is to employ D-amino acids and reversed sequences to mimic a lead peptide conformation, either separately or as the combined retro-inverso peptide. In this work, we examine the conformations of inverse, reverse and retro-inverso peptides of p53(15-29) using implicit solvent molecular dynamics simulation and circular dichroism spectroscopy. In order to obtain converged ensembles for the peptides, we find enhanced sampling is required via the replica exchange molecular dynamics method. From these replica exchange simulations, the D-peptide analogues of p53(15-29) result in a predominantly left-handed helical conformation. When the parent sequence is reversed sequence as either the L-peptide and D-peptide, these peptides display a greater helical propensity, feature reflected by NMR and CD studies in TFE/water solvent. The simulations also indicate that, while approximately similar orientations of the side-chains are possible by the peptide analogues, their ability to mimic the parent peptide is severely compromised by backbone orientation (for D-amino acids) and side-chain orientation (for reversed sequences). A retro-inverso peptide is disadvantaged as a mimic in both aspects, and further chemical modification is required to enable this concept to be used fruitfully in peptidomimetic design. The replica exchange molecular simulation approach adopted here, with its ability to provide detailed conformational insights into modified peptides, has potential as a tool to guide structure-based design of new improved peptidomimetics.

Structural modelling and dynamics of proteins for insights into drug interactions

Proteins are the workhorses of biomolecules and their function is affected by their structure and their structural rearrangements during ligand entry, ligand binding and protein-protein interactions. Hence, the knowledge of protein structure and, importantly, the dynamic behaviour of the structure are critical for understanding how the protein performs its function. The predictions of the structure and the dynamic behaviour can be performed by combinations of structure modelling and molecular dynamics simulations. The simulations also need to be sensitive to the constraints of the environment in which the protein resides. Standard computational methods now exist in this field to support the experimental effort of solving protein structures. This review presents a comprehensive overview of the basis of the calculations and the well-established computational methods used to generate and understand protein structure and function and the study of their dynamic behaviour with the reference to lung-related targets.

ABEEMσπ FLUCTUATING CHARGE FORCE FIELD APPLIED TO ALANINE DIPEPTIDE AND ALANINE DIPEPTIDE–WATER SYSTEMS

Atom-bond electronegativity equalization method at σπ level fused into molecular mechanics (ABEEMσπ/MM) divides the bond regions into σ and π bond regions on the basis of previous ABEEM/MM. It may suitably reflect intramolecular and intermolecular interaction and polarization. The fitting function k H-bond in the hydrogen bond (HB) interaction region increases the capability of ABEEMσπ/MM to simulate the hydration. Hydration of alanine dipeptide (AD) in aqueous solution is determined by the intramolecular and intermolecular HBs and the competition among the molecular packing effects. The acceptor molecule in HB complex contains at least one pair of lone pair electrons, sometimes contains π bonds, whose orientations directly effect the orientation of HBs. Therefore, ABEEMσπ/MM has obviously predominance to discuss the AD and AD–water systems, which contain many lone pair electrons, π bonds, and abundant HB nets. Properties of six AD conformers, clusters AD +( H2O )1–4 obtained from ABEEMσπ/MM agree well with the results of experiments, ab initio and other force fields. Structural and dynamical properties of the hydration water molecules have just embodied that the ABEEMσπ/MM gives correct hydration description relative to other force fields.

Molecular dynamics of ribosomal elongation factors G and Tu

Translation on the ribosome is controlled by external factors. During polypeptide lengthening, elongation factors EF-Tu and EF-G consecutively interact with the bacterial ribosome. EF-Tu binds and delivers an aminoacyl-tRNA to the ribosomal A site and EF-G helps translocate the tRNAs between their binding sites after the peptide bond is formed. These processes occur at the expense of GTP. EF-Tu:tRNA and EF-G are of similar shape, share a common binding site, and undergo large conformational changes on interaction with the ribosome. To characterize the internal motion of these two elongation factors, we used 25 ns long all-atom molecular dynamics simulations. We observed enhanced mobility of EF-G domains III, IV, and V and of tRNA in the EF-Tu:tRNA complex. EF-Tu:GDP complex acquired a configuration different from that found in the crystal structure of EF-Tu with a GTP analogue, showing conformational changes in the switch I and II regions. The calculated electrostatic properties of elongation factors showed no global similarity even though matching electrostatic surface patches were found around the domain I that contacts the ribosome, and in the GDP/GTP binding region.

Development of CHARMM polarizable force field for nucleic acid bases based on the classical Drude oscillator model

A polarizable force field for nucleic acid bases based on the classical Drude oscillator model is presented. Parameter optimization was performed to reproduce crystallographic geometries, crystal unit cell parameters, heats of sublimation, vibrational frequencies and assignments, dipole moments, molecular polarizabilities and quantum mechanical base-base and base-water interaction energies. The training and validation data included crystals of unsubstituted and alkyl-substituted adenine, guanine, cytosine, uracil, and thymine bases, hydrated crystals, and hydrogen bonded base pairs. Across all compounds, the RMSD in the calculated heats of sublimation is 4.1%. This equates to an improvement of more than 2.5 kcal/mol in accuracy compared to the nonpolarizable CHARMM27 force field. However, the level of agreement with experimental molecular volume decreased from 1.7% to 2.1% upon moving from the nonpolarizable to the polarizable model. The representation of dipole moments is significantly improved with the Drude polarizable force field. Unlike in additive force fields, there is no requirement for the gas-phase dipole moments to be overestimated, illustrating the ability of the Drude polarizable force field to treat accurately differently dielectric environments and indicating the improvements in the electrostatic model. Validation of the model was performed on the basis of the calculation of the gas-phase binding enthalpies of base pairs obtained via potential of mean force calculations; the additive and polarizable models both performed satisfactorily with average differences of 0.2 and 0.9 kcal/mol, respectively, and rms differences of 1.3 and 1.7 kcal/mol, respectively. Overall, considering the number of significant improvements versus the additive CHARMM force field, the incorporation of explicit polarizability into the force field for nucleic acid bases represents an additional step toward accurate computational modeling of biological systems.

Modelling proteins: conformational sampling and reconstruction of folding kinetics

In the last decades biomolecular simulation has made tremendous inroads to help elucidate biomolecular processes in-silico. Despite enormous advances in molecular dynamics techniques and the available computational power, many problems involve long time scales and large-scale molecular rearrangements that are still difficult to sample adequately. In this review we therefore summarise recent efforts to fundamentally improve this situation by decoupling the sampling of the energy landscape from the description of the kinetics of the process. Recent years have seen the emergence of many advanced sampling techniques, which permit efficient characterisation of the relevant family of molecular conformations by dispensing with the details of the short-term kinetics of the process. Because these methods generate thermodynamic information at best, they must be complemented by techniques to reconstruct the kinetics of the process using the ensemble of relevant conformations. Here we review recent advances for both types of methods and discuss their perspectives to permit efficient and accurate modelling of large-scale conformational changes in biomolecules. This article is part of a Special Issue entitled: Protein Dynamics: Experimental and Computational Approaches.

Electrostatic solvation energy for two oppositely charged ions in a solvated protein system: salt bridges can stabilize proteins

Born-type electrostatic continuum methods have been an indispensable ingredient in a variety of implicit-solvent methods that reduce computational effort by orders of magnitude compared to explicit-solvent MD simulations and thus enable treatment using larger systems and/or longer times. An analysis of the limitations and failures of the Born approaches serves as a guide for fundamental improvements without diminishing the importance of prior works. One of the major limitations of the Born theory is the lack of a liquidlike description of the response of solvent dipoles to the electrostatic field of the solute and the changes therein, a feature contained in the continuum Langevin-Debye (LD) model applied here to investigate how Coulombic interactions depend on the location of charges relative to the protein/water boundary. This physically more realistic LD model is applied to study the stability of salt bridges. When compared head to head using the same (independently measurable) physical parameters (radii, dielectric constants, etc.), the LD model is in good agreement with observations, whereas the Born model is grossly in error. Our calculations also suggest that a salt bridge on the protein's surface can be stabilizing when the charge separation is < or =4 A.

Modeling of peptides connecting the ligand-binding and transmembrane domains in the GluR2 glutamate receptor

Ligand-gated Glutamate receptors (GluR) mediate synaptic signals in the nervous system. Ionotropic GluRs of AMPA type, the subject of this study, are tetrameric assemblies of monomer subunits, each of which is constructed in a modular fashion from functional subdomains. The extracellular ligand-binding domain (LBD) changes its conformation upon binding of an agonist ligand followed by opening of a transmembrane (TM) ion channel. Peptides connecting the LBD and TM domains facilitate gating of the channel, and their structure and composition are important for the receptor functioning. In this study, we used replica exchange molecular dynamics (REMD) simulations to model S1M1 and S2M3 connecting peptides of the GluR2 receptor in two implicit solvents, water and interfacial water/lipid medium characterized by lower polarity. Propensity of these peptides to form helical structures was analyzed using helicity measure derived from the free energy of the simulated ensembles of structures. The S1M1 and S2M3 connecting peptides were not helical in our simulations in both dielectric environments in the absence of the rest of the protein. The structures of the LBD fragment with known high-resolution alpha-helical structure and of the TM3 helix were successfully predicted in the simulations, which in part validate our results. The S2M3 peptide, which is important in gating, formed a well-defined coil structure and salt-bridges with the S2 domain. The S1M1 peptide formed a loop structure via formation of internal salt-bridges. Potential implications of these structures on function of the receptor are discussed.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Bounding the electrostatic free energies associated with linear continuum models of molecular solvation

The importance of electrostatic interactions in molecular biology has driven extensive research toward the development of accurate and efficient theoretical and computational models. Linear continuum electrostatic theory has been surprisingly successful, but the computational costs associated with solving the associated partial differential equations (PDEs) preclude the theory's use in most dynamical simulations. Modern generalized-Born models for electrostatics can reproduce PDE-based calculations to within a few percent and are extremely computationally efficient but do not always faithfully reproduce interactions between chemical groups. Recent work has shown that a boundary-integral-equation formulation of the PDE problem leads naturally to a new approach called boundary-integral-based electrostatics estimation (BIBEE) to approximate electrostatic interactions. In the present paper, we prove that the BIBEE method can be used to rigorously bound the actual continuum-theory electrostatic free energy. The bounds are validated using a set of more than 600 proteins. Detailed numerical results are presented for structures of the peptide met-enkephalin taken from a molecular-dynamics simulation. These bounds, in combination with our demonstration that the BIBEE methods accurately reproduce pairwise interactions, suggest a new approach toward building a highly accurate yet computationally tractable electrostatic model.

Langevin-Debye model for nonlinear electrostatic screening of solvated ions

Ion-ion interactions are central to numerous phenomena in geology, biology, and material science, but current understanding is based on linear theories of limited physical applicability. Using the Langevin-Debye model, we investigate how nonlinear dielectric saturation alters screening of electrostatic interactions between solvated charges and show that charge screening in liquids strongly depends on the charge magnitudes but negligibly on the ion radii. Qualitatively different universal screening functions are predicted for ions with charges of the same or opposite signs.

Interpreting the Coulomb-field approximation for generalized-Born electrostatics using boundary-integral equation theory

The importance of molecular electrostatic interactions in aqueous solution has motivated extensive research into physical models and numerical methods for their estimation. The computational costs associated with simulations that include many explicit water molecules have driven the development of implicit-solvent models, with generalized-Born (GB) models among the most popular of these. In this paper, we analyze a boundary-integral equation interpretation for the Coulomb-field approximation (CFA), which plays a central role in most GB models. This interpretation offers new insights into the nature of the CFA, which traditionally has been assessed using only a single point charge in the solute. The boundary-integral interpretation of the CFA allows the use of multiple point charges, or even continuous charge distributions, leading naturally to methods that eliminate the interpolation inaccuracies associated with the Still equation. This approach, which we call boundary-integral-based electrostatic estimation by the CFA (BIBEE/CFA), is most accurate when the molecular charge distribution generates a smooth normal displacement field at the solute-solvent boundary, and CFA-based GB methods perform similarly. Conversely, both methods are least accurate for charge distributions that give rise to rapidly varying or highly localized normal displacement fields. Supporting this analysis are comparisons of the reaction-potential matrices calculated using GB methods and boundary-element-method (BEM) simulations. An approximation similar to BIBEE/CFA exhibits complementary behavior, with superior accuracy for charge distributions that generate rapidly varying normal fields and poorer accuracy for distributions that produce smooth fields. This approximation, BIBEE by preconditioning (BIBEE/P), essentially generates initial guesses for preconditioned Krylov-subspace iterative BEMs. Thus, iterative refinement of the BIBEE/P results recovers the BEM solution; excellent agreement is obtained in only a few iterations. The boundary-integral-equation framework may also provide a means to derive rigorous results explaining how the empirical correction terms in many modern GB models significantly improve accuracy despite their simple analytical forms.

Mechanical stability of helical beta-peptides and a comparison of explicit and implicit solvent models

Synthetic beta-peptide oligomers have been shown to form stable folded structures analogous to those encountered in naturally occurring proteins. Literature studies have speculated that the conformational stability of beta-peptides is greater than that of alpha-peptides. Direct measurements of that stability, however, are not available. Molecular simulations are used in this work to quantify the mechanical stability of four helical beta-peptides. This is achieved by subjecting the molecules to tension. The potential of mean force associated with the resulting unfolding process is determined using both an implicit and an explicit solvent model. It is found that all four molecules exhibit a highly stable helical structure. It is also found that the energetic contributions to the potential of mean force do not change appreciably when the molecules are stretched in explicit water. In contrast, the entropic contributions decrease significantly. As the peptides unfold, a loss of intramolecular energy is compensated by the formation of additional water-peptide hydrogen bonds. These entropic effects lead in some cases to a loss of stability upon cooling the peptides, a phenomenon akin to the cold denaturing of some proteins. While the location of the free energy minimum and the structural helicity of the peptides are comparable in the implicit-solvent and explicit-water cases, it is found that, in general, the helical structure of the molecules is more stable in the implicit solvent model than in explicit water.

Prediction of protein loop conformations using multiscale modeling methods with physical energy scoring functions

This article examines ab initio methods for the prediction of protein loops by a computational strategy of multiscale conformational sampling and physical energy scoring functions. Our approach consists of initial sampling of loop conformations from lattice-based low-resolution models followed by refinement using all-atom simulations. To allow enhanced conformational sampling, the replica exchange method was implemented. Physical energy functions based on CHARMM19 and CHARMM22 parameterizations with generalized Born (GB) solvent models were applied in scoring loop conformations extracted from the lattice simulations and, in the case of all-atom simulations, the ensemble of conformations were generated and scored with these models. Predictions are reported for 25 loop segments, each eight residues long and taken from a diverse set of 22 protein structures. We find that the simulations generally sampled conformations with low global root-mean-square-deviation (RMSD) for loop backbone coordinates from the known structures, whereas clustering conformations in RMSD space and scoring detected less favorable loop structures. Specifically, the lattice simulations sampled basins that exhibited an average global RMSD of 2.21 +/- 1.42 A, whereas clustering and scoring the loop conformations determined an RMSD of 3.72 +/- 1.91 A. Using CHARMM19/GB to refine the lattice conformations improved the sampling RMSD to 1.57 +/- 0.98 A and detection to 2.58 +/- 1.48 A. We found that further improvement could be gained from extending the upper temperature in the all-atom refinement from 400 to 800 K, where the results typically yield a reduction of approximately 1 A or greater in the RMSD of the detected loop. Overall, CHARMM19 with a simple pairwise GB solvent model is more efficient at sampling low-RMSD loop basins than CHARMM22 with a higher-resolution modified analytical GB model; however, the latter simulation method provides a more accurate description of the all-atom energy surface, yet demands a much greater computational cost.

Benchmarking implicit solvent folding simulations of the amyloid beta(10-35) fragment

A pathogenetic feature of Alzhemier disease is the aggregation of monomeric beta-amyloid proteins (Abeta) to form oligomers. Usually these oligomers of long peptides aggregate on time scales of microseconds or longer, making computational studies using atomistic molecular dynamics models prohibitively expensive and making it essential to develop computational models that are cheaper and at the same time faithful to physical features of the process. We benchmark the ability of our implicit solvent model to describe equilibrium and dynamic properties of monomeric Abeta(10-35) using all-atom Langevin dynamics (LD) simulations, since Alphabeta(10-35) is the only fragment whose monomeric properties have been measured. The accuracy of the implicit solvent model is tested by comparing its predictions with experiment and with those from a new explicit water MD simulation, (performed using CHARMM and the TIP3P water model) which is approximately 200 times slower than the implicit water simulations. The dependence on force field is investigated by running multiple trajectories for Alphabeta(10-35) using the CHARMM, OPLS-aal, and GS-AMBER94 force fields, whereas the convergence to equilibrium is tested for each force field by beginning separate trajectories from the native NMR structure, a completely stretched structure, and from unfolded initial structures. The NMR order parameter, S2, is computed for each trajectory and is compared with experimental data to assess the best choice for treating aggregates of Alphabeta. The computed order parameters vary significantly with force field. Explicit and implicit solvent simulations using the CHARMM force fields display excellent agreement with each other and once again support the accuracy of the implicit solvent model. Alphabeta(10-35) exhibits great flexibility, consistent with experiment data for the monomer in solution, while maintaining a general strand-loop-strand motif with a solvent-exposed hydrophobic patch that is believed to be important for aggregation. Finally, equilibration of the peptide structure requires an implicit solvent LD simulation as long as 30 ns.

Prediction of Protein Loop Conformations using the AGBNP Implicit Solvent Model and Torsion Angle Sampling

The OPLS-AA all-atom force field and the Analytical Generalized Born plus Non-Polar (AGBNP) implicit solvent model, in conjunction with torsion angle conformational search protocols based on the Protein Local Optimization Program (PLOP), are shown to be effective in predicting the native conformations of 57 9-residue and 35 13-residue loops of a diverse series of proteins with low sequence identity. The novel nonpolar solvation free energy estimator implemented in AGBNP augmented by correction terms aimed at reducing the occurrence of ion pairing are important to achieve the best prediction accuracy. Extended versions of the previously developed PLOP-based conformational search schemes based on calculations in the crystal environment are reported that are suitable for application to loop homology modeling without the crystal environment. Our results suggest that in general the loop backbone conformation is not strongly influenced by crystal packing. The application of the temperature Replica Exchange Molecular Dynamics (T-REMD) sampling method for a few examples where PLOP sampling is insufficient are also reported. The results reported indicate that the OPLS-AA/AGBNP effective potential is suitable for high-resolution modeling of proteins in the final stages of homology modeling and/or protein crystallographic refinement.

Combination of molecular dynamics method and 3D-RISM theory for conformational sampling of large flexible molecules in solution

We have developed an algorithm for sampling the conformational space of large flexible molecules in solution, which combines the molecular dynamics (MD) method and the three-dimensional reference interaction site model (3D-RISM) theory. The solvent-induced force acting on solute atoms was evaluated as the gradient of the solvation free energy with respect to the solute-atom coordinates. To enhance the computation speed, we have applied a multiple timestep algorithm based on the RESPA (Reversible System Propagator Algorithm) to the combined MD/3D-RISM method. By virtue of the algorithm, one can choose a longer timestep for renewing the solvent-induced force compared with that of the conformational update. To illustrate the present MD/3D-RISM simulation, we applied the method to a model of acetylacetone in aqueous solution. The multiple timestep algorithm succeeded in enhancing the computation speed by 3.4 times for this model case. Acetylacetone possesses an intramolecular hydrogen-bonding capability between the hydroxyl group and the carbonyl oxygen atom, and the molecule is significantly stabilized due to this hydrogen bond, especially in gas phase. The intramolecular hydrogen bond was kept intact during almost entire course of the MD simulation in gas phase, while in the aqueous solutions the bond is disrupted in a significant number of conformations. This result qualitatively agrees with the behavior on a free energy barrier lying upon the process for rotating a torsional degree of freedom of the hydroxyl group, where it is significantly reduced in aqueous solution by a cancellation between the electrostatic interaction and the solvation free energy.

Solvation effect on conformations of 1,2:dimethoxyethane: charge-dependent nonlinear response in implicit solvent models

The physical content of and, in particular, the nonlinear contributions from the Langevin-Debye model are illustrated using two applications. First, we provide an improvement in the Langevin-Debye model currently used in some implicit solvent models for computer simulations of solvation free energies of small organic molecules, as well as of biomolecular folding and binding. The analysis is based on the implementation of a charge-dependent Langevin-Debye (qLD) model that is modified by subsequent corrections due to Onsager and Kirkwood. Second, the physical content of the model is elucidated by discussing the general treatment within the LD model of the self-energy of a charge submerged in a dielectric medium for three different limiting conditions and by considering the nonlinear response of the medium. The modified qLD model is used to refine an implicit solvent model (previously applied to protein dynamics). The predictions of the modified implicit solvent model are compared with those from explicit solvent molecular dynamics simulations for the equilibrium conformational populations of 1,2-dimethoxyethane (DME), which is the shortest ether molecule to reproduce the local conformational properties of polyethylene oxide, a polymer with tremendous technological importance and a wide variety of applications. Because the conformational population preferences of DME change dramatically upon solvation, DME is a good test case to validate our modified qLD model. The present analysis of the modified qLD model provides the motivation and tools for studying a wide variety of other interesting systems with heterogeneous dielectric properties and spatial anisotropy.

Dehydration-driven solvent exposure of hydrophobic surfaces as a driving force in peptide folding

Recent work has shown that the nature of hydration of pure hydrophobic surfaces changes with the length scale considered: water hydrogen-bonding networks adapt to small exposed hydrophobic species, hydrating or "wetting" them at relatively high densities, whereas larger hydrophobic areas are "dewetted" [Chandler D (2005), Nature 29:640-647]. Here we determine whether this effect is also present in peptides by examining the folding of a beta-hairpin (the 14-residue amyloidogenic prion protein H1 peptide), using microsecond time-scale molecular dynamics simulations. Two simulation models are compared, one explicitly including the water molecules, which may thus adapt locally to peptide configurations, and the other using a popular continuum approximation, the generalized Born/surface area implicit solvent model. The results obtained show that, in explicit solvent, peptide conformers with high solvent-accessible hydrophobic surface area indeed also have low hydration density around hydrophobic residues, whereas a concomitant higher hydration density around hydrophilic residues is observed. This dewetting effect stabilizes the fully folded beta-hairpin state found experimentally. In contrast, the implicit solvent model destabilizes the fully folded hairpin, tending to cluster hydrophobic residues regardless of the size of the exposed hydrophobic surface. Furthermore, the rate of the conformational transitions in the implicit solvent simulation is almost doubled with respect to that of the explicit solvent. The results suggest that dehydration-driven solvent exposure of hydrophobic surfaces may be a significant factor determining peptide conformational equilibria.

Implicit solvent simulations of DNA and DNA-protein complexes: agreement with explicit solvent vs experiment

Molecular dynamics simulations of biomolecules with implicit solvent reduce the computational cost and complexity of such simulations so that longer time scales and larger system sizes can be reached. While implicit solvent simulations of proteins have become well established, the success of implicit solvent in the simulation of nucleic acids has not been fully established to date. Results obtained in this study demonstrate that stable and efficient simulations of DNA and a protein-DNA complex can be achieved with an implicit solvent model based on continuum dielectric electrostatics. Differences in conformational sampling of DNA with two sets of atomic radii that are used to define the dielectric interface between the solute and the continuum dielectric model of the solvent are investigated. Results suggest that depending on the choice of atomic radii agreement is either closer to experimental data or to explicit solvent simulations. Furthermore, partial conformational transitions toward A-DNA conformations when salt is added within the implicit solvent framework are observed.

Elastic bag model for molecular dynamics simulations of solvated systems: application to liquid water and solvated peptides

The fluctuating elastic boundary (FEB) model for molecular dynamics has recently been developed and validated through simulations of liquid argon. In the FEB model, a flexible boundary which consists of particles connected by springs is used to confine the solvated system, thereby eliminating the need for periodic boundary conditions. In this study, we extend this model to the simulation of bulk water and solvated alanine dipeptide. Both the confining potential and boundary particle interaction functions are modified to preserve the structural integrity of the boundary and prevent the leakage of the solute-solvent system through the boundary. A broad spectrum of structural and dynamic properties of liquid water are computed and compared with those obtained from conventional periodic boundary condition simulations. The applicability of the model to biomolecular simulations is investigated through the analysis of conformational population distribution of solvated alanine dipeptide. In most cases we find remarkable agreement between the two simulation approaches.

Computational sidechain placement and protein mutagenesis with implicit solvent models

Structure prediction and computational protein design should benefit from accurate solvent models. We have applied implicit solvent models to two problems that are central to this area. First, we performed sidechain placement for 29 proteins, using a solvent model that combines a screened Coulomb term with an Accessible Surface Area term (CASA model). With optimized parameters, the prediction quality is comparable with earlier work that omitted electrostatics and solvation altogether. Second, we computed the stability changes associated with point mutations involving ionized sidechains. For over 1000 mutations, including many fully or partly buried positions, we compared CASA and two generalized Born models (GB) with a more accurate model, which solves the Poisson equation of continuum electrostatics numerically. CASA predicts the correct sign and order of magnitude of the stability change for 81% of the mutations, compared to 97% with the best GB. We also considered 140 mutations for which experimental data are available. Comparing to experiment requires additional assumptions about the unfolded protein structure, protein relaxation in response to the mutations, and contributions from the hydrophobic effect. With a simple, commonly-used unfolded state model, the mean unsigned error is 2.1 kcal/mol with both CASA and the best GB. Overall, the electrostatic model is not important for sidechain placement; CASA and GB are equivalent for surface mutations, while GB is far superior for fully or partly buried positions. Thus, for problems like protein design that involve all these aspects, the most recent GB models represent an important step forward. Along with the recent discovery of efficient, pairwise implementations of GB, this will open new possibilities for the computational engineering of proteins.

Evaluation of Poisson solvation models using a hybrid explicit/implicit solvent method

Implicit solvent methods have become popular tools in the field of protein dynamics simulations, yet evaluation of their validity has been primarily limited to comparisons with experimental and theoretical data for small molecules. In this paper, we use a recently developed hybrid explicit/implicit solvent methodology to evaluate the accuracy of several Poisson-based implicit solvent models. Specifically, we focus on the calculation of electrostatic solvation free energies of various fixed conformations for two proteins. We show that, among various dielectric boundary definitions, the Lee-Richards molecular surface has the best agreement with hybrid solvent results. Furthermore, certain modifications of the molecular surface Poisson protocol provide varied results. For instance, simple modifications of atomic radii on charged residues generally improve absolute errors but do not significantly reduce relative errors among conformations. On the other hand, using a water-probe radius of 1.0 A, as opposed to the standard value of 1.4 A, to generate the molecular surface, moderately improves both absolute and relative results.

Protein-Folding Dynamics: Overview of Molecular Simulation Techniques

Molecular dynamics (MD) is an invaluable tool with which to study protein folding in silico. Although just a few years ago the dynamic behavior of a protein molecule could be simulated only in the neighborhood of the experimental conformation (or protein unfolding could be simulated at high temperature), the advent of distributed computing, new techniques such as replica-exchange MD, new approaches (based on, e.g., the stochastic difference equation), and physics-based reduced models of proteins now make it possible to study protein-folding pathways from completely unfolded structures. In this review, we present algorithms for MD and their extensions and applications to protein-folding studies, using all-atom models with explicit and implicit solvent as well as reduced models of polypeptide chains.

Secondary structure provides a template for the folding of nearby polypeptides

Although protein structures are primarily encoded by their sequences, they are also critically dependent on environmental factors such as solvents and interactions with other molecules. Here we investigate how the folding-energy landscape of a short peptide is altered by interactions with another peptide, by performing atomistic replica-exchange molecular dynamics simulations of polyalanines in various environments. We analyzed the free-energy landscapes of Ala7 and Ala8 in isolation, near an alpha-helix template, and near a beta-strand template. The isolated Ala7 and Ala8 at 270 K were mainly in polyproline II helix conformations and in equilibrium between the alpha-helix and polyproline II helix, respectively, in harmony with the experiment. Interestingly, we found remarkably strong secondary-structure "templating"; namely, the alpha-helix template enhanced alpha-helix conformation and the beta-strand template induced beta-strand conformation in the simulated Ala8. The alpha-helix template lowered the nearby dielectric constant, which strengthened hydrogen bonds in the simulated Ala8, leading to alpha-helix stabilization. The beta-strand template provided hydrogen bond positions to the simulated Ala8, sharply inducing beta-strand structure. With or without templates, the energy landscape of Ala8 is always funnel-like and centered at the alpha-helix conformation, whereas entropic contribution disfavors the alpha-helix, leading to subtle competition. Secondary-structure templating may play a critical role in protein conformation dynamics in the cellular environment.

Conformational sampling with implicit solvent models: application to the PHF6 peptide in tau protein

Implicit solvent models approximate the effects of solvent through a potential of mean force and therefore make solvated simulations computationally efficient. Yet despite their computational efficiency, the inherent approximations made by implicit solvent models can sometimes lead to inaccurate results. To test the accuracy of a number of popular implicit solvent models, we determined whether implicit solvent simulations can reproduce the set of potential energy minima obtained from explicit solvent simulations. For these studies, we focus on a six-residue amino-acid sequence, referred to as the paired helical filament 6 (PHF6), which may play an important role in the formation of intracellular aggregates in patients with Alzheimer's disease. Several implicit solvent models form the basis of this work--two based on the generalized Born formalism, and one based on a Gaussian solvent-exclusion model. All three implicit solvent models generate minima that are in good agreement with minima obtained from simulations with explicit solvent. Moreover, free-energy profiles generated with each implicit solvent model agree with free-energy profiles obtained with explicit solvent. For the Gaussian solvent-exclusion model, we demonstrate that a straightforward ranking of the relative stability of each minimum suggests that the most stable structure is extended, a result in excellent agreement with the free-energy profiles. Overall, our data demonstrate that for some peptides like PHF6, implicit solvent can accurately reproduce the set of local energy minimum arising from quenched dynamics simulations with explicit solvent. More importantly, all solvent models predict that PHF6 forms extended beta-structures in solution, a finding consistent with the notion that PHF6 initiates neurofibrillary tangle formation in patients with Alzheimer's disease.

Implicit versus explicit solvent in free energy calculations of enzyme catalysis: Methyl transfer catalyzed by catechol O-methyltransferase

We compare free energy calculations for the methyl transfer reaction catalyzed by catechol O-methyltransferase using the quantum mechanical/molecular mechanical free energy method with implicit and explicit solvents. An analogous methylation reaction in a solution is also studied. For the explicit solvent model, we use the three-point transferable intermolecular potential model, and for the implicit model, we use the generalized Born molecular volume model as implemented in CHARMM. We find that activation and reaction free energies calculated with the two models are very similar, despite some structural differences that exist. A significant change in the polarization of the environment occurs as the reaction proceeds. This is more pronounced for the reaction in a solution than for the enzymatic reaction. For the enzymatic reaction, most of the changes take place in the protein rather than in the solvent, and, hence, the benefit of having an instantaneous relaxation of the solvent degrees of freedom is less pronounced for the enzymatic reaction than for the reaction in a solution. This is a likely reason why energies of the enzyme reaction are less sensitive to the choice of atomic radii than are energies of the reaction in a solution.

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.