Search

Acid dissociation, and thus liberation of excess protons in small water droplets, impacts on diverse fields such as interstellar, atmospheric or environmental chemistry. At cryogenic temperatures below 1 K, it is now well established that as few as four water molecules suffice to dissociate the generic strong acid HCl, yet temperature‐driven recombination sets in simply upon heating that cluster. Here, the fundamental question is posed of how many more water molecules are required to stabilize a hydrated excess proton at room temperature. Ab initio path integral simulations disclose that not five, but six water molecules are needed at 300 K to allow for HCl dissociation independently from nuclear quantum effects. In order to provide the molecular underpinnings of these observations, the classical and quantum free energy profiles were decomposed along the dissociation coordinate in terms of the corresponding internal energy and entropy profiles. What decides in the end about acid dissociation, and thus ion pair formation, in a specific microsolvated water cluster at room temperature is found to be a fierce competition between classical configurational entropy and internal energy, where the former stabilizes the undissociated state whereas the latter favors dissociation. It is expected that these are generic findings with broad implications on acid–base chemistry depending on temperature in small water assemblies., Temperature dependent: Extensive ab initio path integral simulations disclose that the number of water molecules required to allow for acid dissociation, and thus for the generation of hydrated excess protons, in restricted aqueous environments such as finite HCl(H2O)n clusters, is temperature dependent.

Aquation free energy profiles of neutral cisplatin and cationic monofunctional derivatives, including triaminochloroplatinum(II) and cis-diammine(pyridine)chloroplatinum(II), were computed using state of the art thermodynamic integration, for which temperature and solvent were accounted for explicitly using density functional theory-based canonical molecular dynamics (DFT-MD). For all the systems, the "inverse-hydration" where the metal center acts as an acceptor of hydrogen bond has been observed. This has motivated to consider the inversely bonded solvent molecule in the definition of the reaction coordinate required to initiate the constrained DFT-MD trajectories. We found that there exists little difference in free enthalpies of activation, such that these platinum-based anticancer drugs are likely to behave the same way in aqueous media. Detailed analysis of the microsolvation structure of the square-planar complexes, along with the key steps of the aquation mechanism, is discussed.

Peptide-based therapeutics are an alternative to small molecule drugs as they offer superior specificity, lower toxicity, and easy synthesis. Here we present an approach that leverages the dramatic performance increase afforded by the recent arrival of GPU accelerated thermodynamic integration (TI). GPU TI facilitates very fast, highly accurate binding affinity optimization of peptides against therapeutic targets. We benchmarked TI predictions using published peptide binding optimization studies. Prediction of mutations involving charged side-chains was found to be less accurate than for non-charged, and use of a more complex 3-step TI protocol was found to boost accuracy in these cases. Using the 3-step protocol for non-charged side-chains either had no effect or was detrimental. We use the benchmarked pipeline to optimize a peptide binding to our recently discovered cancer target: EME1. TI calculations predict beneficial mutations using both canonical and non-canonical amino acids. We validate these predictions using fluorescence polarization and confirm that binding affinity is increased. We further demonstrate that this increase translates to a significant reduction in pancreatic cancer cell viability.

We present a molecular dynamics simulation study of alkali metal cation transport through the double-helical and the head-to-head conformers of the gramicidin ion channel. Our approach is based on a thermodynamic integration network, which consists of a sequence of transport reactions, absolute free energies of solvation and cycles of alchemical transmutations of the ions. In this manner, we can reliably estimate free energies and their statistical errors via a least-squares method without imposing external forces on the system. Within the double helical channel, we find a free energy surface typical for hopping transport between isoenergetic sites of ion localization, separated by comparatively large activation barriers. For fast transport through the head-to-head conformation, the thermodynamic network scheme starts to break down. © 2018 Wiley Periodicals, Inc.

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD-based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907-atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.

Free-energy simulation (FES) is widely used in science and engineering. Unconstrained FES (UFES, e. g., umbrella sampling with histogram binning) and constrained FES (CFES, e. g., blue-moon sampling with mean-force/thermodynamic integration) are two different types of FES. People prefer UFES to CFES, e. g., people correct constrained free-energy profile (CFEP) to unconstrained FEP (UFEP), not the other way around. But remarkably, herein we show for a 1D free body subject to zero force, from an UFEP we find an infinitely-high energy barrier. By contrast, we find no energy barrier from any CFEP. Besides, by providing in-depth comparisons between UFEP and CFEP from the perspective of Jacobian scale factor, this work may also partially address recent issues on FES: (1) UFEP is preferred, but why do three recent papers conclude we should use CFEP for the transition-state theory and why are these three papers not fully consistent with one another? (2) For relating UFEP to CFEP, why are there two different versions of Fixman term associated with velocity and momentum constraints, respectively? (3) How to normalize CFEP, for which the equation seems to be even unavailable?

GSK3β kinase is a noteworthy target for discovery of the drugs that will be used to treat several diseases. In the effort to identify a new inhibitor lead compound, we utilized thermodynamic integration (TI)-molecular dynamics (MD) simulation and kinase assay to investigate the bindings between GSK3β kinase and five compounds that were analogous to a known inhibitor with an available crystal structure. TI-MD simulations of the first two compounds (analogs 1 and 2) were used for calibration. The computed binding affinities of analogs 1 and 2 agreed well with the experimental results. The rest three compounds (analogs 3-5) were newly obtained from a database search, and their affinity data were newly measured in our labs. TI-MD simulations predicted the binding modes and the computed ΔΔG values have a reasonably good correlation with the experimental affinity data. These newly identified inhibitors appear to be new leads according to our survey of GSK3β inhibitors listed in recent review articles. The predicted binding modes of these compounds should aid in designing new derivatives of these compounds in the future.

