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

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

New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations

In a recent article (Lee, M. S.; Salsbury, F. R. Jr.; Brooks, C. L., III. J Chem Phys 2002, 116, 10606), we demonstrated that generalized Born (GB) theory provides a good approximation to Poisson electrostatic solvation energy calculations if one uses the same definitions of molecular volume for each. In this work, we present a new and improved analytic method for reproducing the Lee-Richards molecular volume, which is the most common volume definition for Poisson calculations. Overall, 1% errors are achieved for absolute solvation energies of a large set of proteins and relative solvation energies of protein conformations. We also introduce an accurate SASA approximation that uses the same machinery employed by our GB method and requires a small addition of computational cost. The combined methodology is shown to yield an efficient and accurate implicit solvent representation for simulations of biopolymers.     

A highly efficient hybrid method for calculating the hydration free energy of a protein

We develop a new method for calculating the hydration free energy (HFE) of a protein with any net charge. The polar part of the energetic component in the HFE is expressed as a linear combination of four geometric measures (GMs) of the protein structure and the generalized Born (GB) energy plus a constant. The other constituents in the HFE are expressed as linear combinations of the four GMs. The coefficients (including the constant) in the linear combinations are determined using the three‐dimensional reference interaction site model (3D‐RISM) theory applied to sufficiently many protein structures. Once the coefficients are determined, the HFE and its constituents of any other protein structure are obtained simply by calculating the four GMs and GB energy. Our method and the 3D‐RISM theory give perfectly correlated results. Nevertheless, the computation time required in our method is over four orders of magnitude shorter.     

Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking

In molecular docking, it is challenging to develop a scoring function that is accurate to conduct high-throughput screenings. Most scoring functions implemented in popular docking software packages were developed with many approximations for computational efficiency, which sacrifices the accuracy of prediction. With advanced technology and powerful computational hardware nowadays, it is feasible to use rigorous scoring functions, such as molecular mechanics/Poisson Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA) in molecular docking studies. Here, we systematically investigated the performance of MM/PBSA and MM/GBSA to identify the correct binding conformations and predict the binding free energies for 98 protein-ligand complexes. Comparison studies showed that MM/GBSA (69.4%) outperformed MM/PBSA (45.5%) and many popular scoring functions to identify the correct binding conformations. Moreover, we found that molecular dynamics simulations are necessary for some systems to identify the correct binding conformations. Based on our results, we proposed the guideline for MM/GBSA to predict the binding conformations. We then tested the performance of MM/GBSA and MM/PBSA to reproduce the binding free energies of the 98 protein-ligand complexes. The best prediction of MM/GBSA model with internal dielectric constant 2.0, produced a Spearman's correlation coefficient of 0.66, which is better than MM/PBSA (0.49) and almost all scoring functions used in molecular docking. In summary, MM/GBSA performs well for both binding pose predictions and binding free-energy estimations and is efficient to re-score the top-hit poses produced by other less-accurate scoring functions.     

FEW: A workflow tool for free energy calculations of ligand binding

In the later stages of drug design projects, accurately predicting relative binding affinities of chemically similar compounds to a biomolecular target is of utmost importance for making decisions based on the ranking of such compounds. So far, the extensive application of binding free energy approaches has been hampered by the complex and time‐consuming setup of such calculations. We introduce the free energy workflow (FEW) tool that facilitates setup and execution of binding free energy calculations with the AMBER suite for multiple ligands. FEW allows performing free energy calculations according to the implicit solvent molecular mechanics (MM‐PB(GB)SA), the linear interaction energy, and the thermodynamic integration approaches. We describe the tool's architecture and functionality and demonstrate in a show case study on Factor Xa inhibitors that the time needed for the preparation and analysis of free energy calculations is considerably reduced with FEW compared to a fully manual procedure. © 2013 Wiley Periodicals, Inc.     

Improved precision and efficiency of free energy calculations for small systems using lambda-scaled atomic masses and separating conformational and transformational sampling

