Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smooth-permittivity finite difference Poisson-Boltzmann method

A fast stable finite difference Poisson-Boltzmann (FDPB) model for implicit solvation in molecular dynamics simulations was developed using the smooth permittivity FDPB method implemented in the OpenEye ZAP libraries. This was interfaced with two widely used molecular dynamics packages, AMBER and CHARMM. Using the CHARMM-ZAP software combination, the implicit solvent model was tested on eight proteins differing in size, structure, and cofactors: calmodulin, horseradish peroxidase (with and without substrate analogue bound), lipid carrier protein, flavodoxin, ubiquitin, cytochrome c, and a de novo designed 3-helix bundle. The stability and accuracy of the implicit solvent simulations was assessed by examining root-mean-squared deviations from crystal structure. This measure was compared with that of a standard explicit water solvent model. In addition we compared experimental and calculated NMR order parameters to obtain a residue level assessment of the accuracy of MD-ZAP for simulating dynamic quantities. Overall, the agreement of the implicit solvent model with experiment was as good as that of explicit water simulations. The implicit solvent method was up to eight times faster than the explicit water simulations, and approximately four times slower than a vacuum simulation (i.e., with no solvent treatment).     

Comparison of various implicit solvent models in molecular dynamics simulations of immunoglobulin G light chain dimer

The present study tests performance of different solvation models applied to molecular dynamics simulation of a large, dimeric protein molecule. Analytical Continuum Electrostatics (ACE) with two different parameter sets, older V98 and new V01, and Effective Energy Function (EEF) are employed in molecular dynamics simulation of immunoglobulin G (IgG) light chain dimer and variable domain of IgG light chain. Results are compared with explicit solvent and distance dependent dielectric constant (DDE) calculations. The overall analysis shows that the EEF method yields results comparable to explicit solvent simulations; however, the stability of simulations is lower. On the other hand, the ACE_V98 model does not seem to achieve the accuracy or stability expected in nanosecond timescale MD simulation for the studied systems. The ACE_V01 model greatly improves stability of the calculation; nonetheless, changes in radius of gyration and solvent accessible surface of the studied systems may indicate that the parameter set still needs to be improved if the method is supposed to be used for simulations of large, polymeric proteins. Additionally, electrostatic contribution to the solvation free energy calculated in the ACE model is compared with a numerical treatment of the dielectric continuum model. Wall clock time of all simulations is compared. It shows that EEF calculation is six times faster than corresponding ACE and 50 times faster than explicit solvent simulations.     

Efficient approximate all-atom solvent accessible surface area method parameterized for folded and denatured protein conformations

Continuing advances in computer hardware and software are permitting atomic-resolution molecular simulations for longer time scales and on larger systems. Despite these advances, routinely performing atomistic simulations with explicit water for even small proteins, which reach the folding time of such proteins, remains intractable for the foreseeable future. An implicit approximation of the solvent environment using a solvent accessible surface area (SASA) term in a molecular mechanics potential function allows exclusion of the explicit water molecules in protein simulations. This reduces the number of particles by approximately an order of magnitude. We present a fast and acceptably accurate approximate all-atom SASA method parameterized using a set of folded and heat-denatured conformations of globular proteins. The parameters are shown to be transferable to folded and heat-denatured conformations for another set of proteins. Calculation of the approximate SASA and the associated derivatives with respect to atomic positions for a 4644 atom protein requires only 1/11th the CPU time required for calculation of the nonbonded interactions for this system. On a per atom basis, this algorithm is three times faster than the fastest previously published approximate SASA method and achieves the same level of accuracy.     

IMiCMO: a new integrated ab initio multicenter molecular orbitals method for molecular dynamics calculations in solvent cluster systems

A new computational scheme integrating ab initio multicenter molecular orbitals for determining forces of individual atoms in large cluster systems is presented. This method can be used to treat systems of many molecules, such as solvents by quantum mechanics. The geometry parameters obtained for three models of water clusters by the present method are compared with those obtained by the full ab initio MO method. The results agree very well. The scheme proposed in this article also intended for use in modeling cluster systems using parallel algorithms. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1107–1112, 2001     

