Free Energy Perturbation Method Research Articles

Rigorous Free Energy Calculations in Structure-Based Drug Design.

Structure-based drug design could benefit greatly from computational methodologies that accurately predict the binding affinity of small compounds to target biomolecules. However, the current scoring functions used to rank compounds in virtual screens by docking are not sufficiently accurate to guide reliably the design of tight binding ligands. Thus, there is strong interest in methodologies based on molecular simulations of protein-ligand complexes which are perceived to be more accurate and, with advances in computing power, amenable to routine use. This report provides an overview of the technical details necessary to understand, execute and analyze binding free energy calculations, using free energy perturbation or thermodynamic integration methods. Examples of possible applications in structure-based drug design are discussed. Current methodological limitations are highlighted as well as a number of ongoing developments to improve the scope, reliability, and practicalities of free energy calculations. These efforts are paving the way for a more common use of free energy calculations in molecular design.

Suitability of umbrella- and overlap-sampling methods for calculation of solid-phase free energies by molecular simulation

We examine the ability of two-stage free-energy perturbation methods to yield solid-phase free energies using a system of harmonically coupled particles as a reference. We consider two ways to construct a reference system, one based on derivatives of the intermolecular potential of the target system of interest (the conventional choice in lattice dynamics), and the other based on analysis of pairwise configurational correlations observed in simulations of the target system. For each case, we consider two perturbation techniques that compute the free energy difference between the target and reference systems while avoiding lengthy thermodynamic integration procedures. The methods are overlap sampling as optimized by Bennett, and umbrella sampling optimized in a similar fashion. Such methods require at most two simulations to yield a result, but they can fail if the target and reference do not share a sufficiently large set of relevant configurations. In particular, failure can be expected for large systems, and we examine the question of how large a system can be before this point is reached. Our test case is a face-centered cubic system of r(-12) soft spheres, and we find that for systems of up to 108 particles the methods are accurate for all temperatures up to melting; for systems of 256 particles the methods begin to break down at about half the melting temperature. Significantly, we observe that the correction to the harmonic reference is only weakly dependent on system size, suggesting an N-hybrid technique in which perturbation is applied to a small system and the result added to a large-system harmonic reference to obtain a good estimate of the correct large-system free energy. We also examine these approaches, along with thermodynamic integration in temperature, with respect to their computational efficiency. We find that Bennett's method using a derivative-based harmonic reference is the most efficient of all those examined, particularly when employed in the N-hybrid method.

Combining ab initio quantum mechanics with a dipole-field model to describe acid dissociation reactions in water: First-principles free energy and entropy calculations

We introduce a novel approach to compute dissociation free energy and entropy values in simulations that employ a density functional theory description of the acidic moiety and of the solvent. The approach consists of utilizing an alchemical transformation of a weak acid A-COOH into the strong acid B-COOH, which makes it practical to employ alchemical free energy perturbation methods in the context of ab initio molecular dynamics simulations. The present alchemical transformation circumvents the need to tackle changes in the total number of electrons and atoms by replacing the chemical residue responsible for the change in acidity with an easily tunable external effective potential. Our investigation demonstrates that (1) a simple but effective class of external potentials that control acidity changes in the acetic/trifluoroacetic acid series can be achieved by replacing the methyl and trifluoromethyl substituents by screened dipoles. Using this dipole-field/quantum-mechanics (DF/QM) approach one can predict gas-phase geometries, proton dissociation energies, total dipole moments, and water binding energies in good agreement with full-QM values. (2) The resulting alchemical perturbation calculations are stable and well converged and allow one to compute absolute pK(a) values whose accuracy is limited primarily by the exchange-correlation functional employed: H-COOH=2.5+/-0.6 (full-QM calculation), 3.7 (exp); F(3)C-COOH=0.4+/-0.6 (DF/QM calculation), 0.5 (exp); H(3)C-COOH=3.1+/-0.7 (DF/QM calculation), 4.7 (exp); 3) Our DF/QM model predicts that the difference in acidity between H-COOH and H(3)C-COOH is dominated by solvent entropy effects, in excellent agreement with experimental observations. The calculated difference between the dissociation energies of these acids is DeltaDelta(d)U=0.0+/-0.26 kcal/mol while the experimental value is 0.0+/-0.1 kcal/mol.

