Free Energy Calculation Methods Research Articles

Driven Metadynamics: Reconstructing Equilibrium Free Energies from Driven Adaptive-Bias Simulations

We present a novel free-energy calculation method that constructively integrates two distinct classes of nonequilibrium sampling techniques, namely, driven (e.g., steered molecular dynamics) and adaptive-bias (e.g., metadynamics) methods. By employing nonequilibrium work relations, we design a biasing protocol with an explicitly time- and history-dependent bias that uses on-the-fly work measurements to gradually flatten the free-energy surface. The asymptotic convergence of the method is discussed, and several relations are derived for free-energy reconstruction and error estimation. Isomerization reaction of an atomistic polyproline peptide model is used to numerically illustrate the superior efficiency and faster convergence of the method compared with its adaptive-bias and driven components in isolation.

Assessing the Accuracy and Performance of Implicit Solvent Models for Drug Molecules: Conformational Ensemble Approaches

The accuracy and performance of implicit solvent methods for solvation free energy calculations were assessed on a set of 20 neutral drug molecules. Molecular dynamics (MD) provided ensembles of conformations in water and water-saturated octanol. The solvation free energies were calculated by popular implicit solvent models based on quantum mechanical (QM) electronic densities (COSMO-RS, MST, SMD) as well as on molecular mechanical (MM) point-charge models (GB, PB). The performance of the implicit models was tested by a comparison with experimental water-octanol transfer free energies (ΔG(ow)) by using single- and multiconformation approaches. MD simulations revealed difficulties in a priori estimation of the flexibility features of the solutes from simple structural descriptors, such as the number of rotatable bonds. An increasing accuracy of the calculated ΔG(ow) was observed in the following order: GB1 ~ PB < GB7 ≪ MST < SMD ~ COSMO-RS with a clear distinction identified between MM- and QM-based models, although for the set excluding three largest molecules, the differences among COSMO-RS, MST, and SMD were negligible. It was shown that the single-conformation approach applied to crystal geometries provides a rather accurate estimate of ΔG(ow) for rigid molecules yet fails completely for the flexible ones. The multiconformation approaches improved the performance, but only when the deformation contribution was ignored. It was revealed that for large-scale calculations on small molecules a recent GB model, GB7, provided a reasonable accuracy/speed ratio. In conclusion, the study contributes to the understanding of solvation free energy calculations for physical and medicinal chemistry applications.

Temperature dependent effective potential method for accurate free energy calculations of solids

We have developed a thorough and accurate method of determining anharmonic free energies, the temperature dependent effective potential technique (TDEP). It is based on \emph{ab initio} molecular dynamics followed by a mapping onto a model Hamiltonian that describes the lattice dynamics. The formalism and the numerical aspects of the technique are described in details. A number of practical examples are given, and results are presented, which confirm the usefulness of TDEP within \emph{ab initio} and classical molecular dynamics frameworks. In particular, we examine from first-principles the behavior of force constants upon the dynamical stabilization of body centered phase of Zr, and show that they become more localized. We also calculate phase diagram for $^4$He modeled with the Aziz \emph{et al.} potential and obtain results which are in favorable agreement both with respect to experiment and established techniques.

Separated topologies—A method for relative binding free energy calculations using orientational restraints

Orientational restraints can improve the efficiency of alchemical free energy calculations, but they are not typically applied in relative binding calculations, which compute the affinity difference been two ligands. Here, we describe a new "separated topologies" method, which computes relative binding free energies using orientational restraints and which has several advantages over existing methods. While standard approaches maintain the initial and final ligand in a shared orientation, the separated topologies approach allows the initial and final ligands to have distinct orientations. This avoids a slowly converging reorientation step in the calculation. The separated topologies approach can also be applied to determine the relative free energies of multiple orientations of the same ligand. We illustrate the approach by calculating the relative binding free energies of two compounds to an engineered site in Cytochrome C Peroxidase.

