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.     

Accuracy of free energies of hydration using CM1 and CM3 atomic charges

Absolute free energies of hydration (DeltaGhyd) have been computed for 25 diverse organic molecules using partial atomic charges derived from AM1 and PM3 wave functions via the CM1 and CM3 procedures of Cramer, Truhlar, and coworkers. Comparisons are made with results using charges fit to the electrostatic potential surface (EPS) from ab initio 6-31G* wave functions and from the OPLS-AA force field. OPLS Lennard-Jones parameters for the organic molecules were used together with the TIP4P water model in Monte Carlo simulations with free energy perturbation theory. Absolute free energies of hydration were computed for OPLS united-atom and all-atom methane by annihilating the solutes in water and in the gas phase, and absolute DeltaGhyd values for all other molecules were computed via transformation to one of these references. Optimal charge scaling factors were determined by minimizing the unsigned average error between experimental and calculated hydration free energies. The PM3-based charge models do not lead to lower average errors than obtained with the EPS charges for the subset of 13 molecules in the original study. However, improvement is obtained by scaling the CM1A partial charges by 1.14 and the CM3A charges by 1.15, which leads to average errors of 1.0 and 1.1 kcal/mol for the full set of 25 molecules. The scaled CM1A charges also yield the best results for the hydration of amides including the E/Z free-energy difference for N-methylacetamide in water.     

Modifying the OPLS-AA force field to improve hydration free energies for several amino acid side chains using new atomic charges and an off-plane charge model for aromatic residues

The hydration free energies of amino acid side chains are an important determinant of processes that involve partitioning between different environments, including protein folding, protein complex formation, and protein-membrane interactions. Several recent papers have shown that calculated hydration free energies for polar and aromatic residues (Trp, His, Tyr, Asn, Gln, Asp, Glu) in several common molecular dynamics force fields differ significantly from experimentally measured values. We have attempted to improve the hydration energies for these residues by modifying the partial charges of the OPLS-AA force field based on natural population analysis of density functional theory calculations. The resulting differences between calculated hydration free energies and experimental results for the seven side chain analogs are less than 0.1 kcal/mol. Simulations of the synthetic Trp-rich peptide Trpzip2 show that the new charges lead to significantly improved geometries for interacting Trp-side chains. We also investigated an off-plane charge model for aromatic rings that more closely mimics their electronic configuration. This model results in an improved free energy of hydration for Trp and a somewhat altered benzene-sodium potential of mean force with a more favorable energy for direct benzene-sodium contact.     

Predicting hydration Gibbs energies of alkyl-aromatics using molecular simulation: a comparison of current force fields and the development of a new parameter set for accurate solvation data

The Gibbs energy of hydration is an important quantity to understand the molecular behavior in aqueous systems at constant temperature and pressure. In this work we review the performance of some popular force fields, namely TraPPE, OPLS-AA and Gromos, in reproducing the experimental Gibbs energies of hydration of several alkyl-aromatic compounds--benzene, mono-, di- and tri-substituted alkylbenzenes--using molecular simulation techniques. In the second part of the paper, we report a new model that is able to improve such hydration energy predictions, based on Lennard Jones parameters from the recent TraPPE-EH force field and atomic partial charges obtained from natural population analysis of density functional theory calculations. We apply a scaling factor determined by fitting the experimental hydration energy of only two solutes, and then present a simple rule to generate atomic partial charges for different substituted alkyl-aromatics. This rule has the added advantages of eliminating the unnecessary assumption of fixed charge on every substituted carbon atom and providing a simple guideline for extrapolating the charge assignment to any multi-substituted alkyl-aromatic molecule. The point charges derived here yield excellent predictions of experimental Gibbs energies of hydration, with an overall absolute average deviation of less than 0.6 kJ mol(-1). This new parameter set can also give good predictive performance for other thermodynamic properties and liquid structural information.     

Solvation properties of N-acetyl-β-glucosamine: molecular dynamics study incorporating electrostatic polarization

