Solvation Free Energy Differences Research Articles

Calculations of Absolute Solvation Free Energies with Transformato─Application to the FreeSolv Database Using the CGenFF Force Field.

We recently introduced transformato, an open-source Python package for the automated setup of large-scale calculations of relative solvation and binding free energy differences. Here, we extend the capabilities of transformato to the calculation of absolute solvation free energy differences. After careful validation against the literature results and reference calculations with the PERT module of CHARMM, we used transformato to compute absolute solvation free energies for most molecules in the FreeSolv database (621 out of 642). The force field parameters were obtained with the program cgenff (v2.5.1), which derives missing parameters from the CHARMM general force field (CGenFF v4.6). A long-range correction for the Lennard-Jones interactions was added to all computed solvation free energies. The mean absolute error compared to the experimental data is 1.12 kcal/mol. Our results allow a detailed comparison between the AMBER and CHARMM general force fields and provide a more in-depth understanding of the capabilities and limitations of the CGenFF small molecule parameters.

Dissolution of amyloid aggregates by direct addition of alkali halides

The development of standard dissolution agents for various amyloid aggregates in vitro is essential to acquire basic information on pharmaceutical treatment strategies and improve the recovery rate for target expressed proteins. Nevertheless, no standard dissolution agents are currently applied to many amyloid aggregates. Thus, to proceed one step further, the dissolution ability of α-synuclein amyloid aggregates (α-synuclein amyloid) formed under mild conditions (pD 7.5 and 310 K) was investigated using the direct addition of alkali halides. Fourier-transform infrared spectroscopy and Congo red assay were used to evaluate the efficacy of the dissolution agents of alkali halides for amyloid aggregates. The effects of alkali halides on α- synuclein amyloid were classified into three types: (i) dissolution effect, (ii) no effect, and (iii) promotion effect. These effects correlate with the disruption ability of cations (the direct effect) and the solvation free energy difference (ΔΔGSolv) between the cations and anions (the indirect effect). Remarkably, lithium iodide (LiI) has the strongest dissolution ability for α-synuclein amyloid among the examined alkali halides. Furthermore, dialysis can recover the monomer state from α-synuclein amyloid dissolved in LiI. In addition, LiI also exhibits a strong dissolution ability for bovine insulin amyloid aggregates formed under extreme conditions (pD 1.6 and 363 K). These findings show that LiI has the potential as a standard dissolution agent applied to amyloid aggregates among the alkali halides.

Alchemical free energy simulations without speed limits. A generic framework to calculate free energy differences independent of the underlying molecular dynamics program.

We describe the theory of the so‐called common‐core/serial‐atom‐insertion (CC/SAI) approach to compute alchemical free energy differences and its practical implementation in a Python package called Transformato. CC/SAI is not tied to a specific biomolecular simulation program and does not rely on special purpose code for alchemical transformations. To calculate the alchemical free energy difference between several small molecules, the physical end‐states are mutated into a suitable common core. Since this only requires turning off interactions, the setup of intermediate states is straightforward to automate. Transformato currently supports CHARMM and OpenMM as back ends to carry out the necessary molecular dynamics simulations, as well as post‐processing calculations. We validate the method by computing a series of relative solvation free energy differences.

Dummy Atoms in Alchemical Free Energy Calculations.

In calculations of relative free energy differences, the number of atoms of the initial and final states is rarely the same. This necessitates the introduction of dummy atoms. These placeholders interact with the physical system only by bonded energy terms. We investigate the conditions necessary so that the presence of dummy atoms does not influence the result of a relative free energy calculation. On the one hand, one has to ensure that dummy atoms only give a multiplicative contribution to the partition function so that their contribution cancels from double-free energy differences. On the other hand, the bonded terms used to attach a dummy atom (or group of dummy atoms) to the physical system have to maintain it in a well-defined position and orientation relative to the physical system. A detailed theoretical analysis of both aspects is provided, illustrated by 24 calculations of relative solvation free energy differences, for which all four legs of the underlying thermodynamic cycle were computed. Cycle closure (or lack thereof) was used as a sensitive indicator to probing the effects of dummy atom treatment on the resulting free energy differences. We find that a naive (but often practiced) treatment of dummy atoms results in errors of up to kBT when calculating the relative solvation free energy difference between two small solutes, such as methane and ammonia. While our analysis focuses on the so-called single topology approach to set up alchemical transformations, similar considerations apply to dual topology, at least many widely used variants thereof.