How Does (E)-2-(Acetamidomethylene)succinate Bind to Its Hydrolase? From the Binding Process to the Final Result

The binding of (E)-2-(acetamidomethylene)succinate (E-2AMS) to E-2AMS hydrolase is crucial for biological function of the enzyme and the last step reaction of vitamin B6 biological degradation. In the present study, several molecular simulation methods, including molecular docking, conventional molecular dynamics (MD), steered MD (SMD), and free energy calculation methods, were properly integrated to investigate the detailed binding process of E-2AMS to its hydrolase and to assign the optimal enzyme-substrate complex conformation. It was demonstrated that the substrate binding conformation with trans-form amide bond is energetically preferred conformation, in which E-2AMS's pose not only ensures hydrogen bond formation of its amide oxygen atom with the vicinal oxyanion hole but also provides probability of the hydrophobic interaction between its methyl moiety and the related enzyme's hydrophobic cavity. Several key residues, Arg146, Arg167, Tyr168, Arg179, and Tyr259, orientate the E-2AMS's pose and stabilize its conformation in the active site via the hydrogen bond interaction with E-2AMS. Sequentially, the binding process of E-2AMS to E-2AMS hydrolase was studied by SMD simulation, which shows the surprising conformational reversal of E-2AMS. Several important intermediate structures and some significant residues were identified in the simulation. It is stressed that Arg146 and Arg167 are two pivotal residues responsible for the conformational reversal of E-2AMS in the binding or unbinding. Our research has shed light onto the full binding process of the substrate to E-2AMS hydrolase, which could provide more penetrating insight into the interaction of E-2AMS with the enzyme and would help in the further exploration on the catalysis mechanism.

Meld: Modeling Peptide-Protein Interactions

We have recently developed Meld, a framework to Model with limited data. This framework combines sparse information coming from different experimental, bioinformatics and even evolution sources into a physics based methodology. The physics are implemented through an atomistic force field. An improved version of a hybrid Hamiltonian/temperature replica exchange [1] procedure allows us to handle the limited data along with the physics in a fashion that obeys detail balance. After our initial use in the CASP experiment (critical assessment of structure prediction) to predict the structures of proteins, we are now interested in showing how this procedure can help us to correctly model cases in which small peptides are interacting with proteins. This procedure has the advantage of having both peptide and protein as flexible units. Input data is introduced in the calculation as restraints. A naïve approach in which all the restraints are imposed at the same time would be a failure. Therefore, we allow for errors in the data by enforcing only a fraction of those restraints. We allow the physics to decide which restraints are the most compatible with the system. In this way, the conformational space is greatly reduced, and in combination with 100x improvements in sampling efficiency coming from GPU, allows us to converge into possible solutions. After our calculation is done, clustering methods allow us to identify the candidate docking regions. We plan to combine the top clustering solutions with all atom free energy methods based on confinement techniques [2] to identify the most likely binding site. 1. Sugita, Y., A. Kitao, and Y. Okamoto, “Multidimensional replica-exchange method for free-energy calculations”. Journal of Chemical Physics, 113(15): p. 6042-6051 (2000). 2. Tyka, M., Clarke, R. and Sessions, R. An efficient, path-independent method for free-energy calculation, J. Phys. Chem. B,110, p. 17212-17220 (2006).

A truncated quasiharmonic method for free energy calculations and finite-temperature applications

Harmonic-based finite-temperature calculation methods play an important role in the study of thermodynamic properties of materials. In this study, we propose a truncated quasiharmonic (TQH) method to approximate the Helmholtz free energy by truncating the high-order terms of finite-temperature vibrational energy. To evaluate the efficacy of the TQH method against other established finite-temperature methods, i.e. the quasiharmonic (QH), the modified local harmonic (MLH) and the local quasiharmonic (LQH) methods, analysis of a homogeneous and vacancy-containing atomic system is performed with each method and compared. We found that the TQH method provides improved accuracy over the MLH and LQH methods for a system containing defects while requiring less computational time than the QH method to achieve convergence.