N-Acetyl-β-glucosamine (NAG) is an important moiety of glycoproteins and is involved in many biological functions. However, conformational and dynamical properties of NAG molecules in aqueous solution, the most common biological environment, remain ambiguous due to limitations of experimental methods. Increasing efforts are made to probe structural properties of NAG and NAG-containing macromolecules, like peptidoglycans and polymeric chitin, at the atomic level using molecular dynamics simulations. In this work, we develop a polarizable carbohydrate force field for NAG and contrast simulation results of various properties using this novel force field and an analogous nonpolarizable (fixed charge) model. Aqueous solutions of NAG and its oligomers are investigated; we explore conformational properties (rotatable bond geometry), electrostatic properties (dipole moment distribution), dynamical properties (self-diffusion coefficient), hydrogen bonding (water bridge structure and dynamics), and free energy of hydration. The fixed-charge carbohydrate force field exhibits deviations from the gas phase relative rotation energy of exocyclic hydroxymethyl side chain and of chair/boat ring distortion. The polarizable force field predicts conformational properties in agreement with corresponding first-principles results. NAG-water hydrogen bonding pattern is studied through radial distribution functions (RDFs) and correlation functions. Intermolecular hydrogen bonding between solute and solvent is found to stabilize NAG solution structures while intramolecular hydrogen bonds define glycosidic linkage geometry of NAG oligomers. The electrostatic component of hydration free energy is highly dependent on force field atomic partial charges, influencing a more favorable free energy of hydration in the fixed-charge model compared to the polarizable model.     

Detailed considerations for a balanced and broadly applicable force field: a study of substituted benzenes modeled with OPLS-AA

Modern classical force fields have been traditionally parameterized by attempting to maximize agreement to any number of experimental and/or quantum mechanical target properties. As these force fields are pushed towards obtaining quantitative estimates of often subtle energetic differences, stringent and consistent parameterization criteria, particularly in regard to charge distributions, are required to ensure that systematic errors cancel, that parameters are transferable between molecules, and that performance does not significantly deteriorate when using more approximate methods, such as with continuum solvent models. Relative free energies of hydration are presented here for 40 mono- and disubstituted benzenes modeled with the OPLS-AA force field; heats of vaporization and pure liquid densities at standard conditions are presented when experimental data is available. Overall agreement between OPLS-AA and experiment is remarkable (average error = 0.5 kcal/mol for DeltaDeltaG(hydration), 1.0 kcal/mol for DeltaH(vap) (0), 0.02 g/mL for densities), yet several functional groups are identified as having consistent and correctable errors (alkyl-, nitro-, and thiobenzenes). Relative free energies of hydration obtained with rigorous free energy perturbations using explicit solvent are also compared with energies from minimizations using a generalized Born model (GB). There is high correlation between these estimates (R = 0.99), and as demonstrated here, reparameterization of the aforementioned groups can be guided with rapid GB calculations.     

Comparison of two simulation methods to compute solvation free energies and partition coefficients

The thermodynamic integration (TI) and expanded ensemble (EE) methods are used here to calculate the hydration free energy in water, the solvation free energy in 1‐octanol, and the octanol‐water partition coefficient for a six compounds of varying functionality using the optimized potentials for liquid simulations (OPLS) all‐atom (AA) force field parameters and atomic charges. Both methods use the molecular dynamics algorithm as a primary component of the simulation protocol, and both have found wide applications in fields such as the calculation of activity coefficients, phase behavior, and partition coefficients. Both methods result in solvation free energies and 1‐octanol/water partition coefficients with average absolute deviations (AAD) from experimental data to within 4 kJ/mol and 0.5 log units, respectively. Here, we find that in simulations the OPLS‐AA force field parameters (with fixed charges) can reproduce solvation free energies of solutes in 1‐octanol with AAD of about half that for the solute hydration free energies using a extended simple point charge (SPC/E) model of water. The computational efficiency of the two simulation methods are compared based on the time (in nanoseconds) required to obtain similar standard deviations in the solvation free energies and 1‐octanol/water partition coefficients. By this analysis, the EE method is found to be a factor of nine more efficient than the TI algorithm. For both methods, solvation free energy calculations in 1‐octanol consume roughly an order of magnitude more CPU hours than the hydration free energy calculations. © 2012 Wiley Periodicals, Inc.     

On the transferability of hydration-parametrized continuum electrostatics models to solvated binding calculations

Using molecular mechanics force field partial atomic charges, we show the nonuniqueness of the parametrization of continuum electrostatics models with respect to solute atomic radii and interior dielectric constant based on hydration (vacuum-to-water transfer) free energy data available for small molecules. Moreover, parameter sets that are optimal and equivalent for hydration free energy calculations lead to large variations of calculated absolute and relative electrostatic binding free energies. Hence, parametrization of solvation effects based on hydration data, although a necessary condition, is not sufficient to guarantee its transferability to the calculation of binding free energies in solution.     

Electrostatic solvation free energy of amino acid side chain analogs: implications for the validity of electrostatic linear response in water