Implicit Solvents for the Polarizable Atomic Multipole AMOEBA Force Field.

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the nonlinear and linearized versions of the Poisson-Boltzmann equation (PBE), the domain-decomposition conductor-like screening model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatics models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least-squares style target function based on a library of 103 small-molecule solvation free energy differences. Mean signed errors for the adaptive Poisson-Boltzmann solver (APBS), ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.58 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields opens the door to their use for folding and design applications.

Bisphenol-S and Bisphenol-F alter mouse pancreatic β-cell ion channel expression and activity and insulin release through an estrogen receptor ERβ mediated pathway

Bisphenol-S (BPS) and Bisphenol-F (BPF) are current Bisphenol-A (BPA) substitutes. Here we used pancreatic β-cells from wild type (WT) and estrogen receptor β (ERβ) knockout (BERKO) mice to investigate the effects of BPS and BPF on insulin secretion, and the expression and activity of ion channels involved in β-cell function. BPS or BPF rapidly increased insulin release and diminished ATP-sensitive K+ (KATP) channel activity. Similarly, 48 h treatment with BPS or BPF enhanced insulin release and decreased the expression of several ion channel subunits in β-cells from WT mice, yet no effects were observed in cells from BERKO mice. PaPE-1, a ligand designed to preferentially trigger extranuclear-initiated ER pathways, mimicked the effects of bisphenols, suggesting the involvement of extranuclear-initiated ERβ pathways. Molecular dynamics simulations indicated differences in ERβ ligand-binding domain dimer stabilization and solvation free energy among different bisphenols and PaPE-1. Our data suggest a mode of action involving ERβ whose activation alters three key cellular events in β-cell, namely ion channel expression and activity, and insulin release. These results may help to improve the hazard identification of bisphenols.

Solubility measurement, model evaluation and molecular simulations of (R)-(-)-phenylephrine hydrochloride in three binary solvents

In this research work, laser monitoring method was employed to investigate the mole-fraction solubility data of α adrenergic receptor agonists (R)-(-)-phenylephrine hydrochloride (R-PEHC) in three binary solvents namely “methanol (MtOH) + ethanol (EtOH), MtOH + isopropanol (IPOH) and EtOH + IPOH” over the temperature of (278.15–323.15) K with temperature range of 5 K under 0.1 MPa. The results demonstrated that experimental solubility was positively related with temperature while inversely associated with mass fraction of anti-solvents. At the same temperature and mass fraction of anti-solvent, the general solubility order of R-PEHC in binary solvents was MtOH + EtOH > MtOH + IPOH > EtOH + IPOH. Meanwhile, Hansen solubility parameter (HSP) was employed to investigate the solubility order, and the results showed solubility behavior could be analyzed well by selected four solubility parameters (Δδp, Δδh, Δδt and Δδ). Besides, the molecular simulation technology was conducted to investigate solvation free energy (Esol) between R-PEHC and selected solvents. The results showed that experimental solubility differences in various solvents were attributed to the differences in solvation free energy. Moreover, mathcad software (Version: 15) was applied to correlate experimental solubility with respect of temperatures by using five thermodynamic models (UNIQUAC, Wilson, Three-Suffix Margules, Jouyban-Acree-modified Apelblat (J-A-A) and NRTL model). All selected models provided satisfactory goodness-of-fit with R-PEHC solubility data. Furthermore, apparent thermodynamic properties of solutions were calculated based on van’t Hoff equation as well as experimental solubility, the results indicated that dissolution process was always endothermic and entropy-driven.

Predicting octanol/water partition coefficients for the SAMPL6 challenge using the SM12, SM8, and SMD solvation models.

Blind predictions of octanol/water partition coefficients at 298 K for 11 kinase inhibitor fragment like compounds were made for the SAMPL6 challenge. We used the conventional, "untrained", free energy based approach wherein the octanol/water partition coefficient was computed directly as the difference in solvation free energy in water and 1-octanol. We additionally proposed and used two different forms of a "trained" approach. Physically, the goal of the trained approach is to relate the partition coefficient computed using pure 1-octanol to that using water-saturated 1-octanol. In the first case, we assumed the partition coefficient using water-saturated 1-octanol and pure 1-octanol are linearly correlated. In the second approach, we assume the solvation free energy in water-saturated 1-octanol can be written as a linear combination of the solvation free energy in pure water and 1-octanol. In all cases here, the solvation free energies were computed using electronic structure calculations in the SM12, SM8, and SMD universal solvent models. In the context of the present study, our results in general do not support the additional effort of the trained approach.

LogP prediction performance with the SMD solvation model and the M06 density functional family for SAMPL6 blind prediction challenge molecules.