We present results showing the importance of appropriate treatment of atomic masses in molecular dynamics (MD)-based single topology free-energy perturbations (FEPs) on small molecule systems. The reversibility of gas phase simulations is significantly improved by scaling the atomic mass of mutated atoms with the lambda variable normally used for the scaling of energy terms. Because this effect is less pronounced for solvated systems, it will not cancel in estimates of the relative hydration free energy difference. The advantage of mass scaling is demonstrated by a null mutation of ethane to ethane and the calculation of the relative hydration free energy difference between ethane and n-propane. Furthermore, it is found that the simulation time necessary for converged MD/FEPs is prohibitively large for relative hydration free energy calculations on cyclic alkanes. Therefore, we explore an alternative free energy pathway including strongly constrained conformations to improve convergence in FEP simulations of flexible molecules.     

The treatment of solvation by a generalized Born model and a self-consistent charge-density functional theory-based tight-binding method

We present a model to calculate the free energies of solvation of small organic compounds as well as large biomolecules. This model is based on a generalized Born (GB) model and a self-consistent charge-density functional theory-based tight-binding (SCC-DFTB) method with the nonelectrostatic contributions to the free energy of solvation modeled in terms of solvent-accessible surface areas (SA). The parametrization of the SCC-DFTB/GBSA model has been based on 60 neutral and six ionic molecules composed of H, C, N, O, and S, and spanning a wide range of chemical groups. Effective atomic radii as parameters have been obtained through Monte Carlo Simulated Annealing optimization in the parameter space to minimize the differences between the calculated and experimental free energies of solvation. The standard error in the free energies of solvation calculated by the final model is 1.11 kcal mol(-1). We also calculated the free energies of solvation for these molecules using a conductor-like screening model (COSMO) in combination with different levels of theory (AM1, SCC-DFTB, and B3LYP/6-31G*) and compared the results with SCC-DFTB/GBSA. To assess the efficiency of our model for large biomolecules, we calculated the free energy of solvation for a HIV protease-inhibitor complex containing 3,204 atoms using the SCC-DFTB/GBSA and the SCC-DFTB/COSMO models, separately. The computed relative free energies of solvation are comparable, while the SCC-DFTB/GBSA model is three to four times more efficient, in terms of computational cost.     

Solvent models for protein–ligand binding: Comparison of implicit solvent poisson and surface generalized born models with explicit solvent simulations

Solvent effects play a crucial role in mediating the interactions between proteins and their ligands. Implicit solvent models offer some advantages for modeling these interactions, but they have not been parameterized on such complex problems, and therefore, it is not clear how reliable they are. We have studied the binding of an octapeptide ligand to the murine MHC class I protein using both explicit solvent and implicit solvent models. The solvation free energy calculations are more than 10³ faster using the Surface Generalized Born implicit solvent model compared to FEP simulations with explicit solvent. For some of the electrostatic calculations needed to estimate the binding free energy, there is near quantitative agreement between the explicit and implicit solvent model results; overall, the qualitative trends in the binding predicted by the explicit solvent FEP simulations are reproduced by the implicit solvent model. With an appropriate choice of reference system based on the binding of the discharged ligand, electrostatic interactions are found to enhance the binding affinity because the favorable Coulomb interaction energy between the ligand and protein more than compensates for the unfavorable free energy cost of partially desolvating the ligand upon binding. Some of the effects of protein flexibility and thermal motions on charging the peptide in the solvated complex are also considered. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 591–607, 2001     

Calculation of hydration free energy for a solute with many atomic sites using the RISM theory: A robust and efficient algorithm

We have developed an algorithm for solving the reference interaction site model (RISM) equations for water near a solute molecule with many atomic sites (interaction sites). It is a hybrid of the Newton–Raphson and Picard methods and is judiciously constructed. Various considerations are given so that the computer time can be saved as much as possible. The robustness and high efficiency of the algorithm has been demonstrated for calculating hydration free energies of Met-enkephalin (a peptide with 75 sites) with different conformations. The Jacobian matrix is treated as part of the input data, and it has been found that the same matrix can be used for a considerably large set of different conformations of the solute molecule. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1320–1326, 1997     

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.     

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.     