Electrostatic free energies of solvation for 15 neutral amino acid side chain analogs are computed. We compare three methods of varying computational complexity and accuracy for three force fields: free energy simulations, Poisson-Boltzmann (PB), and linear response approximation (LRA) using AMBER, CHARMM, and OPLS-AA force fields. We find that deviations from simulation start at low charges for solutes. The approximate PB and LRA produce an overestimation of electrostatic solvation free energies for most of molecules studied here. These deviations are remarkably systematic. The variations among force fields are almost as large as the variations found among methods. Our study confirms that success of the approximate methods for electrostatic solvation free energies comes from their ability to evaluate free energy differences accurately.     

Use of a QM/MM-based FEP method to evaluate the anomalous hydration behavior of simple alkyl amines and amides: application to the design of FBPase inhibitors for the treatment of type-2 diabetes

Standard molecular mechanics (MM) force fields predict a nearly linear decrease in hydration free energy with each successive addition of a methyl group to ammonia or acetamide, whereas a nonadditive relationship is observed experimentally. In contrast, the non-additive hydration behavior is reproduced directly using a quantum mechanics (QM)/MM-based free-energy perturbation (FEP) method wherein the solute partial atomic charges are updated at every window. Decomposing the free energies into electrostatic and van der Waals contributions and comparing the results with the corresponding free energies obtained using a conventional FEP method and a QM/MM method wherein the charges are not updated suggests that inaccuracies in the electrostatic free energies are the primary reason for the inability of the conventional FEP method to predict the experimental findings. The QM/MM-based FEP method was subsequently used to evaluate inhibitors of the diabetes drug target fructose-1,6-bisphosphatase adenosine 5'-monophosphate and 6-methylamino purine riboside 5'-monophosphate. The predicted relative binding free energy was consistent with the experimental findings, whereas the relative binding free energy predicted using the conventional FEP method differed from the experimental finding by an amount consistent with the overestimated relative solvation free energies calculated for alkylamines. Accordingly, the QM/MM-based FEP method offers potential advantages over conventional FEP methods, including greater accuracy and reduced user input. Moreover, since drug candidates often contain either functionality that is inadequately treated by MM (e.g., simple alkylamines and alkylamides) or new molecular scaffolds that require time-consuming development of MM parameters, these advantages could enable future automation of FEP calculations as well as greatly increase the use and impact of FEP calculations in drug discovery.     

Testing and validation of the Automated Topology Builder (ATB) version 2.0: prediction of hydration free enthalpies