Mutation-Induced Loop Opening and Energetics for Binding of Tamiflu to Influenza N8 Neuraminidase

Tamiflu, also known as oseltamivir (OTV), binds to influenza A neuraminidase (H5N1) with very high affinity (0.32 nM). However, this inhibitor binds to other neuraminidases as well. In the present work, a systematic computational study is performed to investigate the mechanism underlying the binding of oseltamivir to N8 neuraminidase (NA) in "open" and "closed" conformations of the 150-loop through molecular dynamics simulations and the popular and well established molecular mechanics Poisson-Boltzmann (MM-PBSA) free energy calculation method. Whereas the closed conformation is stable for wild type N8, it transforms into the open conformation for the mutants Y252H, H274Y, and R292K, indicating that bound to oseltamivir these mutants are preferentially in the open conformation. Our calculations show that the binding of wild type oseltamivir to the closed conformation of N8 neuraminidase is energetically favored compared to the binding to the open conformation. We observe water mediated binding of oseltamivir to the N8 neuraminidase in both conformations which is not seen in the case of binding of the same drug to the H5N1 neuraminidase. The decomposition of the binding free energy reveals the mechanisms underlying the binding and changes in affinity due to mutations. Considering the mutant N8 variants in the open conformation adopted during the simulations, we observe a significant loss in the size of the total binding free energy for the N8(Y252H)-OTV, N8(H274Y)-OTV, and N8(R292K)-OTV complexes compared to N8(WT)-OTV, mainly due to the decrease in the size of the intermolecular electrostatic energy. For R292K, an unfavorable shift in the van der Waals interactions also contributes to the drug resistance. The mutations cause a significant expansion in the active site cavity, increasing its solvent accessible surface compared to the crystal structures of both the open and closed conformations. Our study underscores the need to consider dynamics in rationalizing the structure-function relationships of various antiviral inhibitor-NA complexes.

Efficient and Accurate Free Energy Calculations on Trypsin Inhibitors.

Several combinations of free energy calculation methods have been applied to determine the relative free energies of binding between eight para-substituted benzamidines in complex with the serine protease trypsin. With the aim to improve efficiency and maintain accuracy, the linear response approximation (LRA), linear interaction energy (LIE) and third power fitting (TPF) are combined with the one-step perturbation (OSP) to determine the polar and apolar contributions to the free energy, respectively. It is shown that the combination TPF/OSP gives the most accurate results and is 4.5 times more efficient than the rigorous thermodynamic integration (TI). By projecting the electrostatic preorganization energy from the OSP simulations, an increase in efficiency of a factor 7.5 can even be achieved. Loss of accuracy with respect to the TI data is limited to 3.9 and 5.6 kJ/mol, respectively, making it an attractive approach for lead optimization programs in drug research.

GMCT : A Monte Carlo simulation package for macromolecular receptors

Generalized Monte Carlo titration (GMCT) is a versatile suite of computer programs for the efficient simulation of complex macromolecular receptor systems as for example proteins. The computational model of the system is based on a microstate description of the receptor and an average description of its surroundings in terms of chemical potentials. The receptor can be modeled in great detail including conformational flexibility and many binding sites with multiple different forms that can bind different ligand types. Membrane embedded systems can be modeled including electrochemical potential gradients. Overall properties of the receptor as well as properties of individual sites can be studied with a variety of different Monte Carlo (MC) simulation methods. Metropolis MC, Wang-Landau MC and efficient free energy calculation methods are included. GMCT is distributed as free open source software at www.bisb.uni-bayreuth.de under the terms of the GNU Affero General Public License.

Free Energy Calculation Methods of Biomolecules

Free Energy Calculation Methods of Biomolecules

Efficient calculation of α- and β-nitrogen free energies and coexistence conditions via overlap sampling with targeted perturbation