Molecular dynamics simulation of the hydrocarbon region of a biomembrane using a reduced representation model

The development of a coarse‐grained reduced‐representation model of the hydrocarbon region of a biological membrane is reported. The potential is based on the popular Gay–Berne model of liquid crystals, and involves the linking of individual Gay–Berne ellipsoids by harmonic springs to form each hydrocarbon chain. Diffusion coefficients and order parameters have been calculated by molecular dynamics computer simulations for a range of parameter sets. The results clearly demonstrate the presence of a phase transition from an ordered low‐temperature solid phase reminiscent of the L_β′ phase of phospholipids, to a high‐temperature disordered phase reminiscent of the L_α phase. Order parameters calculated for each layer of the model are consistent with the experimental segmental order parameters reported for dipalmitoyl phosphatidylcholine. The application of this model to the study of small molecule diffusion within the membrane core is proposed. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1622–1633, 2001     

Implicit solvent models: Combining an analytical formulation of continuum electrostatics with simple models of the hydrophobic effect

The recent development of approximate analytical formulations of continuum electrostatics opens the possibility of efficient and accurate implicit solvent models for biomolecular simulations. One such formulation (ACE, Schaefer & Karplus, J. Phys. Chem., 1996, 100:1578) is used to compute the electrostatic contribution to solvation and conformational free energies of a set of small solutes and three proteins. Results are compared to finite-difference solutions of the Poisson equation (FDPB) and explicit solvent simulations and experimental data where available. Small molecule solvation free energies agree with FDPB within 1–1.5 kcal/mol, which is comparable to differences in FDPB due to different surface treatments or different force field parameterizations. Side chain conformation free energies of aspartate and asparagine are in qualitative agreement with explicit solvent simulations, while 74 conformations of a surface loop in the protein Ras are accurately ranked compared to FDPB. Preliminary results for solvation free energies of small alkane and polar solutes suggest that a recent Gaussian model could be used in combination with analytical continuum electrostatics to treat nonpolar interactions. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 322–335, 1999     

A multiscale model for kinetics of formation and disintegration of spherical micelles

Dynamics of self-assembly and structural transitions in surfactant systems often involve a large span of length and time scales. A comprehensive understanding of these processes requires development of models connecting phenomena taking place on different scales. In this paper, we develop a multiscale model for formation and disintegration of spherical nonionic micelles. The study is performed under the assumption that the dominant mechanism of micelle formation (disintegration) is a stepwise addition (removal) of single monomers to (from) a surfactant aggregate. Different scales of these processes are investigated using a combination of coarse-grained molecular dynamics simulations, analytical and numerical solution of stochastic differential equations, and a numerical solution of kinetic equations. The removal of a surfactant from an aggregate is modeled by a Langevin equation for a single reaction coordinate, the distance between the centers of mass of the surfactant and the aggregate, with parameters obtained from a series of constrained molecular dynamics simulations. We demonstrate that the reverse process of addition of a surfactant molecule to an aggregate involves at least two additional degrees of freedom, orientation of the surfactant molecule and micellar microstructure. These additional degrees of freedom play an active role in the monomer addition process and neglecting their contribution leads to qualitative discrepancies in predicted surfactant addition rates. We propose a stochastic model for the monomer addition which takes the two additional degrees of freedom into account and extracts the model parameters from molecular dynamics simulations. The surfactant addition rates are determined from Brownian dynamics simulations of this model. The obtained addition and removal rates are then incorporated into the kinetic model of micellar formation and disintegration.     

Incorporation of solvent effects into density functional predictions of molecular polarizabilities and hyperpolarizabilities

A study of the effect of the field, the basis set, the functional, and the cavity size on molecular polarizabilities and hyperpolarizabilities of substituted benzenes in liquid or solution is reported. The calculations have been performed using the density functional theory (DFT) within the conductor‐like screening model (COSMO). The optimized computational parameters are adopted to calculate molecular polarizabilities and hyperpolarizabilities of substituted benzenes in liquid or solution. The results show good agreement with the experimental values. From comparison of the different theoretical results, it is found that at the same theoretical level, the selection of the different solvation models may play an important role in the calculations of molecular solvation polarizability, and using the same solvation model, the effects of the different theoretical methods are relatively small. © 2001 John Wiley & Sons, Inc. Int J Quantum Chem, 2001     