Evaluating the relative free energy of hydration of new thrombin inhibitor candidates using the finite difference thermodynamic integration (FDTI) method

Thrombin is a serine protease involved in blood coagulation. Since thrombin inhibitors appear to be effective in the treatment and prevention of thrombotic and embolic disorders, considerable attention has been focused on the structure and interactions of the enzyme. In this work, we calculated the relative free energies of hydration of the new thrombin inhibitor candidates, p‐substituted derivatives of benzamidine, a well‐known noncovalent thrombin inhibitor. We used molecular dynamics and the finite difference thermodynamic integration (FDTI) algorithm within the Discover program of MSI. We have shown that the orthogonality problem that occurs in the calculation of intraperturbed‐group contributions to the free energy is treated adequately by the FDTI method. We have also shown that problems of singularity and convergence in free energy calculations can be properly solved using this method. To conclude, the calculated free energies of hydration gave the following order of solvation for the candidates: p‐(2‐oxo‐1‐propyl)benzamidine > p‐methylbenzamidine > p‐ethylbenzamidine > p‐(1‐propyl)benzamidine > benzamidine. © 2001 John Wiley & Sons, Inc. Int J Quantum Chem, 2001     

A new set of atomic radii for accurate estimation of solvation free energy by Poisson–Boltzmann solvent model

The Poisson–Boltzmann implicit solvent (PB) is widely used to estimate the solvation free energies of biomolecules in molecular simulations. An optimized set of atomic radii (PB radii) is an important parameter for PB calculations, which determines the distribution of dielectric constants around the solute. We here present new PB radii for the AMBER protein force field to accurately reproduce the solvation free energies obtained from explicit solvent simulations. The presented PB radii were optimized using results from explicit solvent simulations of the large systems. In addition, we discriminated PB radii for N‐ and C‐terminal residues from those for nonterminal residues. The performances using our PB radii showed high accuracy for the estimation of solvation free energies at the level of the molecular fragment. The obtained PB radii are effective for the detailed analysis of the solvation effects of biomolecules. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Effect of aminoacids on the surface free energy of nitrofurantoin

Changes in nitrofurantoin surface free energy components, Lifshitz-van der Waals, _s ^LW , electron donor, _s ^– , and electron acceptor, _s ⁺ , due to adsorption of the aminoacids: lysine, alanine and glutamic acid, were determined by means of the thin-layer wicking technique. It was found that the aminoacids slightly increase the _s ^– component already at 10^–4 M concentration. They reduce the _s ⁺ component practically to zero, and a very sharp increase of _s ^– was observed when the nitrofurantoin surface was precovered from the solutions at concentrations larger tan 10^–4 M. It is concluded that the adsorption of the aminoacids takes place via hydrogen bonding and electrostatic interactions between nitrofurantoin surface and aminoacid molecules. The increase in the _s ^– parameter is probably caused by the presence of carboxyl groups in the aminoacid molecules.     

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.     

Theoretical studies on the effects of methods and parameterization on the calculated free energy of hydration for small molecules

Free energies of hydration (FEH) have been computed for 13 neutral and nine ionic species as a difference of theoretically calculated Gibbs free energies in solution and in the gas phase. In‐solution calculations have been performed using both SCIPCM and PCM polarizable continuum models at the density functional theory (DFT)/B3LYP and ab initio Hartree–Fock levels with two basis sets (6‐31G* and 6‐311++G**). Good linear correlation has been obtained for calculated and experimental gas‐phase dipole moments, with an increase by ～30% upon solvation due to solute polarization. The geometry distortion in solution turns out to be small, whereas solute polarization energies are up to 3 kcal/mol for neutral molecules. Calculation of free energies of hydration with PCM provides a balanced set of values with 6‐31G* and 6‐311++G** basis sets for neutral molecules and ionic species, respectively. Explicit solvent calculations within Monte Carlo simulations applying free energy perturbation methods have been considered for 12 neutral molecules. Four different partial atomic charge sets have been studied, obtained by a fit to the gas‐phase and in‐solution molecular electrostatic potentials at in‐solution optimized geometries. Calculated FEH values depend on the charge set and the atom model used. Results indicate a preference for the all‐atom model and partial charges obtained by a fit to the molecular electrostatic potential of the solute computed at the SCIPCM/B3LYP/6‐31G* level. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