A recently introduced solid-phase free-energy calculation method that is based upon overlap sampling with targeted free-energy perturbation is further developed and extended to systems with orientational degrees of freedom. Specifically we calculate the absolute free energy of the linear-molecular nitrogen model of Etter et al., examining both the low-temperature low-pressure α-N(2) structure and the orientationally disordered β-N(2) phase. In each perturbation (for the α-N(2) phase) to determine the free-energy difference between systems at adjacent temperatures, harmonic coordinate scaling is applied to both the translational and rotational degrees of freedom in the nitrogen molecule to increase the phase-space overlap of the two perturbing systems and consequently, improve the free-energy difference results. For the plastic β-N(2) phase, a novel method that requires several perturbation paths is introduced to calculate its absolute free energy. Through these methods, the absolute free energies for both the α-N(2) and β-N(2) phase can be accurately and precisely determined. We find again that the anharmonic contribution to the free energy has weak dependence on system size. The transition properties for the α-N(2) and β-N(2) phase are also investigated. The α-β phase transition for the model at atmospheric pressure (0.1 MPa) is found to occur at 40.35 ± 0.01 K with volumetric and entropy changes of 0.44 ± 0.01 cm(3)/mol and 1.99 ± 0.01 cal/mol.K respectively.

Analyzing the robustness of the MM/PBSA free energy calculation method: Application to DNA conformational transitions

The ability to predict and characterize free energy differences associated with conformational equilibria or the binding of biomolecules is vital to understanding the molecular basis of many important biological functions. As biological studies focus on larger molecular complexes and properties of the genome, proteome, and interactome, the development and characterization of efficient methods for calculating free energy becomes increasingly essential. The aim of this study is to examine the robustness of the end-point free energy method termed the molecular mechanics Poisson-Boltzmann solvent accessible surface area (MM/PBSA) method. Specifically, applications of MM/PBSA to the conformational equilibria of nucleic acid (NA) systems are explored. This is achieved by comparing A to B form DNA conformational free energy differences calculated using MM/PBSA with corresponding free energy differences determined with a more rigorous and time-consuming umbrella sampling algorithm. In addition, the robustness of NA MM/PBSA calculations is also evaluated in terms of the sensitivity towards the choice of force field and the choice of solvent model used during conformational sampling. MM/PBSA calculations of the free energy difference between A-form and B-form DNA are shown to be in very close agreement with the PMF result determined using an umbrella sampling approach. Further, it is found that the MM/PBSA conformational free energy differences were also in agreement using either the CHARMM or AMBER force field. The influence of ionic strength on conformational stability was particularly insensitive to the choice of force field. Finally, it is also shown that the use of a generalized Born implicit solvent during conformational sampling results in free energy estimates that deviate slightly from those obtained using explicitly solvated MD simulations in these NA systems.

Mechanism of Ser88 Phosphorylation-Induced Dimer Dissociation in Dynein Light Chain LC8

Dynein light chain LC8 is a highly conserved, dimeric protein involved in a variety of essential cellular events. Phosphorylation at Ser88 was found to promote mammalian cell survival and regulate the dimer to monomer transition at physiological pH. Combining molecular dynamics (MD) simulation and free energy calculation methods, we explored the atomistic mechanism of the phosphorylation-induced dimer dissociation. The MD simulation revealed that phosphorylation/phosphomimetic mutation at Ser88 opens an entrance into the dimer interface for water molecules, which disturb the hydrogen bond network around His55 and is expected to raise the pK(a) value and protonation ratio of His55 as well. The free energy calculations showed that the S88E mutation destabilized the dimer by 6.6 kcal/mol, in good agreement with the experimental value of 8.1 kcal/mol. The calculated destabilization upon phosphorylation is 50.8 kcal/mol, showing that phosphorylation definitely prevents dimer formation under physiological conditions. Further analysis of the calculated free energy changes demonstrated that the electrostatic contribution dominates the impact of phosphorylation on dimer dissociation.

Free energy calculation using molecular dynamics simulation combined with the three dimensional reference interaction site model theory. I. Free energy perturbation and thermodynamic integration along a coupling parameter