This work presents a quantum mechanical model for predicting octanol-water partition coefficients of small protein-kinase inhibitor fragments as part of the SAMPL6 LogP Prediction Challenge. The model calculates solvation free energy differences using the M06-2X functional with SMD implicit solvation and the def2-SVP basis set. This model was identified as dqxk4 in the SAMPL6 Challenge and was the third highest performing model in the physical methods category with 0.49 log Root Mean Squared Error (RMSE) for predicting the 11 compounds in SAMPL6 blind prediction set. We also collaboratively investigated the use of empirical models to address model deficiencies for halogenated compounds at minimal additional computational cost. A mixed model consisting of the dqxk4 physical and hdpuj empirical models found improved performance at 0.34 log RMSE on the SAMPL6 dataset. This collaborative mixed model approach shows how empirical models can be leveraged to expediently improve performance in chemical spaces that are difficult for ab initio methods to simulate.

Solvent effects on a derivative of 1,3,4-oxadiazole tautomerization reaction in water: A reaction density functional theory study

In this work, we investigate the simple and water assisted tautomerization reaction in a derivative of 1,3,4-oxadiazole by means of the proposed multiscale reaction density functional theory (RxDFT). The complete scheme and reaction pathways among four different tautomers are investigated, and the free energy profiles of four corresponding tautomerization reactions in the absence and presence of 1–3 water molecules are calculated. The results show that the presence of liquid water considerably reduces the free energy barriers of the tautomerization reactions, but imposes negligible impact on the reaction free energies. The solvent effects are ascribed primarily to the assistance of water molecules, and then to the difference in solvation free energy. The satisfactory accuracy and low computational cost of RxDFT highlight its great potential for studying solvent effect on chemical reactions in solution compared with ab initio or QM/MM method.

Quadrupole Correction: From Molecular Electrostatic Potential to Free Energies of Halogen Bonding.

Recently several molecular mechanics models of halogen bonding have been published. They describe the electrostatic potential anisotropy near the heavy halogen atom (known as a σ-hole) in different ways, ranging from an all-atom multipole expansion to a single positive extra-point charge. However, the question of a reasonable balance between the accuracy and the simplicity of the model remains open. In this work, we introduce the simplistic RESPQ electrostatics model built on the RESP charges complemented with fixed atomic quadrupoles. We show that it: (1) correctly describes the MEP anisotropy of aromatic halogen atoms, (2) improves the description of the halogen-water interaction energies in both halogen and hydrogen bonding cases, (3) provides an excellent estimation of solvation free energy differences of aromatic halogens, and (4) is compatible (with the help of multipole charge cluster approximation) with contemporary molecular modeling packages.

Genistein Binding to Copper(II)-Solvent Dependence and Effects on Radical Scavenging.

Genistein, but not daidzein, binds to copper(II) with a 1:2 stoichiometry in ethanol and with a 1:1 stoichiometry in methanol, indicating chelation by the 5-phenol and the 4-keto group of the isoflavonoid as demonstrated by the Jobs method and UV-visible absorption spectroscopy. In ethanol, the stability constants had the value 1.12 × 1011 L2∙mol−2 for the 1:2 complex and in methanol 6.0 × 105 L∙mol−1 for the 1:1 complex at 25 °C. Binding was not detected in water, as confirmed by an upper limit for the 1:1 stability constant of K = 5 mol−1 L as calculated from the difference in solvation free energy of copper(II) between methanol and the more polar water. Solvent molecules compete with genistein as demonstrated in methanol where binding stoichiometry changes from 1:2 to 1:1 compared to ethanol and methanol/chloroform (7/3, v/v). Genistein binding to copper(II) increases the scavenging rate of the stable, neutral 2,2-diphenyl-1-picrylhydrazyl radical by more than a factor of four, while only small effects were seen for the short-lived but more oxidizing β-carotene radical cation using laser flash photolysis. The increased efficiency of coordinated genistein is concluded to depend on kinetic rather than on thermodynamic factors, as confirmed by the small change in reduction potential of −0.016 V detected by cyclic voltammetry upon binding of genistein to copper(II) in methanol/chloroform solutions.

Photoisomers of Azobenzene Star with a Flat Core: Theoretical Insights into Multiple States from DFT and MD Perspective.