Acidity of the Aqueous Rutile TiO2(110) Surface from Density Functional Theory Based Molecular Dynamics.

The thermodynamics of protonation and deprotonation of the rutile TiO2(110) water interface is studied using a combination of density functional theory based molecular dynamics (DFTMD) and free energy perturbation methods. Acidity constants are computed from the free energy for chaperone assisted insertion/removal of protons in fully atomistic periodic model systems treating the solid and solvent at the same level of theory. The pKa values we find for the two active surface hydroxyl groups on TiO2(110), the bridge OH (Ti2OH(+)), and terminal H2O adsorbed on a 5-fold Ti site (TiOH2) are -1 and 9, leading to a point of zero proton charge of 4, well within the computational error margin (2 pKa units) from the experimental value (4.5-5.5). The computed intrinsic surface acidities have also been used to estimate the dissociation free energy of adsorbed water giving 0.6 eV, suggesting that water dissociation is unlikely on a perfect aqueous TiO2(110) surface. For further analysis, we compare to the predictions of the MUltiSIte Complexation (MUSIC) and Solvation, Bond strength, and Electrostatic (SBE) models. The conclusion regarding the MUSIC model is that, while there is good agreement for the acidity of an adsorbed water molecule, the proton affinity of the bridging oxygen obtained in the DFTMD calculation is significantly lower (more than 5 pKa units) than the MUSIC model value. Structural analysis shows that there are significant differences in hydrogen bonding, in particular to a bridging oxygen which is assumed to be stronger in the MUSIC model compared to what we find using DFTMD. Using DFTMD coordination numbers as input for the MUSIC model, however, led to a pKa prediction which is inconsistent with the estimates obtained from the DFTMD free energy calculation.

An Off-Lattice Flory-Huggins Approach of the Partitioning of Bulky Solutes between Polymers and Interacting Liquids

The desorption of additives and polymer residues from materials (P) is of significant concern for a wide range of applications, including polymer ageing, off-odors, and the safety of materials in contact with food or water. With diffusion coefficients, partition coefficients between P and hydrogen-bonding liquids (L), denoted Ki,L/P for a solute i, are fundamental quantities to assess the loss of plastics constituents. For a polymer with a crystallinity c, they are defined as ln(Ki,L/P) - ln(1-c) =(uexi,P- uexi,L)/(kBT), where {uexi,k}k=L,P are the excess chemical potentials and where kB is the Boltzmann constant. Our ambition was to relate calculations of {uexi,k}k=L,P at an atomistic scale to measurable partition coefficients obtained for bulky solutes with different stiffness and shape, such as hindered phenolic antioxidants, n-alkanes and n-alcohols. For large solutes in dense and cohesive phases (P or L), promising computation techniques involve free energy perturbation and non-equilibrium methods, thermodynamic integration, and extended ensembles (Gibbs or osmotic ensembles). Because the free energy landscape of the coupled i+k system contains several minima and highly heterogeneous barriers, a low convergence is generally achieved. To overcome such complications, we calculated {uexi,,k}k=L,P at atomistic scale in the framework of the Flory-Huggins approximation. The main advantage is that the mixing energy for a wide range of conformers was derived from pairwise interactions and a continuous approximation of the packing of molecules without representing explicitly large molecular systems. The translational entropy in L was inferred from a mesoscopic representation of the liquid medium matching the radial distributions calculated by isothermal and isobaric molecular dynamics simulations. The corresponding predictions were satisfactory compared with 38 experimental Ki,L/P values between polyethylenes, polystyrene and different alcohols (isopropanol, ethanol).

Predicting water uptake in poly(perfluorosulfonic acids) using force field simulation methods