Many biomolecules undergo conformational changes associated with allostery or ligand binding. Observing these changes in computer simulations is difficult if their timescales are long. These calculations can be accelerated by observing the transition on an auxiliary free energy surface with a simpler Hamiltonian and connecting this free energy surface to the target free energy surface with free energy calculations. Here, we show that the free energy legs of the cycle can be replaced with energy representation (ER) density functional approximations. We compute: (1) The conformational free energy changes for alanine dipeptide transitioning from the right-handed free energy basin to the left-handed basin and (2) the free energy difference between the open and closed conformations of β-cyclodextrin, a “host” molecule that serves as a model for molecular recognition in host-guest binding. β-cyclodextrin contains 147 atoms compared to 22 atoms for alanine dipeptide, making β-cyclodextrin a large molecule for which to compute solvation free energies by free energy perturbation or integration methods and the largest system for which the ER method has been compared to exact free energy methods. The ER method replaced the 28 simulations to compute each coupling free energy with two endpoint simulations, reducing the computational time for the alanine dipeptide calculation by about 70% and for the β-cyclodextrin by > 95%. The method works even when the distribution of conformations on the auxiliary free energy surface differs substantially from that on the target free energy surface, although some degree of overlap between the two surfaces is required. © 2016 Wiley Periodicals, Inc.

Aim of this contribution is to review some recent quantum mechanical approaches used to compute the redox potentials of transition metal complexes, with the emphasis on copper and iron species, which are particularly relevant in inorganic biochemistry and in synthetic chemistry of bio-mimetic compounds. The paper presents also new DFT results obtained on Cu and Fe aquo ions in the framework of the Thermodynamic Integration and Grand Canonical Ensemble approaches. Such results show that without explicit inclusion of water molecules in the external solvation shells (even using a continuum solvation model) also very advanced methodologies fail to predict the redox potential in an acceptable manner. This is a confirmation of some previous studies which however never addressed this specific problem along the aforementioned approaches. Better results are obtained, on the contrary, on a series of Cu(II) complexes with Gly, Ala, en, Im, and water ligands with coordination type N2O2 or N4. In this case, the complexes are surrounded by nine water molecules which may partially alleviate the inadequacy of the continuum solvent models, especially in the case of highly positively charged species.

The evaluation of water binding free energies around solute molecules is important for the thermodynamic characterization of hydration or association processes. Here, a rapid approximate method to estimate water binding free energies around (bio)macromolecules from a single molecular dynamics simulation is presented. The basic idea is that endpoint free-energy calculation methods are applied and the endpoint quantities are monitored on a three-dimensional grid around the solute. Thus, a gridded map of water binding free energies around the solute is obtained, that is, from a single short simulation, a map of favorable and unfavorable water binding sites can be constructed. Among the employed free-energy calculation methods, approaches involving endpoint information pertaining to actual thermodynamic integration calculations or endpoint information as exploited in the linear interaction energy method were examined. The accuracy of the approximate approaches was evaluated on the hydration of a cage-like molecule representing either a nonpolar, polar, or charged water binding site and on α- and β-cyclodextrin molecules. Among the tested approaches, the linear interaction energy method is considered the most viable approach. Applying the linear interaction energy method on the grid around the solute, a semi-quantitative thermodynamic characterization of hydration around the whole solute is obtained. Disadvantages are the approximate nature of the method and a limited flexibility of the solute. © 2016 Wiley Periodicals, Inc.

Thermodynamic integration (TI) can provide accurate binding free energy insights in a lead optimization program, but its high computational expense has limited its usage. In the effort of developing an efficient and accurate TI protocol for FabI inhibitors lead optimization program, we carefully compared TI with different Amber molecular dynamics (MD) engines (sander and pmemd), MD simulation lengths, the number of intermediate states and transformation steps, and the Lennard-Jones and Coulomb Softcore potentials parameters in the one-step TI, using eleven benzimidazole inhibitors in complex with Francisella tularensis enoyl acyl reductase (FtFabI). To our knowledge, this is the first study to extensively test the new AMBER MD engine, pmemd, on TI and compare the parameters of the Softcore potentials in the one-step TI in a protein-ligand binding system. The best performing model, the one-step pmemd TI, using 6 intermediate states and 1 ns MD simulations, provides better agreement with experimental results (RMSD = 0.52 kcal/mol) than the best performing implicit solvent method, QM/MM-GBSA from our previous study (RMSD = 3.00 kcal/mol), while maintaining similar efficiency. Briefly, we show the optimized TI protocol to be highly accurate and affordable for the FtFabI system. This approach can be implemented in a larger scale benzimidazole scaffold lead optimization against FtFabI. Lastly, the TI results here also provide structure-activity relationship insights, and suggest the parahalogen in benzimidazole compounds might form a weak halogen bond with FabI, which is a well-known halogen bond favoring enzyme.