This study focuses on comparing physical properties of photoisomers of an azobenzene star with benzene-1,3,5-tricarboxamide core. Three azobenzene arms of the molecule undergo a reversible trans-cis isomerization upon UV-vis light illumination giving rise to multiple states from the planar all-trans one, via two mixed states to the kinked all-cis isomer. Employing density functional theory, we characterize the structural and photophysical properties of each state indicating a role the planar core plays in the coupling between azobenzene chromophores. To characterize the light-triggered switching of solvophilicity/solvophobicity of the star, the difference in solvation free energy is calculated for the transfer of an azobenzene star from its gas phase to implicit or explicit solvents. For the latter case, classical all-atom molecular dynamics simulations of aqueous solutions of azobenzene star are performed employing the polymer consistent force field to shed light on the thermodynamics of explicit hydration as a function of the isomerization state and on the structuring of water around the star. From the analysis of two contributions to the free energy of hydration, the nonpolar van der Waals and the electrostatic terms, it is concluded that isomerization specificity largely determines the polarity of the molecule and the solute-solvent electrostatic interactions. This convertible hydrophilicity/hydrophobicity together with readjustable occupied volume and the surface area accessible to water, affects the self-assembly/disassembly of the azobenzene star with a flat core triggered by light.

Thermodynamic and Structural Factors That Influence the Redox Potentials of Tungsten–Alkylidyne Complexes

The thermodynamic and structural factors that influence the redox properties of an extensive set of tungsten–alkylidyne complexes (W(CR)L4X) are analyzed by combining synthesis, electrochemistry, and computational modeling based on free energy calculations of oxidation potentials at the density functional theory level. The observed linear correlations among oxidation potentials, HOMO energies, and gas-phase ionization energies are found to be consistent with the approximately constant solvation free energy differences between reduced and oxidized species over the complete set. The W–X bond length, trans to the alkylidyne ligand, is found to be a good descriptor of the positioning of the key frontier orbitals that regulate the redox properties of the complexes.

Cluster expansion of the solvation free energy difference: Systematic improvements in the solvation of single ions.

The cluster expansion method has been used in the imperfect gas theory for several decades. This paper proposes a cluster expansion of the solvation free energy difference. This difference, which results from a change in the solute-solvent potential energy, can be written as the logarithm of a finite series. Similar to the Mayer function, the terms in the series are related to configurational integrals, which makes the integrand relevant only for configurations of the solvent molecules close to the solute. In addition, the terms involve interaction of solute with one, two, and so on solvent molecules. The approach could be used for hybrid quantum mechanical and molecular mechanics methods or mixed cluster-continuum approximation. A simple form of the theory was applied for prediction of pKa in methanol; the results indicated that three explicit methanol molecules and the dielectric continuum lead to a root of mean squared error (RMSE) of only 1.3 pKa units, whereas the pure continuum solvation model based on density method leads to a RMSE of 6.6 pKa units.

The free energy of locking a ring: Changing a deoxyribonucleoside to a locked nucleic acid

Locked nucleic acid (LNA), a modified nucleoside which contains a bridging group across the ribose ring, improves the stability of DNA/RNA duplexes significantly, and therefore is of interest in biotechnology and gene therapy applications. In this study, we investigate the free energy change between LNA and DNA nucleosides. The transformation requires the breaking of the bridging group across the ribose ring, a problematic transformation in free energy calculations. To address this, we have developed a 3‐step (easy to implement) and a 1‐step protocol (more efficient, but more complicated to setup), for single and dual topologies in classical molecular dynamics simulations, using the Bennett Acceptance Ratio method to calculate the free energy. We validate the approach on the solvation free energy difference for the nucleosides thymidine, cytosine, and 5‐methyl‐cytosine. © 2017 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Resolving solvophobic interactions inferred from experimental solvation free energies and evaluated from molecular simulations

Ben-Naim estimated the solvent-mediated interaction between methanes based on experimental solvation free energy differences between chemically similar hydrocarbons. Interactions were predicted to be strongest in water, dominated by characteristic entropic gains. We use molecular simulations in combination with an empirical interpolation procedure that bridges interactions from outside methane’s excluded volume down to overlap to test Ben-Naim’s estimates. While Ben-Naim’s approach captures many distinctive trends, the alchemical differences between methane and a methyl unit play a non-trivial role on the predicted association strength and the sign of enthalpic and entropic components of the interaction free energy in water and ethanol.

Statins in therapy: Understanding their hydrophilicity, lipophilicity, binding to 3-hydroxy-3-methylglutaryl-CoA reductase, ability to cross the blood brain barrier and metabolic stability based on electrostatic molecular orbital studies