Free energy perturbation methods were applied to predict water contents in hydrated poly(perfluorosulfonic acids) (PPFSA). The simulations were based on the TEAM force field which was derived from quantum mechanical data calculated for small molecules using density functional theory (DFT) and thermodynamic data of molecular liquids and crystal. The equilibrium water contents in three PPFSA polymers (Nafion-117, Nafion-115 and Hyflon) were predicted by evaluating excess chemical potentials of water in hydrated polymers and in pure water. High level of precision measured by average uncertainty of ca. 0.1 kcal mol(-1), and accuracy in terms of deviation from experimental data by ca. 0.2 kcal mol(-1) were obtained in the predicted excess chemical potentials. The predicted amounts of water uptake agree well with experimental values. In addition, the equilibrium and dynamic properties of hydrated Nafion-117 were calculated and the results agree well with the existing experimental and computational data. The entropy and enthalpy contributions in the calculated excess chemical potentials are analyzed and the results are consistent with intuition. A linear correlation between the entropies and enthalpies is identified for the systems studied, which indicates that just increasing the interaction energies between water and host materials does not guarantee enhancement of the water uptake.

Hydrophobic interactions in model enclosures from small to large length scales: non-additivity in explicit and implicit solvent models.

The binding affinities between a united-atom methane and various model hydrophobic enclosures were studied through high accuracy free energy perturbation methods (FEP). We investigated the non-additivity of the hydrophobic interaction in these systems, measured by the deviation of its binding affinity from that predicted by the pairwise additivity approximation. While only small non-additivity effects were previously reported in the interactions in methane trimers, we found large cooperative effects (as large as -1.14 kcal mol(-1) or approximately a 25% increase in the binding affinity) and anti-cooperative effects (as large as 0.45 kcal mol(-1)) for these model enclosed systems. Decomposition of the total potential of mean force (PMF) into increasing orders of multi-body interactions indicates that the contributions of the higher order multi-body interactions can be either positive or negative in different systems, and increasing the order of multi-body interactions considered did not necessarily improve the accuracy. A general correlation between the sign of the non-additivity effect and the curvature of the solute molecular surface was observed. We found that implicit solvent models based on the molecular surface area (MSA) performed much better, not only in predicting binding affinities, but also in predicting the non-additivity effects, compared with models based on the solvent accessible surface area (SASA), suggesting that MSA is a better descriptor of the curvature of the solutes. We also show how the non-additivity contribution changes as the hydrophobicity of the plate is decreased from the dewetting regime to the wetting regime.

Differential quantum tunneling contributions in nitroalkane oxidase catalyzed and the uncatalyzed proton transfer reaction

The proton transfer reaction between the substrate nitroethane and Asp-402 catalyzed by nitroalkane oxidase and the uncatalyzed process in water have been investigated using a path-integral free-energy perturbation method. Although the dominating effect in rate acceleration by the enzyme is the lowering of the quasiclassical free energy barrier, nuclear quantum effects also contribute to catalysis in nitroalkane oxidase. In particular, the overall nuclear quantum effects have greater contributions to lowering the classical barrier in the enzyme, and there is a larger difference in quantum effects between proton and deuteron transfer for the enzymatic reaction than that in water. Both experiment and computation show that primary KIEs are enhanced in the enzyme, and the computed Swain-Schaad exponent for the enzymatic reaction is exacerbated relative to that in the absence of the enzyme. In addition, the computed tunneling transmission coefficient is approximately three times greater for the enzyme reaction than the uncatalyzed reaction, and the origin of the difference may be attributed to a narrowing effect in the effective potentials for tunneling in the enzyme than that in aqueous solution.

Free energy perturbation and MM/PBSA studies on inclusion complexes of some structurally related compounds with β-cyclodextrin

Molecular dynamics (MD) simulations have been conducted to explore time-resolved guest–host interactions involving inclusion complex formation between β-cyclodextrin and organic molecules bearing two peripheral benzene rings in aqueous solution. Moreover, free energy perturbation (FEP) and thermodynamic integration (TI) methods at different simulation times have been employed to estimate the relative free energy of complexation. Also, the less computer-time demanding molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) method was used to estimate the free energy of complexation based on only 1-ns MD simulation. Results showed that both FEP and TI methods were able to reasonably reproduce the experimental thermodynamic quantities. However, long simulation times (e.g. 15 ns) were needed for benzoin mutating to benzanilide (BAN), while moderately shorter times were sufficient for BAN mutating to phenyl benzoate and for benzilic acid mutating to diphenylacetic acid. The results have been discussed in the light of the differences in the chemical structural and conformational features of the guest molecules. In general, it was apparent that the TI method requires less time for convergence of results than the FEP method. However, the less expensive MM/PBSA method proved capable of producing results that are in agreement with those of the more expensive TI and FEP methods.