To test and validate the Automated force field Topology Builder and Repository (ATB; http://compbio.biosci.uq.edu.au/atb/) the hydration free enthalpies for a set of 214 drug-like molecules, including 47 molecules that form part of the SAMPL4 challenge have been estimated using thermodynamic integration and compared to experiment. The calculations were performed using a fully automated protocol that incorporated a dynamic analysis of the convergence and integration error in the selection of intermediate points. The system has been designed and implemented such that hydration free enthalpies can be obtained without manual intervention following the submission of a molecule to the ATB. The overall average unsigned error (AUE) using ATB 2.0 topologies for the complete set of 214 molecules was 6.7 kJ/mol and for molecules within the SAMPL4 7.5 kJ/mol. The root mean square error (RMSE) was 9.5 and 10.0 kJ/mol respectively. However, for molecules containing functional groups that form part of the main GROMOS force field the AUE was 3.4 kJ/mol and the RMSE was 4.0 kJ/mol. This suggests it will be possible to further refine the parameters provided by the ATB based on hydration free enthalpies.     

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations

We used blind predictions of the 47 hydration free energies in the SAMPL4 challenge to test multiple partial charge models in the context of explicit solvent free energy simulations with the general AMBER force field. One of the partial charge models, IPolQ-Mod, is a fast continuum solvent-based implementation of the IPolQ approach. The AM1-BCC, restrained electrostatic potential (RESP) and IpolQ-Mod approaches all perform reasonably well (R² > 0.8), while VCharge, though faster, gives less accurate results (R² of 0.5). The AM1-BCC results are more accurate than those of RESP for tertiary amines and nitrates, but the overall difference in accuracy between these methods is not statistically significant. Interestingly, the IPolQ-Mod method is found to yield partial charges in very close agreement with RESP. This observation suggests that the success of RESP may be attributed to its fortuitously approximating the arguably more rigorous IPolQ approach.     

Molecular dynamics study of selective adsorption of PCB on activated carbon

The selectivity of PCB adsorption from fish oil onto activated carbon (AC) was investigated by means of molecular dynamics to determine the importance of molecular planarity. PCB congeners 77 and 118 were selected for comparison purposes due to pronounced differences in mean adsorption efficiency and molecular geometry; triolein, a triacylglycerol of oleic acid (C18:1), was used as the representative fish oil component. Graphitic carbon structure was set up to serve as activated carbon model. Molecular force fields employed in the simulations combined short-range parameters from the OPLS with partial atomic charges obtained via quantum chemical calculations using DFT/B3LYP/6-31**G+ and Solvation Model 6. We modified the dihedral angle potential between the PCB aromatic rings and applied Schrödinger's Jaguar package to evaluate the required force field constants. Our complete system comprised a number of PCB molecules dissolved in triacylglycerol that overlaid and filled the pores of an AC structure. The production run of 4 ns provided strong indications that smaller pores will be conductive to better selectivity though also resulted in certain doubts concerning the estimation and assignment of partial atomic charges on the activated carbon. The majority of PCB molecules trapped in pores were attached via cl-AC “bonding”, leaving the main part of the PCB molecule free to interact with triolein. The cl-AC adsorption energy was found to surpass the energy criteria conventionally used for hydrogen bonds. Planar orientation assumed by a PCB molecule in a very energetically favored position on top of the graphite sheet clearly supported the π-cloud overlap hypothesis.     

Improved interaction potentials for charged residues in proteins

Electrostatic interactions dominate the structure and free energy of biomolecules. To obtain accurate free energies involving charged groups from molecular simulations, OPLS-AA parameters have been reoptimized using Monte Carlo free energy perturbation. New parameters fit a self-consistent, experimental set of hydration free energies for acetate (Asp), propionate (Glu), 4-methylimidazolium (Hip), n-butylammonium (Lys), and n-propylguanidinium (Arg), all resembling charged residue side chains, including beta-carbons. It is shown that OPLS-AA free energies depend critically on the type of water model, TIP4P or TIP3P; i.e., each water model requires specific water-charged molecule interaction potentials. New models (models 1 and 3) are thus described for both water models. Uncertainties in relative free energies of charged residues are approximately 2 kcal/mol with the new parameters, due to variations in system setup (MAEs of ca. 1 kcal/mol) and noise from simulations (ca. 1 kcal/mol). The latter error of approximately 1 kcal/mol contrasts MAEs from standard OPLS-AA of up to 13 kcal/mol for the entire series of charged residues or up to 5 kcal/mol for the cationic series Lys, Arg, and Hip. The new parameters can be used directly in molecular simulations with no modification of neutral residues needed and are envisioned to be particular important in simulations where charged residues change environment.     

Prediction of cyclohexane-water distribution coefficients for the SAMPL5 data set using molecular dynamics simulations with the OPLS-AA force field

All-atom molecular dynamics simulations were used to predict water-cyclohexane distribution coefficients \(D_{cw}\) of a range of small molecules as part of the SAMPL5 blind prediction challenge. Molecules were parameterized with the transferable all-atom OPLS-AA force field, which required the derivation of new parameters for sulfamides and heterocycles and validation of cyclohexane parameters as a solvent. The distribution coefficient was calculated from the solvation free energies of the compound in water and cyclohexane. Absolute solvation free energies were computed by an established protocol using windowed alchemical free energy perturbation with thermodynamic integration. This protocol resulted in an overall root mean square error in \(\log D_{cw}\) of almost 4 log units and an overall signed error of ?3 compared to experimental data. There was no substantial overall difference in accuracy between simulating in NVT and NPT ensembles. The signed error suggests a systematic error but the experimental \(D_{cw}\) data on their own are insufficient to uncover the source of this error. Preliminary work suggests that the major source of error lies in the hydration free energy calculations.     

An all-atom force field for metallocenes

A new all-atom force field, for the molecular modeling of metallocenes was constructed. Quantum chemical calculations were performed to obtain several force field terms not yet defined in the literature. The remainder were transferred from the OPLS-AA/AMBER framework. The parametrization work included the obtention of geometrical parameters, torsion energy profiles, and distributions of atomic charges that blend smoothly with the OPLS-AA specification for a variety of organic molecular fragments. Validation was carried out by comparing simulated and experimental data for five different ferrocene derivatives in the crystalline phase. The present model can be regarded as a step toward a general force field for metallocenes, built in a systematic way, easily integrated with OPLS-AA, and transferable between different metal-ligand combinations.     

Free energies of hydration from thermodynamic integration: Comparison of molecular mechanics force fields and evaluation of calculation accuracy

Four commonly used molecular mechanics force fields, CHARMM22, OPLS, CVFF, and GROMOS87, are compared for their ability to reproduce experimental free energies of hydration (ΔG_hydr) from molecular dynamics (MD) simulations for a set of small nonpolar and polar organic molecules: propane, cyclopropane, dimethylether, and acetone. ΔG_hydr values were calculated by multiconfiguration thermodynamic integration for each of the different force fields with three different sets of partial atomic charges: full charges from an electrostatic potential fit (ESP), and ESP charges scaled by 0.8 and 0.6. All force fields, except for GROMOS87, give reasonable results for ΔG_hydr · if partial atomic charges of appropriate magnitude are assigned. For GROMOS87, the agreement with experiment for hydrocarbons (propane and ethane) was improved considerably by modifying the repulsive part of the carbon-water oxygen Lennard-Jones potential. The small molecules studied are related to the chemical moieties constituting camphor (C₁₀H_l6O). By invoking force-field transferability, we calculated the ΔG_hydr for camphor. With the modified GROMOS force field, a ΔG_hydr within 4 kJ/mol of the experimental value of −14.8 kJ/mol was obtained. Camphor is one of the largest molecules for which an absolute hydration free energy has been calculated by molecular simulation. The accuracy and reliability of the thermodynamic integration calculations were analyzed in detail and we found that, for ΔG_hydr calculations for the set of small molecules in aqueous solution, molecular dynamics simulations of 0.8–1.0 ns in length give an upper statistical error bound of 1.5 kJ/mol, whereas shorter simulations of 0.25 nm in length given an upper statistical error bound of 3.5 kJ/mol. © 1997 by John Wiley & Sons, Inc.     

Alchemical prediction of hydration free energies for SAMPL

Hydration free energy calculations have become important tests of force fields. Alchemical free energy calculations based on molecular dynamics simulations provide a rigorous way to calculate these free energies for a particular force field, given sufficient sampling. Here, we report results of alchemical hydration free energy calculations for the set of small molecules comprising the 2011 Statistical Assessment of Modeling of Proteins and Ligands challenge. Our calculations are largely based on the Generalized Amber Force Field with several different charge models, and we achieved RMS errors in the 1.4-2.2 kcal/mol range depending on charge model, marginally higher than what we typically observed in previous studies (Mobley et al. in J Phys Chem B 111(9):2242-2254, 2007, J Chem Theory Comput 5(2):350-358, 2009, J Phys Chem B 115:1329-1332, 2011; Nicholls et al. in J Med Chem 51:769-779, 2008; Klimovich and Mobley in J Comput Aided Mol Design 24(4):307-316, 2010). The test set consists of ethane, biphenyl, and a dibenzyl dioxin, as well as a series of chlorinated derivatives of each. We found that, for this set, using high-quality partial charges from MP2/cc-PVTZ SCRF RESP fits provided marginally improved agreement with experiment over using AM1-BCC partial charges as we have more typically done, in keeping with our recent findings (Mobley et al. in J Phys Chem B 115:1329-1332, 2011). Switching to OPLS Lennard-Jones parameters with AM1-BCC charges also improves agreement with experiment. We also find a number of chemical trends within each molecular series which we can explain, but there are also some surprises, including some that are captured by the calculations and some that are not.     

Accuracy of free energies of hydration for organic molecules from 6-31g*-derived partial charges

Absolute free energies of hydration have been computed for 13 diverse organic molecules using partial charges derived from ab initio 6-31G* wave functions. Both Mulliken charges and charges fit to the electrostatic potential surface (EPS) were considered in conjunction with OPLS Lennard–Jones parameters for the organic molecules and the TIP4P model of water. Monte Carlo simulations with statistical perturbation theory yielded relative free energies of hydration. These were converted to absolute quantities through perturbations to reference molecules for which absolute free energies of hydration had been obtained previously in TIP4P water. The average errors in the computed absolute free energies of hydration are 1.1 kcal/mol for the 6-31G* EPS charges and 4.0 kcal/mol for the Mulliken charges. For the EPS charges, the largest individual errors are under 2 kcal/mol except for acetamide, in which case the error is 3.7 kcal/mol. The hydrogen bonding between the organic solutes and water has also been characterized. © John Wiley & Sons, Inc.     

Efficient calculation of SAMPL4 hydration free energies using OMEGA,SZYBKI, QUACPAC,and Zap TK

Several submissions for the SAMPL4 hydration free energy set were calculated using OpenEye tools, including many that were among the top performing submissions. All of our best submissions used AM1BCC charges and Poisson–Boltzmann solvation. Three submissions used a single conformer for calculating the hydration free energy and all performed very well with mean unsigned errors ranging from 0.94 to 1.08 kcal/mol. These calculations were very fast, only requiring 0.5–2.0 s per molecule. We observed that our two single-conformer methodologies have different types of failure cases and that these differences could be exploited for determining when the methods are likely to have substantial errors.