Prediction of hydration free energies for aliphatic and aromatic chloro derivatives using molecular dynamics simulations with the OPLS-AA force field

All-atom molecular dynamics computer simulations were used to blindly predict the hydration free energies of a range of chloro-organic compounds as part of the SAMPL3 challenge. All compounds were parameterized within the framework of the OPLS-AA force field, using an established protocol to compute the absolute hydration free energy via a windowed free energy perturbation approach and thermodynamic integration. Three different approaches to deriving partial charge parameters were pursued: (1) using existing OPLS-AA atom types and charges with minor adjustments of partial charges on equivalent connecting atoms; (2) calculation of quantum mechanical charges via geometry optimization, followed by electrostatic potential (ESP) fitting, using Jaguar at the LMP2/cc-pVTZ(-F) level; and (3) via geometry optimization and CHelpG charges (Gaussian03 at the HF/6-31G* level), followed by two-stage RESP fitting. Protocol 3 generated the most accurate predictions with a root mean square (RMS) error of 1.2 kcal mol(-1) for the entire data set. It was found that the deficiency of the standard OPLS-AA parameters, protocol 1 (RMS error 2.4 kcal mol(-1) overall), was mostly due to compounds with more than three chlorine substituents on an aromatic ring. For this latter subset, the RMS errors were 1.4 kcal mol(-1) (protocol 3) and 4.3 kcal mol(-1) (protocol 1), respectively. We propose new OPLS-AA atom types for aromatic carbon and chlorine atoms in rings with ≥4 Cl-substituents that perform better than the best QM-based approach, resulting in an RMS error of 1.2 kcal mol(-1) for these difficult compounds.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Prediction of hydration free energies for the SAMPL4 diverse set of compounds using molecular dynamics simulations with the OPLS-AA force field

All-atom molecular dynamics computer simulations were used to blindly predict the hydration free energies of a range of small molecules as part of the SAMPL4 challenge. Compounds were parametrized on the basis of the OPLS-AA force field using three different protocols for deriving partial charges: (1) using existing OPLS-AA atom types and charges with minor adjustments of partial charges on equivalent connecting atoms and derivation of new parameters for a number of distinct chemical groups (N-alkyl imidazole, nitrate) that were not present in the published force field; (2) calculation of quantum mechanical charges via geometry optimization, followed by electrostatic potential (ESP) fitting, using Jaguar at the LMP2/cc-pVTZ(-F) level; and (3) via geometry optimization and CHelpG charges (Gaussian09 at the HF/6-31G* level), followed by two-stage RESP fitting. The absolute hydration free energy was computed by an established protocol including alchemical free energy perturbation with thermodynamic integration. The use of standard OPLS-AA charges (protocol 1) with a number of newly parametrized charges and the use of histidine derived parameters for imidazole yielded an overall root mean square deviation of the prediction from the experimental data of 1.75 kcal/mol. The precision of our results appears to be mainly limited by relatively poor reproducibility of the Lennard-Jones contribution towards the solvation free energy, for which we observed large variability that could be traced to a strong dependence on the initial system conditions.     

Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations

An efficient Monte Carlo (MC) algorithm including concerted rotations is directly compared to molecular dynamics (MD) in all-atom statistical mechanics folding simulations of small polypeptides. The previously reported algorithm "concerted rotations with flexible bond angles" (CRA) has been shown to successfully locate the native state of small polypeptides. In this study, the folding of three small polypeptides (trpzip2/H1/Trp-cage) is investigated using MC and MD, for a combined sampling time of approximately 10(11) MC configurations and 8 micros, respectively. Both methods successfully locate the experimentally determined native states of the three systems, but they do so at different speed, with 2-2.5 times faster folding of the MC runs. The comparison reveals that thermodynamic and dynamic properties can reliably be obtained by both and that results from folding simulations do not depend on the algorithm used. Similar to previous comparisons of MC and MD, it is found that one MD integration step of 2 fs corresponds to one MC scan, revealing the good sampling of MC. The simplicity and efficiency of the MC method will enable its future use in folding studies involving larger systems and the combination with replica exchange algorithms.     