Relative solvation free energies calculated using an ab initio QM/MM-based free energy perturbation method: dependence of results on simulation length

Molecular dynamics (MD) simulations in conjunction with thermodynamic perturbation approach was used to calculate relative solvation free energies of five pairs of small molecules, namely; (1) methanol to ethane, (2) acetone to acetamide, (3) phenol to benzene, (4) 1,1,1 trichloroethane to ethane, and (5) phenylalanine to isoleucine. Two studies were performed to evaluate the dependence of the convergence of these calculations on MD simulation length and starting configuration. In the first study, each transformation started from the same well-equilibrated configuration and the simulation length was varied from 230 to 2,540 ps. The results indicated that for transformations involving small structural changes, a simulation length of 860 ps is sufficient to obtain satisfactory convergence. In contrast, transformations involving relatively large structural changes, such as phenylalanine to isoleucine, require a significantly longer simulation length (>2,540 ps) to obtain satisfactory convergence. In the second study, the transformation was completed starting from three different configurations and using in each case 860 ps of MD simulation. The results from this study suggest that performing one long simulation may be better than averaging results from three different simulations using a shorter simulation length and three different starting configurations.

λ‐Dynamics free energy simulation methods

Free energy calculations are fundamental to obtaining accurate theoretical estimates of many important biological phenomena including hydration energies, protein-ligand binding affinities and energetics of conformational changes. Unlike traditional free energy perturbation and thermodynamic integration methods, lambda-dynamics treats the conventional "lambda" as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present article, we provide an overview of the theory of lambda-dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how lambda-dynamics has been used to rapidly and reliably compute relative hydration free energies and binding affinities for series of ligands, to accurately identify crystallographically observed binding modes starting from incorrect orientations, and to model the effects of mutations upon protein stability. Finally, we suggest how lambda-dynamics may be extended to facilitate modeling efforts in structure-based drug design.

Computations of Absolute Solvation Free Energies of Small Molecules Using Explicit and Implicit Solvent Model.

Accurate determination of absolute solvation free energy plays a critical role in numerous areas of biomolecular modeling and drug discovery. A quantitative representation of ligand and receptor desolvation, in particular, is an essential component of current docking and scoring methods. Furthermore, the partitioning of a drug between aqueous and nonpolar solvents is one of the important factors considered in pharmacokinetics. In this study, the absolute hydration free energy for a set of 239 neutral ligands spanning diverse chemical functional groups commonly found in drugs and drug-like candidates is calculated using the molecular dynamics free energy perturbation method (FEP/MD) with explicit water molecules, and compared to experimental data as well as its counterparts obtained using implicit solvent models. The hydration free energies are calculated from explicit solvent simulations using a staged FEP procedure permitting a separation of the total free energy into polar and nonpolar contributions. The nonpolar component is further decomposed into attractive (dispersive) and repulsive (cavity) components using the Weeks-Chandler-Anderson (WCA) separation scheme. To increase the computational efficiency, all of the FEP/MD simulations are generated using a mixed explicit/implicit solvent scheme with a relatively small number of explicit TIP3P water molecules, in which the influence of the remaining bulk is incorporated via the spherical solvent boundary potential (SSBP). The performances of two fixed-charge force fields designed for small organic molecules, the General Amber force field (GAFF), and the all-atom CHARMm-MSI, are compared. Because of the crucial role of electrostatics in solvation free energy, the results from various commonly used charge generation models based on the semiempirical (AM1-BCC) and QM calculations [charge fitting using ChelpG and RESP] are compared. In addition, the solvation free energies of the test set are also calculated using Poisson-Boltzmann (PB) and Generalized Born model of solvation (GB), which are two widely used continuum electrostatic implicit solvent models. The protocol for running the absolute solvation free energy calculations used throughout is automated as much as possible, with minimum user intervention, so that it can be used in large-scale analysis and force field optimization.

Hydrogenase/Ferredoxin Charge-Transfer Complexes: Effect of Hydrogenase Mutations on the Complex Association