In a recent article (Bieler et al., J. Chem. Theory Comput. 2014, 10, 3006), we introduced a combination of λ-dynamics and local-elevation umbrella-sampling termed λ-LEUS to calculate free-energy changes associated with alchemical processes using molecular dynamics simulations. This method was suggested to be more efficient than thermodynamic integration (TI), because the dynamical variation of the alchemical variable λ opens up pathways to circumvent barriers in the orthogonal space (defined by the N – 1 degrees of freedom that are not subjected to the sampling enhancement), a feature λ-LEUS shares with Hamiltonian replica-exchange (HR) approaches. However, the mutation considered, hydroquinone to benzene in water, was no real challenge in terms of orthogonal-space properties, which were restricted to solvent-relaxation processes. In the present article, we revisit the comparison between TI and λ-LEUS considering non-trivial mutations of the central residue X of a KXK tripeptide in water (with X = G, E, K, S, F, or Y). Side-chain interactions that may include salt bridges, hydrogen bonds or steric clashes lead to slow relaxation in the orthogonal space, mainly in the two-dimensional subspace spanned by the central φ and ψ dihedral angles of the peptide. The efficiency enhancement afforded by λ-LEUS is confirmed in this more complex test system and can be attributed explicitly to the improved sampling of the orthogonal space. The sensitivity of the results to the nontrivial choices of a mass parameter and of a thermostat coupling time for the alchemical variable is also investigated, resulting in recommended ranges of 50 to 100 u nm2 and 0.2 to 0.5 ps, respectively. © 2015 Wiley Periodicals, Inc.

We show that the solvation free energy and vapor pressure, important thermodynamic properties of pure substances in liquid or solid states, can be obtained from short, about 20 ps, molecular dynamics simulations. The method combines the determination of free energy of a chemical in vacuum using the normal-mode analysis (NMA, energy minimization), and in the condensed phase using the two-phase thermodynamic (2PT) model. We have examined the calculation results for common liquids and solids, including water, alcohol, acid, aromatics, and alkanes. The results, referred to as 2PT-NMA, is comparable to those calculated from thermodynamic integration (TI) for liquids, and is readily applicable to solids, where simple TI is not applicable. Furthermore, the free energy from 2PT-NMA converges (20 ps) much faster than that from TI (1 ns). The new method could be a very useful tool for fast screening of condensed phase pressure from the trajectory of MD simulations. © 2015 American Institute of Chemical Engineers AIChE J, 61: 2298–2306, 2015

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein-specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild-type azurin (ΔEcal) were calculated and compared with experimental values. The ΔEcal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein-specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations.

The amyloid pathology associated with long-term haemodialysis is due to the deposition of β2-microglobulin, the non-polymorphic light chain of class I major histocompatibility complex, that accumulates at bone joints into amyloid fibrils. Several lines of evidence show the relevance of the tryptophan residue at position 60 for the fibrillogenic transition of the protein. A comparative 15N NMR relaxation analysis is presented for wild-type human β2-microglobulin and W60G β2-microglobulin, i.e. the mutant with a glycyne replacing the natural tryptophan residue at position 60. The experimental data, collected at 11.4 T and 310 K, were analyzed by means of the reduced spectral density approach. Molecular dynamics (MD) simulations and corresponding thermodynamic integration, together with hydrodynamic calculations were performed to support data interpretation. The analysis results for the mutant protein are consistent with a reduced aggregation with respect to the wild-type counterpart, as a consequence of an increased conformational rigidity probed by either NMR relaxation and MD simulations. Although dynamics in solution is other than fibrillar competence, the assessed properties of the mutant protein can be related with its reduced ability of forming fibrils when seeded in 20% trifluoroethanol. Copyright © 2013 John Wiley & Sons, Ltd.

The determination of differences in solvation free energies between related drug molecules remains an important challenge in computational drug optimization, when fast and accurate calculation of differences in binding free energy are required. In this study, we have evaluated the performance of five commonly used polarized continuum model (PCM) methodologies in the determination of solvation free energies for 53 typical alcohol and alkane small molecules. In addition, the performance of these PCM methods, of a thermodynamic integration (TI) protocol and of the Poisson-Boltzmann (PB) and generalized Born (GB) methods, were tested in the determination of solvation free energies changes for 28 common alkane-alcohol transformations, by the substitution of an hydrogen atom for a hydroxyl substituent. The results show that the solvation model D (SMD) performs better among the PCM-based approaches in estimating solvation free energies for alcohol molecules, and solvation free energy changes for alkane-alcohol transformations, with an average error below 1 kcal/mol for both quantities. However, for the determination of solvation free energy changes on alkane-alcohol transformation, PB and TI yielded better results. TI was particularly accurate in the treatment of hydroxyl groups additions to aromatic rings (0.53 kcal/mol), a common transformation when optimizing drug-binding in computer-aided drug design.

In this note we review recent progress about fluctuation relations and their applicability to free energy recovery in single molecule experiments. We underline the importance of the operational definition for the mechanical work and the non-invariance of fluctuation relations under Galilean transformations, both aspects currently amenable to experimental test. Finally we describe a generalization of the Crooks fluctuation relation useful to recover free energies of partially equilibrated states and thermodynamic branches., 15 pages, 8 figures