A semi-implicit solvent model for the simulation of peptides and proteins

We present a new model of biomolecules hydration based on macroscopic electrostatic theory, that can both describe the microscopic details of solvent-solute interactions and allow for an efficient evaluation of the electrostatic hydration free energy. This semi-implicit model considers the solvent as an ensemble of polarizable pseudoparticles whose induced dipole describe both the electronic and orientational solvent polarization. In the presented version of the model, there is no mutual dipolar interaction between the particles, and they only interact through short-ranged Lennard-Jones interactions. The model has been integrated into a molecular dynamics code, and offers the possibility to simulate efficiently the conformational evolution of biomolecules. It is able to provide estimations of the electrostatic solvation free energy within short time windows during the simulation. It has been applied to the study of two small peptides, the octaalanine and the N-terminal helix of ribonuclease A, and two proteins, the bovine pancreatic trypsin inhibitor and the B1 immunoglobin-binding domain of streptococcal protein G. Molecular dynamics simulations of these biomolecules, using a slightly modified Amber force field, provide stable and meaningful trajectories in overall agreement with experiments and all-atom simulations. Correlations with respect to Poisson-Boltzmann electrostatic solvation free energies are also presented to discuss the parameterization of the model and its consequences.     

FACTS: Fast analytical continuum treatment of solvation

An efficient method for calculating the free energy of solvation of a (macro)molecule embedded in a continuum solvent is presented. It is based on the fully analytical evaluation of the volume and spatial symmetry of the solvent that is displaced from around a solute atom by its neighboring atoms. The two measures of solvent displacement are combined in empirical equations to approximate the atomic (or self) electrostatic solvation energy and the solvent accessible surface area. The former directly yields the effective Born radius, which is used in the generalized Born (GB) formula to calculate the solvent-screened electrostatic interaction energy. A comparison with finite-difference Poisson data shows that atomic solvation energies, pair interaction energies, and their sums are evaluated with a precision comparable to the most accurate GB implementations. Furthermore, solvation energies of a large set of protein conformations have an error of only 1.5%. The solvent accessible surface area is used to approximate the nonpolar contribution to solvation. The empirical approach, called FACTS (Fast Analytical Continuum Treatment of Solvation), is only four times slower than using the vacuum energy in molecular dynamics simulations of proteins. Notably, the folded state of structured peptides and proteins is stable at room temperature in 100-ns molecular dynamics simulations using FACTS and the CHARMM force field.     

Extracting parameters for base-pair level models of DNA from molecular dynamics simulations

 A method is described to extract a complete set of sequence-dependent energy parameters for a rigid base-pair model of DNA from molecular dynamics (MD) simulations. The method is properly consistent with equilibrium statistical mechanics and leads to effective inertia parameters for the base-pair units as well as stacking and stiffness parameters for the base-pair junctions. We give explicit formulas that yield a complete set of base-pair model parameters in terms of equilibrium averages that can be estimated from a time series generated in an MD simulation. The expressions to be averaged depend strongly both on the choice of coordinates used to describe rigid-body orientations and on the choice of strain measures at each junction. Received: 12 July 2000 / Accepted: 5 January 2001 / Published online: 3 May 2001     

The use of a generalized born model for the analysis of protein conformational transitions: a comparative study with explicit solvent simulations for chemotaxis Y protein (CheY)