The [FeFe]-hydrogenases in the green alga Chlamydomonas reinhardtii utilize photogenerated electrons to reduce protons into hydrogen gas. The electrons are supplied from photosystem I and transferred to the [FeFe]-hydrogenase through specific hydrogenase-ferredoxin association. To understand how structural and kinetic factors control the association better, we used Brownian dynamics simulation methods to simulate the charge-transfer complex formation between both native and in silico mutants of the [FeFe]-hydrogenase HYDA2 and the [2Fe2S]-ferredoxin FDX1 from C. reinhardtii . The changes in binding free energy between different HYDA2 mutants and the native FDX1 were calculated by the free-energy perturbation method. Within the limits of our current models, we found that two HYDA2 mutations, T99K(H) and D102K(H), led to lower binding free energies and higher association rate with FDX1 and are thus promising targets for improving hydrogen production rates in engineered organisms.

Atomistic Insights into Regulatory Mechanisms of the HER2 Tyrosine Kinase Domain: A Molecular Dynamics Study

HER2 (ErbB2/Neu) is a receptor tyrosine kinase belonging to the epidermal growth factor receptor (EGFR)/ErbB family and is overexpressed in 20–30% of human breast cancers. Although several crystal structures of ErbB kinases have been solved, the precise mechanism of HER2 activation remains unknown, and it has been suggested that HER2 is unique in its requirement for phosphorylation of Y877, a key tyrosine residue located in the activation loop. To elucidate mechanistic details of kinase domain regulation, we performed molecular dynamics simulations of a homology-modeled HER2 kinase structure in active and inactive conformations. Principal component analysis of the atomistic fluctuations reveals a tight coupling between the activation loop and catalytic loop that may contribute to alignment of residues required for catalysis in the active kinase. The free energy perturbation method is also employed to predict a role for phosphorylated Y877 in stabilizing the kinase conformations. Finally, simulation results are presented for a HER2/EGFR heterodimer and reveal that the dimeric interface induces a rearrangement of the αC helix toward the active conformation. Elucidation of the molecular regulatory mechanisms in HER2 will help establish structure-function relationships in the wild-type kinase, as well as predict mutations with a propensity for constitutive activation in HER2-mediated cancers.

Free Energy Calculations of Sparsomycin Analogs Binding to the Ribosome with Molecular Dynamics Simulations

The accurate calculation of absolute binding free energy is one of the holy grails of computer-aided drug design. The emerging successes reported in computing the binding free energy of small ligands to proteins using molecular dynamics (MD) simulations indicated that such physics-based approaches hold the promise of expediting the rational drug discovery process. Among numerous receptor-ligand systems, ribosome-antibiotic binding provides an important paradigm for studying the molecular recognition of RNAs by small molecules. The interactions of the 50S bacteria ribosomal subunit with antibiotic sparsomycin and its derivatives have been studied through the calculation of the binding free energy and the characterization of conformational dynamics. The standard binding free energies of the complexes were calculated using free energy perturbation (FEP) method. Restraining potentials affecting the orientational, translational and conformational freedom of the ligand and receptor were applied and then removed during the simulations to enhance the sampling and the convergence. The loss of ligand conformational entropy upon binding was estimated with Umbrella Sampling method by calculating the Potential of Mean Force as a function of the RMSD relative to the reference conformation of the ligand. Due to the large size of the ribosome, the Generalized Solvent Boundary Potential method was used to reduce the computational cost of MD/FEP calculations. For a deeply buried binding pocket in the ribosome, the fluctuation of solvent occupancy during the alchemical free energy calculation was also characterized by combining the MD with Grand Canonical Monte Carlo simulation. This computational study further revealed the mechanism of ribosome-antibiotics interactions and shed light on the design of ribosomal drugs. With the above stated developments, the evaluation of the binding free energies has become computationally more appealing for large systems.

Free-Energy Analysis of Solvation with the Method of Energy Representation

A new theory of solutions, the method of energy representation, is introduced by adopting the solute-solvent interaction energy as the coordinate of distribution functions. The density-functional theory is formulated over the energy coordinate, and an approximate functional for the solvation free energy is given in terms of energy distribution functions in the solution and reference solvent systems. The method of energy representation greatly expands the scope of solution theory and is amenable to supercritical fluid, flexible molecules with intramolecular degrees of freedom, inhomogeneous system, and quantum-mechanical/molecular-mechanical (QM/MM) system. Through the combination with molecular simulation, the functional for the solvation free energy is demonstrated to perform well for nonpolar, polar, and ionic solutes in water over a wide range of thermodynamics conditions, with drastic reduction of the computational demand compared to the standard free-energy perturbation and thermodynamic integration methods. As an application to inhomogeneous system involving flexible species, the molecular binding into micelle and membrane is analyzed by treating micelle and membrane as a mixed solvent system consisting of water and amphiphilic molecule.