Overcoming free energy barriers using unconstrained molecular dynamics simulations

Association of unconstrained molecular dynamics (MD) and the formalisms of thermodynamic integration and average force [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)] have been employed to determine potentials of mean force. When implemented in a general MD code, the additional computational effort, compared to other standard, unconstrained simulations, is marginal. The force acting along a chosen reaction coordinate xi is estimated from the individual forces exerted on the chemical system and accumulated as the simulation progresses. The estimated free energy derivative computed for small intervals of xi is canceled by an adaptive bias to overcome the barriers of the free energy landscape. Evolution of the system along the reaction coordinate is, thus, limited by its sole self-diffusion properties. The illustrative examples of the reversible unfolding of deca-L-alanine, the association of acetate and guanidinium ions in water, the dimerization of methane in water, and its transfer across the water liquid-vapor interface are examined to probe the efficiency of the method.     

Precision of free energies calculated by molecular dynamics simulations of peptides in solution

Errors in free energies for molecular replacement and for conformation change of a small model peptide have been determined empirically by repeated simulations from different starting points. All calculations have been done using thermodynamic integration, in which the system's potential energy is coupled to a parameter λ, that is increased or decreased by a small amount at each step of the simulation. The effects of several factors that may alter the precision are evaluated. These factors include: the length of the simulation, the dependence of the potential energy on λ, the use of conformational restraints, and their magnitude and form. The methods used for restraint and conformational forcing are described in detail. The free energy change, calculated as the mean from several successive simulations with alternately increasing and decreasing λ, is found to be independent of the length of the simulations. As expected, longer simulations produce more precise results. The variation of the calculated free energies is found to consist of two parts, a random error and a systematic hysteresis, i.e., a dependence on the direction in which λ changes. The hysteresis varies as the inverse of the length of the simulation and the random error as the inverse square root The advantage of the use of a different (nonlinear) dependence of the attractive and repulsive parts of the nonbonded potential energy on the coupling parameter when “creating” particles in solution is found to be very large. This nonlinear coupling was found to be superior to the use of linear coupling and a nonlinear change of the coupling parameter with the simulation time. The hysteresis in conformational free energy calculations is found to increase markedly if too weak a forcing restraint is chosen. It is shown how to deconvolute the contribution of a torsional restraint from the dependence of the free energy on a torsion angle.     

Predicting hydration free energies with chemical accuracy: the SAMPL4 challenge

An implicit solvent model described by a non-simple dielectric medium is used for the prediction of hydration free energies on the dataset of 47 molecules in the SAMPL4 challenge. The solute is represented by a minimal parameter set model based on a new all atom force-field, named the liquid simulation force-field. The importance of a first solvation shell correction to the hydration free energy prediction is discussed and two different approaches are introduced to address it: either with an empirical correction to a few functional groups (alcohol, ether, ester, amines and aromatic nitrogen), or an ab initio correction based on the formation of a solute/explicit water complex. Both approaches give equally good predictions with an average unsigned error <1 kcal/mol. Chemical accuracy is obtained.     

Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation

The raw ionic solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system shape, and size) and treatment of electrostatic interactions (Coulombic, lattice sum, or cutoff based) used during these simulations. In the present article, it is shown that correction terms can be derived for the effect of (A) an incorrect solvent polarization around the ion due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site and the possible presence of a liquid-vacuum interface in the simulated system. Taking the hydration free energy of the sodium cation as a test case, it is shown that the raw solvation free energies obtained using seven different types of boundary conditions and electrostatic schemes commonly used in explicit-solvent simulations (for a total of 72 simulations differing in the corresponding simulation parameters) can be corrected so as to obtain a consistent value for this quantity.     

Predicting hydration energies for multivalent ions