The performances of Bennett's acceptance ratio method and thermodynamic integration (TI) for the calculation of free energy differences in protein simulations are compared. For the latter, the standard trapezoidal rule, Simpson's rule, and Clenshaw-Curtis integration are used as numerical integration methods. We evaluate the influence of the number and definition of intermediate states on the precision, accuracy, and efficiency of the free energy calculations. Our results show that non-equidistantly spaced intermediate states are in some cases beneficial for the TI methods. Using several combinations of softness parameters and the λ power dependence, it is shown that these benefits are strongly dependent on the shape of the integrand. Although TI is more user-friendly due to its simplicity, it was found that Bennett's acceptance ratio method is the more efficient method. It is also the least dependent on the choice of the intermediate states, making it more robust than TI.

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.

Silencing in RNAi is strongly affected by guide-strand/target-mRNA mismatches. Target nucleation is thought to occur at positions 2–8 of the guide (“seed region”); successful hybridization in this region is the primary determinant of target binding affinity and hence target cleavage. To define a molecular basis for the target sequence selectivity in RNAi, we studied all possible distinct single mismatches in seven positions of the seed region — a total of 21 substitutions. We report results from soft-core thermodynamic integration simulations to determine changes in target binding free energies to Argonaute due to single mismatches in the guide strand, which arise during binding of an imperfectly matched target mRNA. In agreement with experiment, most mismatches impair target binding, consistent with a prominent role for binding affinity changes in RNAi sequence selectivity. Individual Argonaute residues located near the mismatched base pair are found to contribute significantly to binding affinity changes. We also employ this methodology to analyze the mismatch-dependent free energy changes for dissociation of a DNA·RNA hybrid from Argonaute, as a model for the escape of miRNAs from the silencing pathway. Several mismatched sequences of the miRNA have increased affinity to Argonaute, implying that some mismatches may reduce the probability for escape. Furthermore, calculations of base-substitution-dependent free energy changes for binding ssDNA reveal mild sequence sensitivity as expected for guide strand binding to Argonaute. Our findings give a thermodynamic basis for RNAi target sequence selectivity and suggest that miRNA mismatches may increase silencing effectiveness and thus could be evolutionarily advantageous.

The performance of enveloping distribution sampling (EDS) simulations to estimate free enthalpy differences associated with seven alchemical transformations of A-T into G-C base pairs at the netropsin binding site in the minor groove of a 13-base pair DNA duplex in aqueous solution is evaluated. It is demonstrated that sufficient sampling can be achieved with a two-state EDS Hamiltonian even for large perturbations such as the simultaneous transformation of up to three A-T into three G-C base pairs. The two parameters required to define the EDS reference state Hamiltonian are obtained automatically using a modified version of a scheme presented in earlier work. The sensitivity of the configurational sampling to a variation of these parameters is investigated in detail. Although for relatively small perturbations, that is, one base pair, the free enthalpy estimate depends only weakly on the EDS parameters, the sensitivity is stronger for the largest perturbation. Yet, EDS offers various convenient measures to evaluate the degree of sampling and thus the reliability of the free enthalpy estimate and appears to be an efficient alternative to the conventional thermodynamic integration methodology to obtain free energy differences for molecular systems.

The conformations and relative stabilities of folded and extended 3-fluoro-γ-aminobutyric acid (3F-GABA) conformers were studied using explicit solvation models. Geometry optimisations in the gas phase with one or two explicit water molecules favour folded and neutral structures containing intramolecular NH···O-C hydrogen bonds. With three or five explicit water molecules zwitterionic minima are obtained, with folded structures being preferred over extended conformers. The stability of folded versus extended zwitterionic conformers increases on going from a PCM continuum solvation model to the microsolvated complexes, though extended structures become less disfavoured with the inclusion of more water molecules. Full explicit solvation was studied with a hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme and molecular dynamics simulations, including more than 6000 TIP3P water molecules. According to free energies obtained from thermodynamic integration at the PM3/MM level and corrected for B3LYP/MM total energies, the fully extended conformer is more stable than folded ones by about -4.5 kJ mol(-1). B3LYP-computed (3)J(F,H) NMR spin-spin coupling constants, averaged over PM3/MM-MD trajectories, agree best with experiment for this fully extended form, in accordance with the original NMR analysis. The seeming discrepancy between static PCM calculations and experiment noted previously is now resolved. That the inexpensive semiempirical PM3 method performs so well for this archetypical zwitterion is encouraging for further QM/MM studies of biomolecular systems.

In the past analyses of the so-called van der Waals end point problem focused on thermodynamic integration. Here we investigate which of the recommendations, such as the need for soft-core potentials, are still valid when Bennett's acceptance ratio method is used. We show that in combination with Bennett's acceptance ratio method intermediate states characterized by the coupling parameter λ can be replaced by intermediate states in which Lennard-Jones interactions are turned on or off on an “atom by atom” basis. By doing so, there is no necessity to use soft-core potentials. In fact, one can compute free energy differences without dedicated code, making it possible to use any molecular dynamics program to compute alchemical free energy differences. Such an approach, which we illustrate by several examples, makes it possible to exploit the tremendous computational power of the graphics processing unit. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011