The atomic electrostatic potentials calculated by the CHELPG method have been shown to be sensitive indicators of the gas phase and solution properties of the statins. Solvation free energies in water, n-octanol and n-octane have been determined using the SMD solvent model. The percentage hydrophilicity and hydrophobicity (or lipophilicity) of the statins in solution have been determined using (a) the differences in solvation free energies between n-octanol and n-octane as a measure of hydrophilicity, and the solvation energy in octane as a measure of hydrophobicity (b) the sum of the atomic electrostatic charges on the hydrogen bonding and polar bonding nuclei of the common pharmacophore combined with a solvent measure of hydrophobicity, and (c) using the buried surface areas after statin binding to HMGCR to calculate the hydrophobicity of the bound statins. The data suggests that clinical definitions of statins as either “hydrophilic” or “lipophilic” based on experimental partition coefficients are misleading.An estimate of the binding energy between rosuvastatin and HMGCR has been made using: (a) a coulombic electrostatic interaction model, (b) the calculated desolvation and resolvation of the statin in water, and (c) the first shell transfer solvation energy as a proxy for the restructuring of the water molecules immediately adjacent to the active binding site of HMGCR prior to binding. Desolvation and resolvation of the statins before and after binding to HMGCR are major determinants of the energetics of the binding process.An analysis of the amphiphilic nature of lovastatin anion, acid and lactone and fluvastatin anion and their abilities to cross the blood brain barrier has indicated that this process may be dominated by desolvation and resolvation effects, rather than the statin molecular size or statin–lipid interactions within the bilayer.The ionization energy and electron affinity of the statins are sensitive physical indicators of the ease that the various statins can undergo endogenous oxidative metabolism. The absolute chemical hardness is also an indicator of the stability of the statins, and may be a useful indicator for drug design.

Multiscale Free Energy Simulations: An Efficient Method for Connecting Classical MD Simulations to QM or QM/MM Free Energies Using Non-Boltzmann Bennett Reweighting Schemes.

The reliability of free energy simulations (FES) is limited by two factors: (a) the need for correct sampling and (b) the accuracy of the computational method employed. Classical methods (e.g., force fields) are typically used for FES and present a myriad of challenges, with parametrization being a principle one. On the other hand, parameter-free quantum mechanical (QM) methods tend to be too computationally expensive for adequate sampling. One widely used approach is a combination of methods, where the free energy difference between the two end states is computed by, e.g., molecular mechanics (MM), and the end states are corrected by more accurate methods, such as QM or hybrid QM/MM techniques. Here we report two new approaches that significantly improve the aforementioned scheme; with a focus on how to compute corrections between, e.g., the MM and the more accurate QM calculations. First, a molecular dynamics trajectory that properly samples relevant conformational degrees of freedom is generated. Next, potential energies of each trajectory frame are generated with a QM or QM/MM Hamiltonian. Free energy differences are then calculated based on the QM or QM/MM energies using either a non-Boltzmann Bennett approach (QM-NBB) or non-Boltzmann free energy perturbation (NB-FEP). Both approaches are applied to calculate relative and absolute solvation free energies in explicit and implicit solvent environments. Solvation free energy differences (relative and absolute) between ethane and methanol in explicit solvent are used as the initial test case for QM-NBB. Next, implicit solvent methods are employed in conjunction with both QM-NBB and NB-FEP to compute absolute solvation free energies for 21 compounds. These compounds range from small molecules such as ethane and methanol to fairly large, flexible solutes, such as triacetyl glycerol. Several technical aspects were investigated. Ultimately some best practices are suggested for improving methods that seek to connect MM to QM (or QM/MM) levels of theory in FES.

Is the Isodesmic Reaction Approach a Better Model for Accurate Calculation of pK a of Organic Superbases? A Computational Study

The acid-base dissociation constant (p<i>K</i> _a) can be related to the solubility and binding of drugs. However, measuring accurate p<i>K</i> _a values is a challenging task. In this study, we have examined the p<i>K</i> _a of various organic superbases: naphthalenes, cyclic guanidines, vinamidines, and acyclic guanidines computationally. We have calculated the p<i>K</i> _a of such superbases by employing two methods: a conventional thermodynamic cycle and a second method based on an isodesmic reaction. The thermodynamic cycle involves computation of solvation free energy by using gas-phase free energy and the difference in solvation free energies (∆G_solv) between products and reactants. Calculations performed with the isodesmic reaction approach do not use the free energy of solvation; hence, the accuracy of the approach is less sensitive to solvent molecules and global charges of the calculated species. The root-mean-square errors (RMSE) predict that the p<i>K</i> _a of the studied organic superbases are more accurate when calculated with the isodesmic reaction approach.