We have predicted the free energy of hydration for 40 monovalent and multivalent cations and anions using density functional theory and the implicit solvent model COnductor like Screening MOdel for Real Solvents (COSMO‐RS) at the Becke‐Perdew (BP)/Triple zeta valence with polarization functions (TZVP) level. Agreement with experimental data for monovalent and divalent ions is good and shows no significant systematic errors. Predictions are noticeably better than with standard COSMO. The agreement with experimental data for trivalent and tetravalent ions is slightly worse and shows systematic errors. Our results indicate that quantum chemical calculations combined with COSMO‐RS solvent treatment is a reliable method for treating multivalent ions in solution, provided one hydration shell of explicit water molecules is included for metal cations. The accuracy is not high enough to allow absolute predictions of hydration energies but could be used to investigate trends for several ions, thanks to the low computational cost, in particular for ligand exchange reactions. © 2014 Wiley Periodicals, Inc.     

Efficient calculation of many stacking and pairing free energies in DNA from a few molecular dynamics simulations

Through the use of the one-step perturbation approach, 130 free energies of base stacking and 1024 free energies of base pairing in DNA have been calculated from only five simulations of a nonphysical reference state. From analysis of a diverse set of 23 natural and unnatural bases, it appears that stacking free energies and stacking conformations play an important role in pairing of DNA nucleotides. On the one hand, favourable pairing free energies were found for bases that do not have the possibility to form canonical hydrogen bonds, while on the other hand, good hydrogen-bonding possibilities do not guarantee a favourable pairing free energy if the stacking of the bases dictates an unfavourable conformation. In this application, the one-step perturbation approach yields a wealth of both energetic and structural information at minimal computational cost.     

Absolute free energies of binding of peptide analogs to the HIV-1 protease from molecular dynamics simulations

The constants of binding of five peptide analogs to the active site of the HIV-1 aspartic-protease are calculated based on a novel sampling scheme that is efficient and does not introduce any approximations in addition to the energy function used to describe the system. The results agree with experiments. The squared correlation coefficient of the calculated vs. the measured values is 0.79. The sampling scheme consists of a series of molecular dynamics integrations with biases. The biases are selected based on an estimate of the probability density function of the system in a way to explore the conformational space and to reduce the statistical error in the calculated binding constants. The molecular dynamics integrations are done with a vacuum potential using a short cutoff scheme. To estimate the probability density of the simulated system, the results of the molecular dynamics integrations are combined using an extension of the weighted histogram analysis method (C. Bartels, Chem. Phys. Letters 331 (2000) 446-454). The probability density of the solvated ligand-protein system is obtained by applying a correction for the use of the short cutoffs in the simulations and by taking into account solvation with an electrostatic term and a hydrophobic term. The electrostatic part of the solvation is determined by finite difference Poisson-Boltzmann calculations; the hydrophobic part of the solvation is set proportional to the solvent accessible surface area. Setting the hydrophobic surface tension parameter equal to 8 mol(-1) K(-1) A(-2), absolute binding constants are in the muM to nM range. This is in agreement with experiments. The standard errors determined from eight repeated binding constant determinations are a factor of 14 to 411. A single determination of a binding constant is done with 499700 steps of molecular dynamics integration and 4500 finite difference Poisson-Boltzmann calculations. The simulations can be analyzed with respect to conformational changes of the active site of the HIV-1 protease or the ligands upon binding and provide information that complements experiments and can be used in the drug development process.     

Estimating the temperature dependence of peptide folding entropies and free enthalpies from total energies in molecular dynamics simulations

The temperature dependence of thermodynamic quantities, such as heat capacity, entropy and free enthalpy, may be obtained by using well-known equations that relate these quantities to the enthalpy of the molecular system of interest at a range of temperatures. In turn, the enthalpy of a molecular system can be estimated from molecular dynamics simulations of an appropriate model. To demonstrate this, we have investigated the temperature dependence of the enthalpy, heat capacity, entropy and free enthalpy of a system that consists of a beta-heptapeptide in methanol and have used the statistical mechanics relationships to describe the thermodynamics of the folding/unfolding equilibrium of the peptide. The results illustrate the power of current molecular simulation force fields and techniques in establishing the link between thermodynamic quantities and conformational distributions.     

Calculation of the entropy and free energy of peptides by molecular dynamics simulations using the hypothetical scanning molecular dynamics method