QM/MM−PBSA Method To Estimate Free Energies for Reactions in Proteins

We have developed a method to estimate free energies of reactions in proteins, called QM/MM-PBSA. It estimates the internal energy of the reactive site by quantum mechanical (QM) calculations, whereas bonded, electrostatic, and van der Waals interactions with the surrounding protein are calculated at the molecular mechanics (MM) level. The electrostatic part of the solvation energy of the reactant and the product is estimated by solving the Poisson-Boltzmann (PB) equation, and the nonpolar part of the solvation energy is estimated from the change in solvent-accessible surface area (SA). Finally, the change in entropy is estimated from the vibrational frequencies. We test this method for five proton-transfer reactions in the active sites of [Ni,Fe] hydrogenase and copper nitrite reductase. We show that QM/MM-PBSA reproduces the results of a strict QM/MM free-energy perturbation method with a mean absolute deviation (MAD) of 8-10 kJ/mol if snapshots from molecular dynamics simulations are used and 4-14 kJ/mol if a single QM/MM structure is used. This is appreciably better than the original QM/MM results or if the QM energies are supplemented with a point-charge model, a self-consistent reaction field, or a PB model of the protein and the solvent, which give MADs of 22-36 kJ/mol for the same test set.

Free Energy Calculations on the Binding of Colchicine and Its Derivatives with the α/β-Tubulin Isoforms

Tubulin is the target for numerous small molecule ligands which alter microtubule dynamics leading to cell cycle arrest and apoptosis. Many of these ligands are currently used clinically for the treatment of several types of cancer, and they bind to one of three distinct binding sites within beta-tubulin (paclitaxel, vinca, and colchicine), all of which have been identified crystallographically. Unfortunately, serious side effects always accompany chemotherapy since these drugs bind to tubulin indiscriminately, leading to the death of both cancerous and healthy cells. However, the existence and distribution of divergent tubulin isoforms provide a platform upon which we may build novel chemotherapeutic drugs that can differentiate between different cell types and therefore reduce undesirable side effects. We report results of computational analysis that aims at predicting differences between the binding energies of a family of colchicine derivatives against 10 human alpha/beta-tubulin isoforms. Free energy perturbation method has been used in our calculations and the results provide a proof of principle by indicating significant differences both among the derivatives and between tubulin isoforms.

Response to “Comment on ‘Free energy simulations of single and double ion occupancy in gramicidin A’ ” [J. Chem. Phys. 128, 227101 (2008)

We respond to the criticism that one-dimensional (1D) construction of the potential of mean force (PMF) of ions in channels is flowed. Comparison of the 1D PMF results in the gramicidin A channel with independent free energy difference calculations obtained by using the free energy perturbation and thermodynamic integration methods shows complete agreement, thus providing a justifications for the 1D PMF approximation.

Accurate calculations of free-energy differences by the distribution method

We employ the strategy used in the successive umbrella sampling method [P. Virnau and M. Muller, J. Chem. Phys. 120, 10925 (2004)] to obtain the energy-difference distribution over its desired range. This is very helpful in calculating free-energy differences, where the source of the error is well recognized as the insufficient sampling over the relevant tail region in the energy-difference distribution. The distribution method proposed here employs the idea of restricting the sampling within an appropriate energy range, as was presented by Shing and Gubbins in their restricted umbrella sampling method [Mol. Phys. 46, 1109 (1982)]. We demonstrate the efficiency of the distribution method by calculating the free-energy difference of a model of harmonic oscillators where the systems exhibit nonoverlap features in their important phase spaces through the original Metropolis sampling. For this particular case, we show that the distribution method outperforms the free-energy perturbation method and even the Bennett's acceptance ratio method [J. Comput. Phys. 22, 245 (1976)] with the fastest convergence and the smallest relative errors. We further demonstrate the application of the distribution method with a simple point charge water model.