Quantitative determination of surface energy using atomic force microscopy: the case of hydrophobic/hydrophobic contact and hydrophilic/hydrophilic contact

Atomic force microscopy (AFM) has been used to determine the surface energy of chemically modified surfaces at a local scale. In order to achieve this aim, it was necessary to graft both the AFM tip and the substrate with the same chemical functional groups. Two different organothiols terminated either by hydrophilic or hydrophobic chemical functionalities were used. Grafting process classically reported shows that after UV/ozone treatment for 30 min, the tip is coated by thermal deposition with 4‐5‐nm‐thick titanium layer followed by a 30‐nm‐thick gold layer. Finally, the tip is grafted by organothiols. The thickness of the layer deposited on the tip is of the same order of magnitude as the tip radius. To avoid the use of Ti and to decrease the thickness of the gold layer, we have developed a new way of grafting by using organic molecules like (3‐mercaptopropyl)triethoxysilane (MPS) as a linkage agent. Then this way of grafting was checked. Finally, AFM force‐distance curves, between grafted tips and chemically modified surface, were carried out in contact mode. Calibration of the various parts of the apparatus and especially of the cantilever (spring constant and tip radius) is of major importance to reach quantitative data. Finally, by applying a suitable theory of contact, we were able to determine the surface energy of our system. Copyright © 2005 John Wiley & Sons, Ltd.     

The application of three approximate free energy calculations methods to structure based ligand design: Trypsin and its complex with inhibitors

Three approximate free energy calculation methods are examined and applied to an example ligand design problem. The first of the methods uses a single simulation to estimate the relative binding free energies for related ligands that are not simulated. The second method is similar, except that it uses only first derivatives of free energy with respect to atomic parameters (most often charge, van der Waals equilibrium distance, and van der Waals well depth) to calculate free energy differences. The last method PROFEC (Pictorial Representation of Free Energy Components), generates contour maps that show how binding free energy changes when additional particles are added near the ligand. These three methods are applied to a benzamidine/trypsin complex. They each reproduce the general trends in the binding free energies, indicating that they might be useful for suggesting how ligands could be modified to improve binding and, consequently, useful in structure-based drug design.     

Convergence and error estimation in free energy calculations using the weighted histogram analysis method

The weighted histogram analysis method (WHAM) has become the standard technique for the analysis of umbrella sampling simulations. In this article, we address the challenges (1) of obtaining fast and accurate solutions of the coupled nonlinear WHAM equations, (2) of quantifying the statistical errors of the resulting free energies, (3) of diagnosing possible systematic errors, and (4) of optimally allocating of the computational resources. Traditionally, the WHAM equations are solved by a fixed-point direct iteration method, despite poor convergence and possible numerical inaccuracies in the solutions. Here, we instead solve the mathematically equivalent problem of maximizing a target likelihood function, by using superlinear numerical optimization algorithms with a significantly faster convergence rate. To estimate the statistical errors in one-dimensional free energy profiles obtained from WHAM, we note that for densely spaced umbrella windows with harmonic biasing potentials, the WHAM free energy profile can be approximated by a coarse-grained free energy obtained by integrating the mean restraining forces. The statistical errors of the coarse-grained free energies can be estimated straightforwardly and then used for the WHAM results. A generalization to multidimensional WHAM is described. We also propose two simple statistical criteria to test the consistency between the histograms of adjacent umbrella windows, which help identify inadequate sampling and hysteresis in the degrees of freedom orthogonal to the reaction coordinate. Together, the estimates of the statistical errors and the diagnostics of inconsistencies in the potentials of mean force provide a basis for the efficient allocation of computational resources in free energy simulations.