Hypothetical scanning (HS) is a method for calculating the absolute entropy S and free energy F from a sample generated by any simulation technique. With this approach each sample configuration is reconstructed with the help of transition probabilities (TPs) and their product leads to the configuration's probability, hence to the entropy. Recently a new way for calculating the TPs by Monte Carlo (MC) simulations has been suggested, where all system interactions are taken into account. Therefore, this method--called HSMC--is in principle exact where the only approximation is due to insufficient sampling. HSMC has been applied very successfully to liquid argon, TIP3P water, self-avoiding walks on a lattice, and peptides. Because molecular dynamics (MD) is considered to be significantly more efficient than MC for a compact polymer chain, in this paper HSMC is extended to MD simulations as applied to peptides. Like before, we study decaglycine in vacuum but for the first time also a peptide with side chains, (Val)(2)(Gly)(6)(Val)(2). The transition from MC to MD requires implementing essential changes in the reconstruction process of HSMD. Results are calculated for three microstates, helix, extended, and hairpin. HSMD leads to very stable differences in entropy TDeltaS between these microstates with small errors of 0.1-0.2 kcal/mol (T=100 K) for a wide range of calculation parameters with extremely high efficiency. Various aspects of HSMD and plans for future work are discussed.     

Revealing binding selectivity of inhibitors toward bromodomain-containing proteins 2 and 4 using multiple short molecular dynamics simulations and free energy analyses

ABSTRACT