The one-step free energy perturbation approach can be applied to obtain conformational state-specific free energy differences (FEDs) associated with changes in force field parameters, and thus offers the possibility to consider conformational equilibria during force field parameterization. In this work, using the alanine decapeptide in explicit water solution as a model, the α-helical and β-hairpin state-specific FEDs associated with force field changes between two widely used parameter sets of the GROMOS force field, namely, 43A1 and 53A6, were determined using one-step perturbation. The results mostly deviated by only 1 kJ mol−1 in absolute or a few percent in relative values from thermodynamic integration results, suggesting that the convergence ranges of one-step perturbation were large enough to cover the substantial changes in nonbonded parameters between the two parameter sets. It was also found that one-step perturbation may give larger errors when the changes from the reference state include a large decrease in van der Waals radius, as indicated by the result for the β-hairpin state-specific free energy change going from 53A6 to 43A1. According to the free energy results, the α-helical state of the alanine decapeptide is destabilized by 15 kJ mol−1, i.e., 1.5 kJ mol−1 per residue, relative to the β-hairpin state when going from 43A1 to 53A6, in agreement with previous direct simulations in which native α-helices were often found to be unstable in simulations using 53A6, despite that the 53A6 parameters better reproduce a range of thermodynamic properties of small molecular systems. By applying one-step perturbation to analyze the effects of perturbing individual parameters, the differential stabilization of the two secondary structure states can be traced to the changes in van der Waals parameters, especially a van der Waals parameter involved in third-neighbor interactions. This study provides an example of the efficiency of one-step perturbation in force field development, reducing the computational cost by orders of magnitude. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the “double-wide” variant of the exponential formula. Use of the exponential formula in only strict “forward” or “backward” mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011

We attempt to optimize the efficiency of thermodynamic integration, as defined by the minimal number of unphysical intermediate states required for the computation of accurate and precise free energy differences. The suitability of various numerical quadrature methods is tested. In particular, we compare the trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Kronrod-Patterson, and Clenshaw-Curtis integration, as well as integration based on a cubic spline approximation of the integrand. We find that Simpson's rule and spline integration are already significantly more efficient that the trapezoidal rule, i.e., correct free energy differences can be obtained using fewer λ-states. We demonstrate that Simpson's rule can be used advantageously with nonequidistant values of the abscissa, which increases the flexibility of the method. Efficiency is enhanced even further if higher order methods, such as Gauss-Legendre, Gauss-Kronrod-Patterson, or Clenshaw-Curtis integration, are used; no more than seven λ-states, which in the case of Clenshaw-Curtis integration include the physical end states, were required for accurate results in all test problems studied. Thus, the performance of thermodynamic integration can equal that of Bennett's acceptance ratio method. We also show, however, that the high efficiency found here relies on the particular functional form of the soft-core potential used; overall, thermodynamic integration is more susceptible to the details of the hybrid Hamiltonian used than Bennett's acceptance ratio method. Therefore, we recommend Bennett's acceptance ratio method as the most robust method to compute alchemical free energy differences; nevertheless, scenarios when thermodynamic integration may be preferable are discussed.

Free-energy differences govern the equilibrium between bound and unbound states of a host and its guest molecules. The understanding of the underlying entropic and enthalpic contributions, and their complex interplay are crucial for the design of new drugs and inhibitors. In this study, molecular dynamics (MD) simulations were performed with inclusion complexes of α-cyclodextrin (αCD) and three monosubstituted benzene derivatives to investigate host–guest binding. αCD Complexes are an ideal model system, which is experimentally and computationally well-known. Thermodynamic integration (TI) simulations were carried out under various conditions for the free ligands in solution and bound to αCD. The two possible orientations of the ligand inside the cavity were investigated. Agreement with experimental data was only found for the more stable orientation, where the substituent resides inside the cavity. The better stability of this conformation results from stronger Van der Waals interactions and a favorable antiparallel host–guest dipole–dipole alignment. To estimate the entropic contributions, simulations were performed at three different temperatures (250, 300, and 350 K) and using positional restraints for the host. The system was found to be insensitive to both factors, due to the large and symmetric cavity of αCD, and the nondirectional nature of the host–guest interactions.

This report details an approach to improve the accuracy of free energy difference estimates using thermodynamic integration data (slope of the free energy with respect to the switching variable λ) and its application to calculating solvation free energy. The central idea is to utilize polynomial fitting schemes to approximate the thermodynamic integration data to improve the accuracy of the free energy difference estimates. Previously, we introduced the use of polynomial regression technique to fit thermodynamic integration data (Shyu and Ytreberg, J Comput Chem, 2009, 30, 2297). In this report we introduce polynomial and spline interpolation techniques. Two systems with analytically solvable relative free energies are used to test the accuracy of the interpolation approach. We also use both interpolation and regression methods to determine a small molecule solvation free energy. Our simulations show that, using such polynomial techniques and nonequidistant λ values, the solvation free energy can be estimated with high accuracy without using soft-core scaling and separate simulations for Lennard-Jones and partial charges. The results from our study suggest that these polynomial techniques, especially with use of nonequidistant λ values, improve the accuracy for ΔF estimates without demanding additional simulations. We also provide general guidelines for use of polynomial fitting to estimate free energy. To allow researchers to immediately utilize these methods, free software and documentation is provided via http://www.phys.uidaho.edu/ytreberg/software. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010

Molecular dynamics simulations were carried out to calculate free energy differences between the folded and unfolded states of wild type and mutant collagen model peptides. The calculated stability of the collagen models was compared with the severity of osteogenesis imperfecta. Free energy differences of Gly Xaa (Xaa: Ser, Cys, Glu, and Asp) mutations between the wild type and the mutants at position 15 of the model peptide were 3.8, 4.2, 5.6, and 8.8 kcal/mol, respectively. The corresponding free energy differences of a second Gly mutation at the same position in different chains were, on average, 1.3, 1.5, 2.9, and 5.4 kcal/mol, respectively. Free energy simulations were also performed to estimate the relative stability between an oxidized form and a reduced form of the mutants containing two Cys residues, which indicated that the mutant of the collagen-like peptide containing an intramolecular disulfide bond was more stable than the mutant containing one Cys residue but less stable than the wild type. The calculated free energy differences between an oxidized and a reduced form of the mutants containing two Cys residues are 0.8 and 2.6 kcal/mol for the disulfide bonds between Chains A and B and between Chains A and C, respectively. © 2010 Wiley Periodicals, Inc. Biopolymers 95: 182–193, 2011.