To investigate whether implicit solvent models are appropriate for mechanistic studies of conformational transition in proteins, a recently developed generalized Born model (GBSW) was applied to a small signaling protein, chemotaxis protein Y (CheY), with different combinations of the phosphorylation state and conformation of the system; the results were compared to explicit solvent simulations using a stochastic boundary condition. The subtle but distinct conformational transitions involved in CheY activation makes the system ideally suited for comparing implicit and explicit solvent models because these conformational transitions are potentially accessible in both types of simulations. The structural and dynamical properties analyzed include not only those localized to the active site region but also throughout the protein, such as sidechain methyl group order parameters, backbone hydrogen bonding lifetime and occupancy as well as principal components of the trajectories. Overall, many properties were well reproduced by the GBSW simulations when compared with the explicit solvent calculations, although a number of observations consistently point to the suggestion that the current parameterization of the GBSW model tends to overestimate hydrogen-bonding interactions involving both charged groups and (charge-neutral) backbone atoms. This deficiency led to overstabilization of certain secondary structural motifs and more importantly, qualitatively different behaviors for the active site groups (Thr 87, Ala 88, the beta4-alpha4 loop) in response to phosphorylation, when compared with explicit solvent simulations. The current study highlights the value of carrying out both explicit and implicit solvent simulations for complementary mechanistic insights in the analysis of conformational transition in biomolecules.     