Emerging evidences indicate bromodomain-containing proteins 2 and 4 (BRD2 and BRD4) play critical roles in cancers, inflammations, cardiovascular diseases and other pathologies. Multiple short molecular dynamics (MSMD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method were applied to investigate the binding selectivity of three inhibitors 87D, 88M and 89G towards BRD2 over BRD4. The root-mean-square fluctuation (RMSF) analysis indicates that the structural flexibility of BRD4 is stronger than that of BRD2. Moreover the calculated distances between the C_α atoms in the centres of the ZA_loop and BC_loop of BRD4 are also bigger than that of BRD2. The rank of binding free energies calculated using MM-GBSA method agrees well with that determined by experimental data. The results show that 87D can bind more favourably to BRD2 than BRD4, while 88M has better selectivity on BRD4 over BRD2. Residue-based free-energy decomposition method was utilized to estimate the inhibitor-residue interaction spectrum and the results not only identify the hot interaction spots of inhibitors with BRD2 and BRD4, but also demonstrate that several common residues, including (W370, W374), (P371, P375), (V376, V380) and (L381, L385) belonging to (BRD2, BRD4), generate significant binding difference of inhibitors to BRD2 and BRD4.     

Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, lambda, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter lambda is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the lambda free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (lambda = 0 and lambda = 1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature T(lambda) is assigned to lambda and a fictitious mass m(lambda) is introduced as a means of creating an adiabatic separation between lambda and the rest of the system. Under these conditions, it is shown that the lambda free energy profile can be directly computed from the adiabatic probability distribution function of lambda without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new lambda adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories.     

Comparison between self-guided Langevin dynamics and molecular dynamics simulations for structure refinement of protein loop conformations

This article presents a comparative analysis of two replica‐exchange simulation methods for the structure refinement of protein loop conformations, starting from low‐resolution predictions. The methods are self‐guided Langevin dynamics (SGLD) and molecular dynamics (MD) with a Nosé–Hoover thermostat. We investigated a small dataset of 8‐ and 12‐residue loops, with the shorter loops placed initially from a coarse‐grained lattice model and the longer loops from an enumeration assembly method (the Loopy program). The CHARMM22 + CMAP force field with a generalized Born implicit solvent model (molecular‐surface parameterized GBSW2) was used to explore conformational space. We also assessed two empirical scoring methods to detect nativelike conformations from decoys: the all‐atom distance‐scaled ideal‐gas reference state (DFIRE‐AA) statistical potential and the Rosetta energy function. Among the eight‐residue loop targets, SGLD out performed MD in all cases, with a median of 0.48 Å reduction in global root‐mean‐square deviation (RMSD) of the loop backbone coordinates from the native structure. Among the more challenging 12‐residue loop targets, SGLD improved the prediction accuracy over MD by a median of 1.31 Å, representing a substantial improvement. The overall median RMSD for SGLD simulations of 12‐residue loops was 0.91 Å, yielding refinement of a median 2.70 Å from initial loop placement. Results from DFIRE‐AA and the Rosetta model applied to rescoring conformations failed to improve the overall detection calculated from the CHARMM force field. We illustrate the advantage of SGLD over the MD simulation model by presenting potential‐energy landscapes for several loop predictions. Our results demonstrate that SGLD significantly outperforms traditional MD in the generation and populating of nativelike loop conformations and that the CHARMM force field performs comparably to other empirical force fields in identifying these conformations from the resulting ensembles. Published 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Solvation free energies and hydration structure of N-methyl-p-nitroaniline

Solvation Gibbs energies of N-methyl-p-nitroaniline (MNA) in water and 1-octanol are calculated using the expanded ensemble molecular dynamics method with a force field taken from the literature. The accuracy of the free energy calculations is verified with the experimental Gibbs free energy data and found to reproduce the experimental 1-octanol∕water partition coefficient to within ±0.1 in log unit. To investigate the hydration structure around N-methyl-p-nitroaniline, an independent NVT molecular dynamics simulation was performed at ambient conditions. The local organization of water molecules around the solute MNA molecule was investigated using the radial distribution function (RDF), the coordination number, and the extent of hydrogen bonding. The spatial distribution functions (SDFs) show that the water molecules are distributed above and below the nitrogen atoms parallel to the plane of aromatic ring for both the methylamino and nitro functional groups. It is found that these groups have a significant effect on the hydration of MNA with water molecules forming two weak hydrogen bonds with both the methylamino and nitro groups. The hydration structures around the functional groups in MNA in water are different from those that have been found for methylamine, nitrobenzene, and benzene in aqueous solutions, and these differences together with weak hydrogen bonds explain the lower solubility of MNA in water. The RDFs together with SDFs provide a tool for the understanding the hydration of MNA (and other molecules) and therefore their solubility.     

Predicting hydration free energies of neutral compounds by a parametrization of the polarizable continuum model

A parametrization of the polarizable continuum model (PCM) is presented having the experimental hydration free energies of 215 neutral molecules as target. The cavitation and dispersion contributions were based on the Tu?on-Silla-Pascual-Ahuir (Tu?on; et al. Chem. Phys. Lett. 1993, 203, 289) and Floris-Tomasi (Floris, F.; Tomasi, J. J. Comput. Chem. 1989, 10, 616) expressions, respectively. Both the polar and nonpolar contributions were evaluated on the same solvent-excluding molecular surface that used unscaled Bondi atomic radii. The parametrization was provided for the HF, Xalpha, LSDA, B3LYP, and mPW1PW91 methods at the 6-31G(d) basis set, and the results are in fair agreement with the experimental data. For the sake of comparison, the PCM(UAHF) and our parametrization (PCM2), both at HF level, have produced DeltaG(PCM(UAHF)) = aDeltaGexp (a = 1.02 +/- 0.02, r = 0.945, sd = 0.987, Ftest = 1778) and DeltaG(PCM2) = aDeltaGexp (a = 0.95 +/- 0.02, r = 0.952, sd = 0.843, Ftest = 2070), respectively. The mean absolute deviations from experimental data were 0.67 and 0.68 kcal/mol for PCM(UAHF) and PCM2, respectively.     

Determination of solvation free energies by adaptive expanded ensemble molecular dynamics

A new method of calculating absolute free energies is presented. It was developed as an extension to the expanded ensemble molecular dynamics scheme and uses probability density estimation to continuously optimize the expanded ensemble parameters. The new method is much faster as it removes the time-consuming and expertise-requiring step of determining balancing factors. Its efficiency and accuracy are demonstrated for the dissolution of three qualitatively very different chemical species in water: methane, ionic salts, and benzylamine. A recently suggested optimization scheme by Wang and Landau [Phys. Rev. Lett. 86, 2050 (2001)] was also implemented and found to be computationally less efficient than the proposed adaptive expanded ensemble method.