The accurate calculation of the free energy of interaction of protein-water-ligand systems has an important role in molecular recognition and drug design that is often not fully considered. We report free energy thermodynamic integration calculations used to evaluate the effects of inclusion, neglect, and targeting and removal (i.e., systematic substitution by ligand functional groups) of an important, tightly bound, water molecule in the SH3 domain of Abl tyrosine kinase. The effects of this water molecule on the free energies of interaction of several Abl-SH3 domain-ligand systems reveal that there is an unfavourable free energy change associated with its removal into the bulk solvent. Only three substitutions by an additional functional group (out of methyl, ethyl, hydroxyl, amino, and amide groups) in the phenyl ring of a tyrosine in the peptide ligand resulted in a favourable change in the free energy of binding upon replacement of the ordered water molecule. This computational approach provides a direct route to the systematic and rigorous prediction of the thermodynamic influence of ordered, structural water molecules on ligand modification and optimization in drug design by calculating free energy changes in protein-water-ligand systems.

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.

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, λ-dynamics treats the conventional "λ" as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present paper, we provide an overview of the theory of λ-dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how λ-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 λ-dynamics may be extended to facilitate modeling efforts in structure-based drug design.

The molecular-dynamics-based calculation of accurate free energy differences for biomolecular systems is a challenging task. Accordingly, convergence and accuracy of established equilibrium methods has been subject of many studies, often focusing at small test systems. In contrast, the potential of more recently proposed nonequilibrium methods, derived from the Jarzynski and Crooks equalities, has not yet fully been explored. Here, we compare the performance of these methods by calculating free energy differences for test systems at different levels of complexity and varying extent of the involved perturbations. We consider the interconversion of ethane into methanol, the switching of a tryptophane-sidechain in a tripeptide, and the binding of two different ligands to the globular protein snurportin 1. On the basis of our results, we suggest and assess a new nonequilibrium free energy method, Crooks Gaussian Intersection (CGI), which combines the advantages of existing methods. CGI is highly parallelizable and, for the test systems considered here, is shown to outperform the other studied equilibrium and nonequilibrium methods. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009

When estimating binding affinities of a ligand, which can exists in multiple forms, for a target molecule, one must consider all possible competing equilibria. Here, a method is presented that estimates the contribution of the protonation equilibria of a ligand in solution to the measured or calculated binding affinity. The method yields a correction to binding constants that are based on the total concentration of inhibitor (the sum of all ionized forms of the inhibitor in solution) to account for the complexed form of the inhibitor only. The method is applied to the calculation of the difference in binding affinity of two inhibitors, 2-phosphoglycolate (PGA) and its phoshonate analog 3-phosphonopropionate (3PP), for the glycolytic enzyme triosephosphate isomerase. Both inhibitors have three titrating sites and exist in solution as a mixture of different forms. In this case the form that actually binds to the enzyme is present at relative low concentrations. The contributions of the alternative forms to the difference in binding energies is estimated by means of molecular dynamics simulations and corrections. The inhibitors undergo a pKa shift upon binding that is estimated by ab initio calculations. An interesting finding is that the affinity difference of the two inhibitors is not due to different interactions in the active site of the enzyme, but rather due to the difference in the solvation properties of the inhibitors. Protein 2009. © 2008 Wiley-Liss, Inc.

All proteins, nucleic acids, and other biomolecules contain residues capable of exchanging protons with their environment. These proton transfer phenomena lead to pH sensitivity of many molecular processes underlying biological phenomena. In the course of biological evolution, Nature has invented some mechanisms to use pH gradients to regulate biomolecular processes inside cells or in interstitial fluids. Therefore, an ability to model protonation equilibria in molecular systems accurately would be of enormous value for our understanding of biological processes and for possible rational influence on them, like in developing pH dependent drugs to treat particular diseases. This work presents a derivation, by thermodynamic and statistical mechanical methods, of an expression for the free energy of a complex molecular system at arbitrary ionization state of its titratable residues. This constitutes one of the elements of modeling protonation equilibria. Starting from a consideration of a simple acid–base equilibrium of a model compound with a single tritratable group, we arrive at an expression which is of general validity for complex systems. The only approximation used in this derivation is the postulating that the interaction energy between any pair of titratable sites does not depend on the protonation states of all the remaining ionizable groups. © 2007 Wiley Periodicals, Inc. Biopolymers 89: 262–269, 2008. This article was originally published online as an accepted preprint. The “Published Online” date corresponds to the preprint version. You can request a copy of the preprint by emailing the Biopolymers editorial office at biopolymers@wiley.com