Surveying implicit solvent models for estimating small molecule absolute hydration free energies

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1-BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0-1.7 kcal/mol and R(2)=0.81-0.91) and with experimental hydration free energies (average unsigned errors=1.1-1.4 kcal/mol and R(2)=0.66-0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.     

The SGB/NP hydration free energy model based on the surface generalized born solvent reaction field and novel nonpolar hydration free energy estimators

The development and parameterization of a solvent potential of mean force designed to reproduce the hydration thermodynamics of small molecules and macromolecules aimed toward applications in conformation prediction and ligand binding free energy prediction is presented. The model, named SGB/NP, is based on a parameterization of the Surface Generalized Born continuum dielectric electrostatic model using explicit solvent free energy perturbation calculations and a newly developed nonpolar hydration free energy estimator motivated by the results of explicit solvent simulations of the thermodynamics of hydration of hydrocarbons. The nonpolar model contains, in addition to the more commonly used solvent accessible surface area term, a component corresponding to the attractive solute-solvent interactions. This term is found to be important to improve the accuracy of the model, particularly for cyclic and hydrogen bonding compounds. The model is parameterized against the experimental hydration free energies of a set of small organic molecules. The model reproduces the experimental hydration free energies of small organic molecules with an accuracy comparable or superior to similar models employing more computationally demanding estimators and/or a more extensive set of parameters.     

Strategies to model the near-solute solvent molecular density/polarization

The solvent molecular distribution significantly affects the behavior of the solute molecules and is thus important in studying many biological phenomena. It can be described by the solvent molecular density distribution, g, and the solvent electric dipole distribution, p. The g and p can be computed directly by counting the number of solvent molecules/dipoles in a microscopic volume centered at r during a simulation or indirectly from the mean force F and electrostatic field E acting on the solvent molecule at r, respectively. However, it is not clear how the g and p derived from simulations depend on the solvent molecular center or the solute charge and if the g(F) and p(E) computed from the mean force and electric field acting on the solvent molecule, respectively, could reproduce the corresponding g and p obtained by direct counting. Hence, we have computed g, p, g(F), and p(E) using different water centers from simulations of a solute atom of varying charge solvated in TIP3P water. The results show that g(F) and p(E) can reproduce the g and p obtained using a given count center. This implies that rather than solving the coordinates of each water molecule by MD simulations, the distribution of water molecules could be indirectly obtained from analytical formulas for the mean force F and electrostatic field E acting on the solvent molecule at r. Furthermore, the dependence of the g and p distributions on the solute charge revealed provides an estimate of the change in g and p surrounding a biomolecule upon a change in its conformation.     

Anomalous mixing behavior of polyisobutylene/polypropylene blends: molecular dynamics simulation study

The unusual mixing behavior of polyisobutylene (PIB) with head-to-head (hhPP) and head-to-tail polypropylene (PP) is studied using large-scale molecular dynamics (MD). The heats of mixing and Flory chi parameters were computed from MD simulations of both blends using a united atom model. The chi parameters from the simulations were estimated from the structure factors using the random phase approximation in analogy with neutron scattering (SANS) experiments. MD simulations for syndiotactic hhPP/PIB predicted a lower critical solution temperature with a chi parameter in very good agreement with SANS experiments on the atactic hhPP/PIB blend. MD simulations also predicted that the isotactic PP/PIB blend was immiscible at high molecular weight in qualitative agreement with cloud point measurements on atactic PP/PIB.     

Stability of ion triplets in ionic liquid/lithium salt solutions: Insights from implicit and explicit solvent models and molecular dynamics simulations

Binding energies of ion triplets formed in ionic liquids by Li⁺ with two anions have been studied using quantum‐chemical calculations with implicit and explicit solvent supplemented by molecular dynamics (MD) simulations. Explicit solvent approach confirms variation of solute‐ionic liquid interactions at distances up to 2 nm, resulting from structure of solvation shells induced by electric field of the solute. Binding energies computed in explicit solvent and from the polarizable continuum model approach differ largely, even in sign, but relative values generally agree between these two models. Stabilities of ion triplets obtained in quantum‐chemical calculations for some systems disagree with MD results; the discrepancy is attributed to the difference between static optimized geometries used in quantum chemical modeling and dynamic structures of triplets in MD simulations. © 2015 Wiley Periodicals, Inc.     

Structure-activity relationships of pyrethroid insecticides. Part 2. The use of molecular dynamics for conformation searching and average parameter calculation

Summary Molecular dynamics simulations have been performed on a number of conformationally flexible pyrethroid insecticides. The results indicate that molecular dynamics is a suitable tool for conformational searching of small molecules given suitable simulation parameters. The structures derived from the simulations are compared with the static conformation used in a previous study. Various physicochemical parameters have been calculated for a set of conformations selected from the simulations using multivariate analysis. The averaged values of the parameters over the selected set (and the factors derived from them) are compared with the single conformation values used in the previous study.     

Elucidating the vibrational spectra of hydrogen-bonded aggregates in solution: electronic structure calculations with implicit solvent and first-principles molecular dynamics simulations with explicit solvent for 1-hexanol in n-hexane

Fourier transform infrared spectroscopy is a popular method for the experimental investigation of hydrogen-bonded aggregates, but linking spectral information to microscopic information on aggregate size distribution and aggregate architecture is an arduous task. Static electronic structure calculations with an implicit solvent model, Car-Parrinello molecular dynamics (CPMD) using the Becke-Lee-Yang-Parr (BLYP) exchange and correlation energy functionals and classical molecular dynamics simulations for the all-atom version of the optimized parameters for liquid simulations (OPLS-AA) force field were carried out for an ensemble of 1-hexanol aggregates solvated in n-hexane. The initial configurations for these calculations were size-selected from a distribution of aggregates obtained from a large-scale Monte Carlo simulation. The vibrational spectra computed from the static electronic structure calculations for monomers and dimers and from the CPMD simulations for aggregates up to pentamers demonstrate the extent of the contribution of dangling or nondonating hydroxyl groups found in linear and branched aggregates to the "monomeric" peak. Furthermore, the computed spectra show that there is no simple relationship between peak shift and aggregate size nor architecture, but the effect of hydrogen-bond cooperativity is shown to differentiate polymer-like (cooperative) and dimer-like (noncooperative) hydrogen bonds in the vibrational spectrum. In contrast to the static electronic structure calculations and the CPMD simulations, the classical molecular dynamics simulations greatly underestimate the vibrational peak shift due to hydrogen bonding.     

MOCASSIN: Metropolis and kinetic Monte Carlo for solid electrolytes

The combination of density functional theory and Monte Carlo simulations is a powerful approach for the atomistic modeling of defect transport in solid electrolytes. The present contribution introduces the MOCASSIN software (Monte Carlo for Solid State Ionics) for kinetic and Metropolis Monte Carlo simulations of crystalline materials. MOCASSIN combines model building, visualization, and simulation, aiming to provide accessible MC for end users. Developed for the investigation of solid electrolytes, MOCASSIN is ideal for screening common variation parameters, such as temperature and doping fraction. The input effort is minimized using space groups for processing symmetry. The graphical interface for model building allows complex model input, including multiple mobile species, multiple migration paths, small polaron hopping, vehicle movements, multiple complex migration mechanisms, and custom interaction clusters. The software is provided free of charge for noncommercial usage.