This article proposes a free energy calculation method based on the molecular dynamics simulation combined with the three dimensional reference interaction site model theory. This study employs the free energy perturbation (FEP) and the thermodynamic integration (TDI) along the coupling parameters to control the interaction potential. To illustrate the method, we applied it to a complex formation process in aqueous solutions between a crown ether molecule 18-Crown-6 (18C6) and a potassium ion as one of the simplest model systems. Two coupling parameters were introduced to switch the Lennard-Jones potential and the Coulomb potential separately. We tested two coupling procedures: one is a "sequential-coupling" to couple the Lennard-Jones interaction followed by the Coulomb coupling, and the other is a "mixed-coupling" to couple both the Lennard-Jones and the Coulomb interactions together as much as possible. The sequential-coupling both for FEP and TDI turned out to be accurate and easily handled since it was numerically well-behaved. Furthermore, it was found that the sequential-coupling had relatively small statistical errors. TDI along the mixed-coupling integral path was to be carried out carefully, paying attention to a numerical behavior of the integrand. The present model system exhibited a nonmonotonic behavior in the integrands for TDI along the mixed-coupling integral path and also showed a relatively large statistical error. A coincidence within a statistical error was obtained among the results of the free energy differences evaluated by FEP, TDI with the sequential-coupling, and TDI with the mixed-coupling. The last one is most attractive in terms of the computer power and is accurate enough if one uses a proper set of windows, taking the numerical behavior of the integrands into account. TDI along the sequential-coupling integral path would be the most convenient among the methods we tested, since it seemed to be well-balanced between the computational load and the accuracy. The numerical results reported in this article qualitatively agree with the experimental data for the potassium ion recognition by the 18C6 in aqueous solution.

Metadynamics As a Tool for Mapping the Conformational and Free-Energy Space of Peptides — The Alanine Dipeptide Case Study

There is a need for a fast, accurate, and reliable method of sampling the conformational space of peptides and proteins in order to obtain a balanced free-energy profile which can lead to our understanding of protein structure. We have utilized metadynamics for the conformational study of the solvated alanine dipeptide molecule, and our results show that the method has proven to be competent as a fast, robust, and reliable method for the conformation free-energy calculations of peptides in an explicit solvent, surpassing traditional methods such as umbrella sampling. We have also addressed the issue of the influence of different water models on the resulting free-energy profile in order to consistently decompose the setting of our simulation. All of the explicit water models for the simulation of biomolecules TIP3P, TIP4P, TIP4P/Ew, TIP5P, and SPCE have exhibited similar effects on the conformational preferences of alanine dipeptide with no significant differences. On the other hand, by comparing the potential energy surface in the gas phase and the free-energy surface in a water environment, we have shown that the interaction with water molecules is one of the most important structure-driving elements, with a great influence on the free-energy surface (FES) of the solvated peptide and the conformational preferences of the peptide backbone. All of the tested force fields (ff03, ff99SB, opls-aa, and charmm27) appreciably differ in the population of the individual conformers and the barriers between them. Significant divergence was found on both the potential energy surface (PES) as well as free-energy surface (FES) calculated by charmm27. We have therefore concluded that the differences originate dominantly from the parametrization of the peptide backbone in the given force field rather than from a noncovalent interaction with water molecules.

Free energy estimation of short DNA duplex hybridizations