The structures and relative energies of the most stable conformers of both naked and microsolvated phenylalanine, Phe.(H(2)O)(n)(n=0-3), are calculated by density functional theory. For selected structures, coordination-constrained ab initio molecular dynamics simulations determine the proton-transfer mechanism connecting neutral and zwitterionic forms of Phe. The associated free-energy profiles are calculated by thermodynamic integration. While no zwitterionic free-energy minimum is found for naked Phe, microsolvation is found to stabilize the zwitterionic form. For cluster sizes n > or = 3, the proton-transfer equilibrium shifts towards the zwitterionic structure for specific proton-transfer pathways. The energetically most favourable interconversion path between the neutral and zwitterionic forms is through a H(2)O bridge with free-energy barriers as low as 14.4 kJ mol(-1) for Phe.(H(2)O)(3). The free energy required for breaking a carboxylic OH bond involved in intramolecular hydrogen bonding is typically lower than in the water-assisted case. However, the resulting zwitterion turns out to be unstable with respect to the backward proton-transfer reaction.

The theoretical study has been performed to refine the procedure for calculations of Gibbs free energy with a relative accuracy of less than 1 kcal/mol. Three benchmark intermolecular complexes are examined via several quantum-chemical methods, including the second-order Moller-Plesset perturbation (MP2), coupled cluster (CCSD(T)), and density functional (BLYP, B3LYP) theories augmented by Dunnings correlation-consistent basis sets. The effects of electron correlation, basis set size, and anharmonicity are systematically analyzed, and the results are compared with available experimental data. The results of the calculations suggest that experimental accuracy can be reached only by extrapolation of MP2 and CCSD(T) total energies to the complete basis set. The contribution of anharmonicity to the zero point energy and TDeltaSint values is fairly small. The new, economic way to reach chemical accuracy in the calculations of the thermodynamic parameters of intermolecular interactions is proposed. In addition, interaction energy (De) and free energy change (DeltaA) for considered species have been evaluated by Carr-Parrinello molecular dynamics (CPMD) simulations and static BLYP-plane wave calculations. The free energy change along the reaction paths were determined by the thermodynamic integration/"Blue Moon Ensemble" technique. Comparison between obtained values, and available experimental and conventional ab initio results has been made. We found that the accuracy of CPMD simulations is affected by several factors, including statistical uncertainty and convergence of constrained forces (TD integration), and the nature of DFT (density functional theory) functional. The results show that CPMD technique is capable of reproducing interaction and free energy with an accuracy of 1 kcal/mol and 2-3 kcal/mol respectively.

The relative stabilities of homochiral and heterochiral forms of selected dipeptides, AA, AS, AC, AV, AF, AD, AK, tripeptides, AAA, AVA, and an acetylpentapeptide, AcGLSFA, have been calculated using thermodynamic integration protocols and the GROMOS 53A6 force field. Integration pathways have been designed that produce minimal disturbance to the system, including the use of soft atoms, low-energy intermediates, and chiral inversion of the smaller amino acid in the peptide. Comparison of the results obtained by thermodynamic integration between the diastereomeric forms (in explicit water, at 300 K) and from exhaustive global minimum-energy searches for the individual dipeptides (implicit water, ε = 78, 0 K) suggests that entropic contributions to the relative stability of the chiral forms are important. This conclusion is supported by the results of explicit calculation of the effect of temperature on the relative stability of alanylvalylalanine diastereomers. The Gibbs free energy calculations predict that at ambient temperature and pressure homochiral dipeptides with small side chains or polar groups in the vicinity of the peptide backbone, AA, AS, and AD, are more stable than their heterochiral counterparts by fractions of a kJ/mol. For bigger side chains, AC, AV, AF, and AK, the heterochiral diastereomers appear to be more stable. Predicted relative stabilities are in line with observations reported in the literature for AE and YY. Excellent agreement is found for the calculated and experimentally determined relative stabilities of the diastereomers of the dipeptide AA and of all-L, AcGLSFA and its diastereomer containing D-serine in the central position. Addition of counterions to the solvent box has no significant effects on charged and neutral forms. From the present findings it would appear unlikely that the intrinsic stability difference between homo- and heterochiral dipeptides has been a driving force in a primordial selection process leading to the incorporation of amino acids with a single enantiomeric configuration in natural proteins. © 2006 Wiley Periodicals, Inc.

Density functional theory together with Car-Parrinello ab initio molecular dynamics simulation has been used to investigate the free energy profiles (FEP) of monomer capture in Grubbs- and SHOP-type olefin polymerization catalysts. The FEPs along the reaction coordinates at 300 K were determined directly by a point wise thermodynamic integration technique. Comparison between potential energy profile (PEP) and the FEP has been made. The results show that, for both catalysts, the PEP for the monomer ethylene uptake by the metal center is a typical Morse curve without energy barrier. However, a small barrier (1.8 kcal/mol for Grubbs catalyst and 2.4 kcal/mol for SHOP catalyst) exists on the FEP. The pi complexation energy on the FES at 300 K is higher by 10-12 kcal/mol over that on the PES. The differences between FES and PES are due to entropy contribution. Slow growth simulations on the ethylene capture process show that the ethylene attacks the metal center by an asynchronous mode. This indicates that the forming of the pi-bonding between the metal and ethylene is initiated by electrophilic attack of the metal to one of the ethylene carbons.

