Hydration Free Energy Calculations Research Articles

Calculated hydration free energies become less accurate with increases in molecular weight

In order for computer-aided drug design to fulfil its long held promise of delivering new medicines faster and cheaper, extensive development and validation work must be done first. This pertains particularly to molecular dynamics force fields where one important aspect–the hydration free energy (HFE) of small molecules–is often insufficiently analyzed. While most benchmarking studies report excellent accuracies of calculated hydration free energies–usually within 2 kcal/mol of experimental values–we find that deeper analysis reveals significant shortcomings. Herein, we report a dependence of HFE prediction errors on ligand molecular weight–the higher the weight, the bigger the prediction error and the higher the probability the calculated result is erroneous by a large amount. We show that in the drug-like molecular weight region, HFE predictions can easily be off by 5 kcal/mol or more. This is likely to be highly problematic in a drug discovery and development setting. We make our HFE results and molecular descriptors freely and fully available in order to encourage deeper analysis of future molecular dynamics results and facilitate development of the next generation of force fields.

Kartograf: A Geometrically Accurate Atom Mapper for Hybrid-Topology Relative Free Energy Calculations.

Relative binding free energy (RBFE) calculations have emerged as a powerful tool that supports ligand optimization in drug discovery. Despite many successes, the use of RBFEs can often be limited by automation problems, in particular, the setup of such calculations. Atom mapping algorithms are an essential component in setting up automatic large-scale hybrid-topology RBFE calculation campaigns. Traditional algorithms typically employ a 2D subgraph isomorphism solver (SIS) in order to estimate the maximum common substructure. SIS-based approaches can be limited by time-intensive operations and issues with capturing geometry-linked chemical properties, potentially leading to suboptimal solutions. To overcome these limitations, we have developed Kartograf, a geometric-graph-based algorithm that uses primarily the 3D coordinates of atoms to find a mapping between two ligands. In free energy approaches, the ligand conformations are usually derived from docking or other previous modeling approaches, giving the coordinates a certain importance. By considering the spatial relationships between atoms related to the molecule coordinates, our algorithm bypasses the computationally complex subgraph matching of SIS-based approaches and reduces the problem to a much simpler bipartite graph matching problem. Moreover, Kartograf effectively circumvents typical mapping issues induced by molecule symmetry and stereoisomerism, making it a more robust approach for atom mapping from a geometric perspective. To validate our method, we calculated mappings with our novel approach using a diverse set of small molecules and used the mappings in relative hydration and binding free energy calculations. The comparison with two SIS-based algorithms showed that Kartograf offers a fast alternative approach. The code for Kartograf is freely available on GitHub (https://github.com/OpenFreeEnergy/kartograf). While developed for the OpenFE ecosystem, Kartograf can also be utilized as a standalone Python package.

Identifying Systematic Force Field Errors Using a 3D-RISM Element Counting Correction.

Hydration free energies of small molecules are commonly used as benchmarks for solvation models. However, errors in predicting hydration free energies are partially due to the force fields used and not just the solvation model. To address this, we have used the 3D reference interaction site model (3D-RISM) of molecular solvation and existing benchmark explicit solvent calculations with a simple element count correction (ECC) to identify problems with the non-bond parameters in the general AMBER force field (GAFF). 3D-RISM was used to calculate hydration free energies of all 642 molecules in the FreeSolv database, and a partial molar volume correction (PMVC), ECC, and their combination (PMVECC) were applied to the results. The PMVECC produced a mean unsigned error of 1.01±0.04kcal/mol and root mean squared error of 1.44±0.07kcal/mol, better than the benchmark explicit solvent calculations from FreeSolv, and required less than 15 s of computing time per molecule on a single CPU core. Importantly, parameters for PMVECC showed systematic errors for molecules containing Cl, Br, I, and P. Applying ECC to the explicit solvent hydration free energies found the same systematic errors. The results strongly suggest that some small adjustments to the Lennard-Jones parameters for GAFF will lead to improved hydration free energy calculations for all solvent models.

The Ion-Dipole Correction of the 3DRISM Solvation Model to Accurately Compute Water Distributions around Negatively Charged Biomolecules.

The 3D reference interaction site model (3DRISM) provides an efficient grid-based solvation model to compute the structural and thermodynamic properties of biomolecules in aqueous solutions. However, it remains challenging for existing 3DRISM methods to correctly predict water distributions around negatively charged solute molecules. In this paper, we first show that this challenge is mainly due to the orientation of water molecules in the first solvation shell of the negatively charged solute molecules. To properly consider this orientational preference, position-dependent two-body intramolecular correlations of solvent need to be included in the 3DRISM theory, but direct evaluations of these position-dependent two-body intramolecular correlations remain numerically intractable. To address this challenge, we introduce the Ion-Dipole Correction (IDC) to the 3DRISM theory, in which we incorporate the orientation preference of water molecules via an additional solute-solvent interaction term (i.e., the ion-dipole interaction) while keeping the formulism of the 3DRISM equation unchanged. We prove that this newly introduced IDC term is equivalent to an effective direct correlation function which can effectively consider the orientation effect that arises from position dependent two-body correlations. We first quantitatively validate our 3DRISM-IDC theory combined with the PSE3 closure on Cl-, [ClO]- (a two-site anion), and [NO2]- (a three-site anion). For all three anions, we show that our 3DRISM-IDC theory significantly outperforms the 3DRISM theory in accurately predicting the solvation structures in comparison to MD simulations, including RDFs and 3D water distributions. Furthermore, we have also demonstrated that the 3DRISM-IDC can improve the accuracy of hydration free-energy calculation for Cl-. We further demonstrate that our 3DRISM-IDC theory yields significant improvements over the 3DRISM theory when applied to compute the solvation structures for various negatively charged solute molecules, including adenosine triphosphate (ATP), a short peptide containing 19 residues, a DNA hairpin containing 24 nucleotides, and a riboswitch RNA molecule with 77 nucleotides. We expect that our 3DRISM-IDC-PSE3 solvation model holds great promise to be widely applied to study solvation properties for nucleic acids and other biomolecules containing negatively charged functional groups.

RestraintMaker: a graph-based approach to select distance restraints in free-energy calculations with dual topology

The calculation of relative binding free energies (RBFE) involves the choice of the end-state/system representation, of a sampling approach, and of a free-energy estimator. System representations are usually termed “single topology” or “dual topology”. As the terminology is often used ambiguously in the literature, a systematic categorization of the system representations is proposed here. In the dual-topology approach, the molecules are simulated as separate molecules. Such an approach is relatively easy to automate for high-throughput RBFE calculations compared to the single-topology approach. Distance restraints are commonly applied to prevent the molecules from drifting apart, thereby improving the sampling efficiency. In this study, we introduce the program RestraintMaker, which relies on a greedy algorithm to find (locally) optimal distance restraints between pairs of atoms based on geometric measures. The algorithm is further extended for multi-state methods such as enveloping distribution sampling (EDS) or multi-site lambda -dynamics. The performance of RestraintMaker is demonstrated for toy models and for the calculation of relative hydration free energies. The Python program can be used in script form or through an interactive GUI within PyMol. The selected distance restraints can be written out in GROMOS or GROMACS file formats. Additionally, the program provides a human-readable JSON format that can easily be parsed and processed further. The code of RestraintMaker is freely available on GitHub https://github.com/rinikerlab/restraintmaker.

Efficient Alchemical Intermediate States in Free Energy Calculations Using λ-Enveloping Distribution Sampling.

Alchemical free energy calculations generally require intermediate states along a coupling parameter λ to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced λ-enveloping distribution sampling (λ-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional λ-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of λ-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules.

Scaffold Hopping Transformations Using Auxiliary Restraints for Calculating Accurate Relative Binding Free Energies.

In silico screening of drug-target interactions is a key part of the drug discovery process. Changes in the drug scaffold via contraction or expansion of rings, the breaking of rings, and the introduction of cyclic structures from acyclic structures are commonly applied by medicinal chemists to improve binding affinity and enhance favorable properties of candidate compounds. These processes, commonly referred to as scaffold hopping, are challenging to model computationally. Although relative binding free energy (RBFE) calculations have shown success in predicting binding affinity changes caused by perturbing R-groups attached to a common scaffold, applications of RBFE calculations to modeling scaffold hopping are relatively limited. Scaffold hopping inevitably involves breaking and forming bond interactions of quadratic functional forms, which is highly challenging. A novel method for handling ring opening/closure/contraction/expansion and linker contraction/expansion is presented here. To the best of our knowledge, RBFE calculations on linker contraction/expansion have not been previously reported. The method uses auxiliary restraints to hold the atoms at the ends of a bond in place during the breaking and forming of the bonds. The broad applicability of the method was demonstrated by examining perturbations involving small-molecule macrocycles and mutations of proline in proteins. High accuracy was obtained using the method for most of the perturbations studied. The rigor of the method was isolated from the force field by validating the method using relative and absolute hydration free energy calculations compared to standard simulation results. Unlike other methods that rely on λ-dependent functional forms for bond interactions, the method presented here can be employed using modern molecular dynamics software without modification of codes or force field functions.

First-principles hydration free energies of oxygenated species at water-platinum interfaces.

The hydration free energy of atoms and molecules adsorbed at liquid-solid interfaces strongly influences the stability and reactivity of solid surfaces. However, its evaluation is challenging in both experiments and theories. In this work, a machine learning aided molecular dynamics method is proposed and applied to oxygen atoms and hydroxyl groups adsorbed on Pt(111) and Pt(100) surfaces in water. The proposed method adopts thermodynamic integration with respect to a coupling parameter specifying a path from well-defined non-interacting species to the fully interacting ones. The atomistic interactions are described by a machine-learned inter-atomic potential trained on first-principles data. The free energy calculated by the machine-learned potential is further corrected by using thermodynamic perturbation theory to provide the first-principles free energy. The calculated hydration free energies indicate that only the hydroxyl group adsorbed on the Pt(111) surface attains a hydration stabilization. The observed trend is attributed to differences in the adsorption site and surface morphology.

Alchemical transformations for concerted hydration free energy estimation with explicit solvation.

We present a family of alchemical perturbation potentials that enable the calculation of hydration free energies of small- to medium-sized molecules in a single concerted alchemical coupling step instead of the commonly used sequence of two distinct coupling steps for Lennard-Jones and electrostatic interactions. The perturbation potentials we employ are non-linear functions of the solute-solvent interaction energy designed to focus sampling near entropic bottlenecks along the alchemical pathway. We present a general framework to optimize the parameters of alchemical perturbation potentials of this kind. The optimization procedure is based on the λ-function formalism and the maximum-likelihood parameter estimation procedure we developed earlier to avoid the occurrence of multi-modal distributions of the coupling energy along the alchemical path. A novel soft-core function applied to the overall solute-solvent interaction energy rather than individual interatomic pair potentials critical for this result is also presented. Because it does not require modifications of core force and energy routines, the soft-core formulation can be easily deployed in molecular dynamics simulation codes. We illustrate the method by applying it to the estimation of the hydration free energy in water droplets of compounds of varying size and complexity. In each case, we show that convergence of the hydration free energy is achieved rapidly. This work paves the way for the ongoing development of more streamlined algorithms to estimate free energies of molecular binding with explicit solvation.

Ligand Affinities within the Open-Boundary Molecular Mechanics/Coarse-Grained Framework (I): Alchemical Transformations within the Hamiltonian Adaptive Resolution Scheme.

Our recently developed Open-Boundary Molecular Mechanics/Coarse Grained (OB-MM/CG) framework predicts ligand poses in important pharmaceutical targets, such as G-protein Coupled Receptors, even when experimental structural information is lacking. The approach, which is based on GROMOS and AMBER force fields, allows for grand-canonical simulations of protein-ligand complexes by using the Hamiltonian Adaptive Resolution Scheme (H-AdResS) for the solvent. Here, we present a key step toward the estimation of ligand binding affinities for their targets within this approach. This is the implementation of the H-AdResS in the GROMACS code. The accuracy of our implementation is established by calculating hydration free energies of several molecules in water by means of alchemical transformations. The deviations of the GROMOS- and AMBER-based H-AdResS results from the reference fully atomistic simulations are smaller than the accuracy of the force field and/or they are in the range of the published results. Importantly, our predictions are in good agreement with experimental data. The current implementation paves the way to the use of the OB-MM/CG framework for the study of large biological systems.

Optimization of Alchemical Pathways Using Extended Thermodynamic Integration.

Thermodynamic integration (TI) is a commonly used method to determine free-energy differences. One of its disadvantages is that many intermediate λ-states need to be sampled in order to be able to integrate accurately over ⟨∂H/∂λ⟩. Here, we use the recently introduced extended TI to study alternative parameterizations of H(λ) and its influence on the smoothness of the ⟨∂H/∂λ⟩ curves as well as the efficiency of the simulations. We find that the extended TI approach can be used to select curves of low curvature. An optimal parameterization is suggested for the calculation of hydration free energies. For calculations of relative binding free energies, we show that optimized parameterizations of the Hamiltonian in the unbound state also effectively lower the curvature in the bound state of the ligand.

Dialing in Direct Air Capture of CO2 by Crystal Engineering of Bisiminoguanidines.

Direct air capture (DAC) technologies that extract carbon dioxide from the atmosphere via chemical processes have the potential to restore the atmospheric CO2 concentration to an optimal level. This study elucidates structure-property relationships in DAC by crystallization of bis(iminoguanidine) (BIG) carbonate salts. Their crystal structures are analyzed by X-ray and neutron diffraction to accurately measure key structural parameters including molecular conformations, hydrogen bonding, and π-stacking. Experimental measurements of key properties, such as aqueous solubilities and regeneration energies and temperatures, are complemented by first-principles calculations of lattice and hydration free energies, as well as free energies of reactions with CO2, and BIG regenerations. Minor structural modifications in the molecular structure of the BIGs are found to result in major changes in the crystal structures and the aqueous solubilities within the series, leading to enhanced DAC.

Accurate and rapid calculation of hydration free energy and its physical implication for biomolecular functions.

Here we review a new method for calculating a hydration free energy (HFE) of a solute and discuss its physical implication for biomolecular functions in aqueous environments. The solute hydration is decomposed into processes 1 and 2. A cavity matching the geometric characteristics of the solute at the atomic level is created in process 1. Solute-water van der Waals and electrostatic interaction potentials are incorporated in process 2. The angle-dependent integral equation theory combined with our morphometric approach is applied to process 1, and the three-dimensional reference interaction site model theory is employed for process 2. Molecular models are adopted for water. The new method is characterized by the following. Solutes with various sizes including proteins can be treated in the same manner. It is almost as accurate as the molecular dynamics simulation despite its far smaller computational burden. It enables us to handle a solute possessing a significantly large total charge without difficulty. The HFE can be decomposed into a variety of physically insightful, energetic, and entropic components. It is best suited to the elucidation of mechanisms of protein folding, pressure and cold denaturation of a protein, and different types of molecular recognition.

Alchemical Hydration Free-Energy Calculations Using Molecular Dynamics with Explicit Polarization and Induced Polarity Decoupling: An On-the-Fly Polarization Approach.

We present a methodology using fixed charge force fields for alchemical solvation free energy calculations which accounts for the change in polarity that the solute experiences as it transfers from the gas-phase to the condensed phase. We update partial charges using QM/MM snapshots, decoupling the electric field appropriately when updating the partial charges. We also show how to account for the cost of self-polarization. We test our methodology on 30 molecules ranging from small polar to large druglike molecules. We use Minimum Basis Iterative Stockholder (MBIS), Restrained Electrostatic Potential (RESP), and AM1-BCC partial charge methodologies. Using our method with MP2/cc-pVTZ and MBIS partial charges yields an average absolute deviation (AAD) of 6.3 kJ·mol-1 in comparison with the AM1-BCC result of 8.6 kJ·mol-1. AM1-BCC is within experimental uncertainty on 10% of the data compared to 30% with our method. We conjecture that results can be further improved by using Lennard-Jones and torsional parameters refitted to MBIS and RESP partial charge methods that use high levels of theory.

Molecular dynamics simulations suggest conformational and hydration difference between zwitterionic poly (carboxybetaine methacrylate) and poly (ethylene glycol)

Molecular dynamics simulations were performed to characterize the conformational features and hydration dynamics of zwitterionic poly(carboxybetaine methacrylate) (pCBMA). Poly(ethylene glycol)(PEG) with comparable molecular weight was taken as referential counterpart. Conformational analyses showed that PEG chain was totally flexible and adopted a curled conformation. Massive segment-to-segment interactions in PEG could result in a hindered hydration effect. While pCBMA had hierarchic flexibility, the fixed backbone and relatively free sidechains, which provided an exclusive hydration volume of each zwitterionic sidechain. In hydration dynamics, pCBMA could largely delay and strongly hydrogen-bond with the surrounding waters. While PEG could suppress the reorientation dynamics of surficial waters due to the enclosed sites and the overlapped hydration volume. Hydration free energy calculations using umbrella sampling revealed the better water affinity of pCBMA against PEG. PEG had high cohesive energy density and needed energy to fully relax upon hydration.

Implicit Model for the Hydration Free Energy Calculation in the Task of the Supercomputer Docking

An implicit solvent model for task of the supercomputer docking is developed. Model is parameterized using the set of the 321 organic molecules of all chemical classes. The accuracy of the present model is higher than the accuracy of the previous model, used in the supercomputer docking.

Thermodynamic Stability of Human γD-Crystallin Mutants Using Alchemical Free-Energy Calculations.

γD-Crystallin (HγDC) is a key structural protein in the human lens, whose aggregation has been associated with the development of cataracts. Single-point mutations and post-translational modifications destabilize HγDC interactions, forming partially folded intermediates, where hydrophobic residues are exposed and thus triggering its aggregation. In this work, we used alchemical free-energy calculations to predict changes in thermodynamic stability (ΔΔG) of 10 alanine-scanning variants and 12 HγDC mutations associated with the development of congenital cataract. Our results show that W42R is the most destabilizing mutation in HγDC. This has been corroborated through experimental determination of ΔΔG employing differential scanning calorimetry. Calculations of hydration free energies from the HγDC WT and the W42R mutant suggested that the mutant has a higher aggregation propensity. Our combined theoretical and experimental results contribute to understand HγDC destabilization and aggregation mechanisms in age-onset cataracts.

An accurate and rapid method for calculating hydration free energies of a variety of solutes including proteins.

A new method is developed for calculating hydration free energies (HFEs) of polyatomic solutes. The solute insertion is decomposed into the creation of a cavity in water matching the geometric characteristics of the solute at the atomic level (process 1) and the incorporation of solute-water van der Waals and electrostatic interactions (process 2). The angle-dependent integral equation theory combined with our morphometric approach and the three-dimensional interaction site model theory are applied to processes 1 and 2, respectively. Neither a stage of training nor parameterization is necessitated. For solutes with various sizes including proteins, the HFEs calculated by the new method are compared to those obtained using a molecular dynamics simulation based on solution theory in energy representation (the ER method developed by Matubayasi and co-workers), currently the most reliable tool. The agreement is very good especially for proteins. The new method is characterized by the following: The calculation can rapidly be finished; a solute possessing a significantly large total charge can be handled without difficulty; and since it yields not only the HFE but also its many physically insightful energetic and entropic components, it is best suited to the elucidation of mechanisms of diverse phenomena such as the receptor-ligand binding, different types of molecular recognition, and protein folding, denaturation, and association.

Performance of 3D-RISM-KH in Predicting Hydration Free Energy: Effect of Solute Parameters.

The three-dimensional reference interaction site model molecular solvation theory with the Kovalenko-Hirata closure relation has been shown to produce excellent solvation characteristics for a large class of (bio)chemical systems in solution. Correct calculation of hydration free energy is central to successful application of any solvation model. In order to find out the best possible force-field parameters to be used for hydration free energy calculation with the aforementioned theory, we have developed an extended database containing a large number of experimental solvation free energies available in the current literature and used a plethora of theoretical models for assessment. The general Amber force field was found to perform satisfactorily, whereas special care should be taken in solute charge assignment with the universal force field.

Evaluation of the SAFT-γ Mie force field with solvation free energy calculations

Solvation free energies can be essential in the process of evaluating and developing force fields. In addition, they can be used to obtain a diversity of other thermodynamic properties such as Henry's law constants and partition coefficients, for instance. In this paper, we study the solvation of asphaltene-like molecules (polyaromatic hydrocarbons) in both aqueous and organic solvents by using a coarse-grained model known as the SAFT-γ Mie force field. This model relies on a top-down approach in which the force-field parameters are obtained using an equation of state. The use of solvation free energy calculations to evaluate this force field can help in improving the model and increasing the scale in which these simulations can be applied. The results presented here were obtained by carrying out molecular dynamics simulations, with the expanded ensemble method being applied to sample an alchemical path. For this, we employed a softcore variant of the Mie potential. From the output of these simulations, we estimated free-energy differences by means of the Multistate Bennett Acceptance Ratio method. The results with organic solvents exhibited low absolute deviations from experimental data. In turn, hydration free energy calculations required a binary interaction parameter to be estimated from molecular dynamics data, which displays a flaw in the top-down parameterization approach. Fortunately, however, a single binary parameter could be used for all pairs of polyaromatic hydrocarbons with water, proving that the SAFT-γ Mie force field exhibits a suitable transferability property that deserves further investigation.