BackgroundEstimation of DNA duplex hybridization free energy is widely used for predicting cross-hybridizations in DNA computing and microarray experiments. A number of software programs based on different methods and parametrizations are available for the theoretical estimation of duplex free energies. However, significant differences in free energy values are sometimes observed among estimations obtained with various methods, thus being difficult to decide what value is the accurate one.ResultsWe present in this study a quantitative comparison of the similarities and differences among four published DNA/DNA duplex free energy calculation methods and an extended Nearest-Neighbour Model for perfect matches based on triplet interactions. The comparison was performed on a benchmark data set with 695 pairs of short oligos that we collected and manually curated from 29 publications. Sequence lengths range from 4 to 30 nucleotides and span a large GC-content percentage range. For perfect matches, we propose an extension of the Nearest-Neighbour Model that matches or exceeds the performance of the existing ones, both in terms of correlations and root mean squared errors. The proposed model was trained on experimental data with temperature, sodium and sequence concentration characteristics that span a wide range of values, thus conferring the model a higher power of generalization when used for free energy estimations of DNA duplexes under non-standard experimental conditions.ConclusionsBased on our preliminary results, we conclude that no statistically significant differences exist among free energy approximations obtained with 4 publicly available and widely used programs, when benchmarked against a collection of 695 pairs of short oligos collected and curated by the authors of this work based on 29 publications. The extended Nearest-Neighbour Model based on triplet interactions presented in this work is capable of performing accurate estimations of free energies for perfect match duplexes under both standard and non-standard experimental conditions and may serve as a baseline for further developments in this area of research.

Analysis of α4 β1integrin specific antagonists binding modes: structural insights by molecular docking, molecular dynamics and linear interaction energy method for free energy calculations

The VLA-4 antigen (α4β1 integrin) is involved in the pathophysiology of a variety of diseases including asthma, multiple sclerosis, rheumatoid arthritis and diabetes. The ligand selectivity toward this integrin remains a difficult problem, mainly due to the fact that 3D structures of most integrins are still unknown. We initially built a 3D computational model of the α4β1 ligand binding site, taking the crystal structure of the integrin αVβ3 as template. Then, we performed a computational study on a set of seven α4β1 antagonists, evaluating the binding modes of 4-[N'-(2-methylphenyl)ureido]phenylacetyl and derivatives by molecular docking. Molecular dynamics simulations were used to improve the receptor-ligand energy landscape exploration by the docking algorithm. The compounds were systematically arranged in two main binding modes, and in all cases, pointed out that these antagonists preferably bind to the α4β1 integrin active site in an extended conformation that resembles the one in solution. LIE (linear interaction energy) calculations also confirmed this statement given that the most prevailing binding mode is also the energetically most favored one. This study benefits the comprehension of the mechanism of this family of antagonists and may provide useful information for rational drug design.

3P070 Multi-scale free energy calculation method with the combination of coarse-grained and all-atom models.(Protein: Property,The 48th Annual Meeting of the Biophysical Society of Japan)

3P070 Multi-scale free energy calculation method with the combination of coarse-grained and all-atom models.(Protein: Property,The 48th Annual Meeting of the Biophysical Society of Japan)

Free energy of a⟨110⟩dumbbell interstitial defect in bcc Fe: Harmonic and anharmonic contributions

The stability of interstitial defect and dislocation structures in bcc Fe as a function of temperature is believed to play a crucial role in determining defect evolution under irradiation conditions. The vibrational properties of defects constitute one contribution to the corresponding energetics and much work has been done within the harmonic approximation to determine the vibrational formation free energy and formation entropy of such defects. Defects can however exhibit strong local strain fields that break the cubic symmetry of the bcc lattice leading to large anharmonicities and a breakdown of the harmonic approximation as an accurate means to calculate vibrational thermodynamic quantities. Moreover, if defect diffusion is active at a time scale comparable to an atomic vibration, strong anharmonicities will always be present at any finite temperature. The current work investigates the vibrational free energy and entropy of the $⟨110⟩$ self-interstitial dumbbell defect in bcc Fe using both harmonic and anharmonic free-energy calculation methods for a range of modern empirical potentials. It is found that depending on the empirical potential and for temperatures where diffusion is limited, the harmonic approximation is justified especially for empirical potentials that have been fitted to third-order elastic constants. The unique applicability range of such calculations for bcc Fe is also discussed given that with rising temperature spin fluctuations become increasingly important ultimately leading to a softening of the 110 shear modulus and to the $\ensuremath{\alpha}\text{-bcc}/\ensuremath{\gamma}\text{-bcc}$ structural phase transformation.