We show that equilibrium monomer and centre-of-mass concentration profiles of lattice polymers in a gravitational field, computed by Monte-Carlo simulations, provide an accurate and efficient road to the osmotic equation-of-state of polymer solutions, via a straightforward application of the hydrostatic equilibrium condition. The method yields the full equation of state over a wide range of concentrations from a single simulation, and does not suffer from significant finite size effects. It has been applied to self-avoiding walk polymer chains with nearest neighbour monomer attractions, from the good solvent to the theta solvent regimes. The consistency of the method has been carefully checked by varying the strength of the gravitational field., Submitted to ChemPhysChem. Special issue in honour of Prof. Michele Parrinello

We present results from molecular dynamics simulations of water near structured hydrophobic surfaces. The surface structures reported herein are a planar alkane crystal as a reference and crystals with a hole and a protrusion of approximately 2.5 nm diameter and 0.5 nm depth or height. All indicators show that surface structuring increases the hydrophobicity: The water density is reduced near the structure elements, and the number of residual contacts between water and the surface decreases by about 40 % with respect to the planar surface. Thermodynamic integration shows that the interfacial energy of the structured surfaces is about 7 mJ m(-2) higher for structured surfaces than for the planar surface. The hydrophobicity increases by a similar amount for the hole and the protrusion geometries compared to the planar surface.

The incorporation of explicit ions to mimic the effect of ionic strength or to neutralize the overall charge on a system in free energy calculations using molecular dynamics simulations is investigated. The difference in the free energy of hydration between two triosephosphate isomerase inhibitors calculated at five different ion concentrations is used as an example. We show that the free energy difference can be highly sensitive to the presence of explicit ions even in cases where the mutation itself does not involve a change in the overall charge. The effect is most significant if the molecule carries a net charge close to the site mutated. Furthermore, it is shown that the introduction of a small number of ions can lead to very severe sampling problems suggesting that in practical calculations convergence can best be achieved by incorporating either no counterions or by simulating at high ionic strength to ensure sufficient sampling of the ion distribution. (C) 2004 Wiley Periodicals, Inc.

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.

Thrombin inhibition is an important strategy for the treatment of thrombotic and embolic disorders. Despite high in vitro affinity for the enzyme, synthetic thrombin inhibitors have limited oral bioavailability due to low intestinal permeability. In a previous work, we obtained the following order for thrombin binding by using molecular dynamics (MD) simulations and the finite difference thermodynamic integration (FDTI) method: p-(2-oxo-1-propyl)benzamidine (POPBz) > p-ethylbenzamidine (PEBz) > p-(1-propyl)benzamidine (PPBz) > p-methylbenzamidine (PMBz) > benzamidine (Bz). As the octanol/water partition coefficient (log Po/w) is a thermodynamic property that may be related to drug transfer across biological membranes, we now turn our attention to the calculation by the FDTI method of relative log Po/w for the benzamidine derivatives. To examine the reliability of the method, log Po/w values for Bz and PMBz were determined experimentally. Experimental log Po/w of −1.02 and −0.89 were obtained for Bz and PMBz, and theoretical log Po/w of −0.12, −0.01, and −1.00 were obtained for PEBz, PPBz, and POPBz. Thus, we improved the hydrophobicity of Bz by adding the methyl, ethyl, propyl, and 2-oxo-propyl groups, but further substitutions should be conducted to raise log Po/w above zero and into a region where reasonable oral absorption might be expected. © 2005 Wiley Periodicals, Inc. Int J Quantum Chem, 2005

The incorporation of explicit ions to mimic the effect of ionic strength or to neutralize the overall charge on a system in free energy calculations using molecular dynamics simulations is investigated. The difference in the free energy of hydration between two triosephosphate isomerase inhibitors calculated at five different ion concentrations is used as an example. We show that the free energy difference can be highly sensitive to the presence of explicit ions even in cases where the mutation itself does not involve a change in the overall charge. The effect is most significant if the molecule carries a net charge close to the site mutated. Furthermore, it is shown that the introduction of a small number of ions can lead to very severe sampling problems suggesting that in practical calculations convergence can best be achieved by incorporating either no counterions or by simulating at high ionic strength to ensure sufficient sampling of the ion distribution. (C) 2004 Wiley Periodicals, Inc.

We calculated the free energy of solvation of the neutral analogs of 18 amino acid side-chains (not including glycine and proline) using the OPLS all-atom force field in TIP4P water, SPC water, and cyclohexane by molecular dynamics simulation and thermodynamic integration. The average unsigned errors in the free energies of solvation in TIP4P, SPC, and cyclohexane are 4.4, 4.9, and 2.1 kJ/mol respectively. Most of the calculated hydration free energies are not favorable enough compared to experiment. The largest errors are found for tryptophan, histidine, glutamic acid, and glutamine. The average unsigned errors in the free energy of transfer from TIP4P to cyclohexane and from SPC to cyclohexane are 4.0 and 4.1 kJ/mol, respectively. The largest errors, of more than 7.5 kJ/mol, are found for histidine, glutamine, and glutamatic acid.

The ability of the GROMOS96 force field to reproduce partition constants between water and two less polar solvents (cyclohexane and chloroform) for analogs of 18 of the 20 naturally occurring amino acids has been investigated. The estimations of the solvation free energies in water, in cyclohexane solution, and chloroform solution are based on thermodynamic integration free energy calculations using molecular dynamics simulations. The calculations show that while the force field reproduces the experimental solvation free energies of nonpolar analogs with reasonable accuracy the solvation free energies of polar analogs in water are systematically overestimated (too positive). The dependence of the calculated free energies on the atomic partial charges was also studied. (C) 2002 Wiley Periodicals, Inc.