Free Energy Perturbation Approach Research Articles

Multiscale Molecular Dynamics of Protein Aggregation

The 60's gave birth to the practical implementation of classical mechanics to unravel the dynamics and energetics of biomolecules. In the 70's the use of generalized force fields and more advanced integrative solutions to the microscopic understanding of nature (like hybrid QM/MM) were introduced. During the 80's, algorithms to obtain free energy values were further developed and in the 90's practical integration schemes of molecular mechanics force fields with other levels of detail (QM on one extreme and advances in implicit solvation on the other) were implemented in widely spread software. In the first decade of the XXIst century a considerable effort has been put in two seemingly discordant models for the simulation of biomolecules. On the one hand, extraordinary advances in computing technologies (both in terms of processor power and of new efficient parallel and distributed computing schemas) have allowed researchers to deal with bigger systems and longer simulations, reaching molecular processes including millions of particles or lying in the microsecond scale. On the other hand, the realization that the relevant answers to many biomolecular problems are not homogeneously distributed through the molecular structure, something already envisioned by the QM/MM pioneers more than three decades ago, has led researchers to find smart ways of putting different emphases on different ranges of the spatial or system time scale. In this context, e.g., molecular aggregation represents a paradigm for multiscalability, as molecular recognition can be understood with simple (semi-)macroscopic terms when the two fragments are far apart, while the atomic interactions need to be considered in full detail upon close distances. In this manuscript the current status of the techniques that use multiple scale representations of biomolecules are reviewed, and the findings are synthesized in a modular schema that can be extensively used when studying aggregation processes. It is shown that a smart alternative to brute force and massive computation of uninteresting regions in the all atom potential energy surface is the consideration of a simplified reference potential, explored thoroughly in the relevant regions, combined with a free energy perturbation approach that transforms this simple representation to a full atom representation.

Escorted free energy simulations

We describe a strategy to improve the efficiency of free energy estimates by reducing dissipation in nonequilibrium Monte Carlo simulations. This strategy generalizes the targeted free energy perturbation approach [C. Jarzynski, Phys. Rev. E 65, 046122 (2002)] to nonequilibrium switching simulations, and involves generating artificial, "escorted" trajectories by coupling the evolution of the system to updates in external work parameter. Our central results are: (1) a generalized fluctuation theorem for the escorted trajectories, and (2) estimators for the free energy difference ΔF in terms of these trajectories. We illustrate the method and its effectiveness on model systems.

Free energy perturbation approach for the rational engineering of the antibody for human hepatitis B virus

HzKR127 is the humanized monoclonal antibody effective for the neutralization of human hepatitis B virus. By means of the free energy perturbation (FEP) calculations based on molecular dynamic (MD) simulations, we examine the mutation-induced variations in the energetic and structural features associated with the interactions between HzKR127 and its antigen. N58A, Y96A, D97A, and D97A/Y102A mutants of HzKR127 are taken in account in this study for which the experimental data for relative efficacies with respect to the wild-type antibody are available. The results of the present MD–FEP simulation studies show that in order to enhance the affinity for the antigen, the engineering of HzKR127 should be made in such a way as to promote the dynamic stability of the overall protein conformation and that of the translational motion of the antigen in the antibody–antigen complex. The relative binding free energies of the four mutant antibodies obtained from MD–FEP calculations compare pretty well with the experimental mutagenesis data with the associated squared correlation coefficient of 0.96. This indicates that MD–FEP calculations may serve as a useful computational tool for rational antibody engineering. Discussed in detail are the differences in the structural features of antibody–antigen interactions between the wild-type and the mutant antibodies that are responsible for the change in binding affinities for the antigen.

Prediction of Relative Stability between TACE/Gelastatin and TACE/Gelastatin Hydroxamate

A gelastatins (1), natural MMP inhibitors, and their hydroxamate analogues (2) in TACE enzyme evaluated for discovery of potent TACE inhibitors. We have employed molecular dynamics simulations to compute the relative free energy of hydration and binding to TACE for gelastatin (1) and its hydroxamate analogue (2). The relative free energy difference is directly described in this article using the free energy perturbation approach as a means to accurately predict the TACE inhibitor of gelastatin analogues. The results show that the good agreement between the experimental and theoretical relative free energies of binding, gelastatin hydroxamate (2) binds stronger to TACE by -3.37 kcal/mol. The desolvation energy costs significantly reduced binding affinity, hydroxamate group associated with high desolvation energy formed strong favorable interactions with TACE with more than compensated for the solvation costs and therefore led to an improvement in relative binding affinity.

Molecular simulations of structures and solvation free energies of passivated gold nanoparticles in supercritical CO2

The interfacial structures and solvation free energies of gold nanoparticles passivated by self-assembled monolayer (SAM) of fluorinated alkanethiols in supercritical CO(2) (scCO(2)) have been studied by using classical molecular dynamics simulation. A fragment-based free-energy perturbation approach was developed here, in which the solvation free energy of passivated metal nanoparticles was partitioned into the contributions from the inner metal core and the outer SAM. This is a first-time attempt to directly simulate the solvation free energy of nano-objects in supercritical fluids. The simulation result suggests that the nanoparticles can be thermodynamically soluble at lower scCO(2) density but insoluble at higher density. We have demonstrated that this density dependence of solvation free energy can be ascribed to the effect of the surface SAM in scCO(2). The presence of solvent molecules greatly affects the morphology of SAM on nanoparticle. It was observed that increasing the chain length in SAM makes nanoparticles more solvophilic at lower scCO(2) density or more solvophobic at higher density. This solvation thermodynamics behavior has been correlated with the specific solvation structure of scCO(2) around the passivated nanoparticles.

Rapid Prediction of Solvation Free Energy. 1. An Extensive Test of Linear Interaction Energy (LIE).

The present study provides a comprehensive systematic analysis on the applicability of the linear interaction energy (LIE) approximation to the prediction of gas-to-water transfer (hydration) free energy. The study is based on molecular dynamics simulations in explicit solvent for an extensive and diverse hydration data set comprising 564 neutral compounds with measured hydration free energies, including a "traditional" data set and the more challenging drug-like SAMPL1 data set. A highly correlative LIE model was achieved without empirical scaling of the solute-solvent interaction energy terms along with a cavity term calibrated to the experiment. This model was particularly accurate for the "traditional" data set and of acceptable accuracy for the SAMPL1 data set, with mean-unsigned-errors below 1 kcal/mol and slightly above 2 kcal/mol, respectively. We have analyzed the sensitivity of the LIE model to several parameters such as continuum correction terms applied outside the explicit water shell, the impact of various charging methods, the applicability of single-conformer representation of the solute, and the inclusion of internal energy terms. The parameters with the greatest sensitivity are the charging methods used, with AM1BCC-SP (without AM1 geometry optimization) charges favored over AM1BCC-OPT and RESP charges. The inclusion of the change in intramolecular van der Waals and electrostatic energies between the solution and gas phases can also lead to improved prediction accuracies. Functional group based error analysis identified several chemical classes as minor outliers with systematic errors. A direct comparison of the LIE and free energy perturbation (FEP) approaches using the same force field and charging method shows that the LIE approximation is at least as accurate as the FEP approach with a reduction of computing time by at least 1 order of magnitude.

On the direct calculation of the free energy of quantization for molecular systems in the condensed phase

Using the path integral formalism or the Feynman-Hibbs approach, various expressions for the free energy of quantization for a molecular system in the condensed phase can be derived. These lead to alternative methods to directly compute quantization free energies from molecular dynamics computer simulations, which were investigated with an eye to their practical use. For a test system of liquid neon, two methods are shown to be most efficient for a direct evaluation of the excess free energy of quantization. One of them makes use of path integral simulations in combination with a single-step free energy perturbation approach and was previously reported in the literature. The other method employs a Feynman-Hibbs effective Hamiltonian together with the thermodynamic integration formalism. However, both methods are found to give less accurate results for the excess free energy of quantization than the estimate obtained from explicit path integral calculations on the excess free energy of the neon liquid in the classical and quantum mechanical limit. Suggestions are made to make both methods more accurate.

Optimization of cutting schemes for the evaluation of molecular electrostatic potentials in proteins via Moving-Domain QM/MM

This work presents new developments of the moving-domain QM/MM (MoD-QM/MM) method for modeling protein electrostatic potentials. The underlying goal of the method is to map the electronic density of a specific protein configuration into a point-charge distribution. Important modifications of the general strategy of the MoD-QM/MM method involve new partitioning and fitting schemes and the incorporation of dynamic effects via a single-step free energy perturbation approach (FEP). Selection of moderately sized QM domains partitioned between C (alpha) and C (from C=O), with incorporation of delocalization of electrons over neighboring domains, results in a marked improvement of the calculated molecular electrostatic potential (MEP). More importantly, we show that the evaluation of the electrostatic potential can be carried out on a dynamic framework by evaluating the free energy difference between a non-polarized MEP and a polarized MEP. A simplified form of the potassium ion channel protein Gramicidin-A from Bacillus brevis is used as the model system for the calculation of MEP.

Free energies of CO2∕H2 capture by p-tert-butylcalix[4]arene: A molecular dynamics study

The interactions of CO(2)H(2) with p-tert-butylcalix[4]arene (TBC4) were studied using potential of mean force (PMF) and free energy perturbation approaches. To the best of our knowledge, the present work is one of the first to employ the constrained mean force approach to evaluate solute selectivity by the TBC4 molecule. The computed PMFs for the interaction of CO(2)H(2) with a single TBC4 molecule establish that the interaction of CO(2) with the open end of the cage structure is attractive while the interaction with H(2) is repulsive. Free energy perturbation calculations were performed for the same two guest molecules with a pair of facing TBC4 molecules used as a representative model as found in the TBC4 molecular solid. At low temperature, both CO(2)H(2) have favorable interactions with the TBC4 pair, with the CO(2) interaction being considerably greater. These results are in agreement with recent experimental data showing considerable CO(2) uptake by TBC4 at moderate pressures.

Computational Study of the Effects of Mutations A156T, D168V, and D168Q on the Binding of HCV Protease Inhibitors.

The effect of the resistance mutations A156T, D168V, and D168Q in HCV protease on the binding of SCH 6, SCH 503034, VX-950, BILN-2061, and compound 1 was evaluated using the free energy perturbation (FEP) approach. All the inhibitors are highly potent against the wild-type enzyme, but their activity was affected differently by the mutants. A156T reduced the activity of SCH 503034, BILN-2061, and VX950 drastically (200-1000-fold) but that of SCH 6 only moderately (27-fold). SCH 503034, SCH 6, and VX-950 were not affected by either mutation D168V or D168Q, but these mutations conferred a high level of resistance to BILN-2061. Comparison of BILN-2061 with its acyclic analogue compound 1 emphasized the importance of inhibitor flexibility in overcoming drug resistance arising from the D168Q mutation. The results from FEP calculations compared well with experimental binding potencies within an error of <1 kcal/mol. Structural analysis was carried out to relate the resistance profiles to the atomic changes in the mutants.

Imaging the Migration Pathways for O2, CO, NO, and Xe Inside Myoglobin

Myoglobin (Mb) is perhaps the most studied protein, experimentally and theoretically. Despite the wealth of known details regarding the gas migration processes inside Mb, there exists no fully conclusive picture of these pathways. We address this deficiency by presenting a complete map of all the gas migration pathways inside Mb for small gas ligands (O2, NO, CO, and Xe). To accomplish this, we introduce a computational approach for studying gas migration, which we call implicit ligand sampling. Rather than simulating actual gas migration events, we infer the location of gas migration pathways based on a free-energy perturbation approach applied to simulations of Mb’s dynamical fluctuations at equilibrium in the absence of ligand. The method provides complete three-dimensional maps of the potential of mean force of gas ligand placement anywhere inside a protein-solvent system. From such free-energy maps we identify each gas docking site, the pathways between these sites, to the heme and to the external solution. Our maps match previously known features of these pathways in Mb, but also point to the existence of additional exits from the protein matrix in regions that are not easily probed by experiment. We also compare the pathway maps of Mb for different gas ligands and for different animal species.

Using a Charging Coordinate in Studies of Ionization Induced Partial Unfolding

The ionization of groups in proteins may sometimes involve a partial unfolding and/or water penetration. Unfortunately the corresponding structural changes might not be captured by microscopic free energy perturbation (FEP) approaches due to activation barriers that are not surmounted in nanosecond FEP simulations. This problem is apparent, for example, from mutation experiments that introduced ionizable groups in protein interiors and from the difficulties to reproduce the corresponding pKa changes by microscopic approaches. Here we develop a new general approach for treating such challenging cases. Our approach drives the protein structural change by increasing the charge of the ionized group beyond its physical value and thus overcoming the barriers for the partial unfolding by a physically consistent process. The potential of our approach is illustrated by the evaluation of the pKa of the Val66Glu mutant of staphylococcal nuclease. In this case it is first demonstrated that standard FEP approaches give extremely disappointing results for this pKa. On the other hand, our "overcharging" approach gives a much more realistic result. We believe that the present approach represents a breakthrough in FEP studies of ionizable residues in proteins, and expect this strategy to be useful in studies of a wide range of challenging problems including simulations of hydrogen exchange processes.

Free Energy Simulations of Uncatalyzed DNA Replication Fidelity: Structure and Stability of T·G and dTTP·G Terminal DNA Mismatches Flanked by a Single Dangling Nucleotide

A reference system for DNA replication fidelity was studied by free energy perturbation (FEP) and linear interaction energy (LIE) methods. The studied system included a hydrated duplex DNA with the 5'-CG dangling end of the templating strand, and dCTP4-.Mg2+ or dTTP4-.Mg2+ inserted opposite the dangling G to form a correct (i.e., Watson-Crick) or incorrect (i.e., wobble) base pair, respectively. The average distance between the 3'-terminal oxygen of the primer strand and the alpha-phosphorus of dNTP was found to be 0.2 A shorter for the correct base pair than for the incorrect base pair. Binding of the incorrect dNTP was found to be disfavored by 0.4 kcal/mol relative to the correct dNTP. We estimated that improved binding and more near-attack configurations sampled by the correct base pair should translate in aqueous solution and in the absence of DNA polymerase into a six times faster rate for the incorporation of the correct dNTP into DNA. The accuracy of the calculated binding free energy difference was verified by examining the relative free energy for melting duplex DNA containing GC and GT terminal base pairs flanked by a 5' dangling C. The calculated LIE and FEP free energies of 1.7 and 1.1 kcal/mol, respectively, compared favorably with the experimental estimate of 1.4 kcal/mol obtained using the nearest neighbor parameters. To decompose the calculated free energies into additive electrostatic and van der Waals contributions and to provide a set of rigorous theoretical data for the parametrization of the LIE method, we suggested a variant of the FEP approach, for which we coined a binding-relevant free energy (BRFE) acronym. BRFE approach is characterized by its unique perturbation pathway and by its exclusion of the intramolecular energy of a rigid part of the ligand from the total potential energy.

Molecular-dynamics evaluation of fluid-phase equilibrium properties by a novel free-energy perturbation approach: Application to gas solubility and vapor pressure of liquid hexane

A novel free-energy perturbation method is developed for the computation of the free energy of transferring a molecule between fluid phases. The methodology consists in drawing a free-energy profile of the target molecule moving across a binary-phase structure built in the computer. The novelty of the method lies in the difference of the definition of the free-energy profile from the common definition. As an important element of the method, the process of making a correction to the transfer free energy with respect to the cutoff of intermolecular forces is elucidated. In order to examine the performance of the method in the application to fluid-phase equilibrium properties, molecular-dynamics computations are carried out for the evaluation of gas solubility and vapor pressure of liquid n-hexane at 298.15 K. The gas species treated are methane, ethane, propane, and n-butane, with the gas solubility expressed as Henry's constant. It is shown that the method works fine and calculated results are generally in good agreement with experiments. It is found that the cutoff correction is strikingly large, constituting a dominant part of the calculated transfer free energy at the cutoff of 8 A.

Through the channel and around the channel: Validating and comparing microscopic approaches for the evaluation of free energy profiles for ion penetration through ion channels.

Microscopic calculations of free energy profiles for ion transport through biological ion channels present a very serious challenge to modern simulation approaches. The main problem is due to the major convergence problems associated with the heterogeneous landscape of the electrostatic environment in ion channels and with the need to evaluate the profile associated with the transfer of the ion from bulk water to the channel environment. This problem is compounded by the lack of reliable and relevant benchmarks that can discriminate between alternative approaches. The present study is aimed at reducing the above problems by defining benchmarks that are directly relevant to ion channels and can also give converging results. This is done by constructing a series of models of a truncated gramicidin channel with different numbers of water molecules and by comparing the profiles for going around the channel and through the channel. These discriminating models are then used to validate and compare the adiabatic charging free energy perturbation (FEP) approach combined with an umbrella sampling approach (Warshel, A. J. Phys. Chem. 1982, 86, 2218) and the potential of mean force (PMF) approach used frequently in studies of ion channels. It is found that both approaches work quite well until one moves to the case of the fully solvated channel. In this limit, the PMF approach may give different results for the overall work of going through the channel and around the channel, while the FEP approach gives physically consistent results. The present benchmark also indicates that the weighted histogram analysis method (WHAM) approach does not offer a significant advantage over earlier approaches at least as much as studies of ion channels are concerned. Finally, it is concluded that the FEP approach may be more useful in evaluating the overall barrier for moving ions from water to ion channels and that in some cases it might be beneficial to use the FEP approach for selective points along the channel and then to connect these points by PMF calculations.

Free energy perturbation approach to the critical assessment of selective cyclooxygenase-2 inhibitors

The discovery of selective cyclooxygenase-2 (COX-2) inhibitors represents a major achievement of the efforts over the past few decades to develop therapeutic treatments for inflammation. To gain insights into designing new COX-2-selective inhibitors, we address the energetic and structural basis for the selective inhibition of COX isozymes by means of a combined computational protocol involving docking experiment, force field design for the heme prothetic group, and free energy perturbation (FEP) simulation. We consider both COX-2- and COX-1-selective inhibitors taking the V523I mutant of COX-2 to be a relevant structural model for COX-1 as confirmed by a variety of experimental and theoretical evidences. For all COX-2-selective inhibitors under consideration, we find that free energies of binding become less favorable as the receptor changes from COX-2 to COX-1, due to the weakening and/or loss of hydrogen bond and hydrophobic interactions that stabilize the inhibitors in the COX-2 active site. On the other hand, COX-1-selective oxicam inhibitors gain extra stabilization energy with the change of residue 523 from valine to isoleucine because of the formations of new hydrogen bonds in the enzyme-inhibitor complexes. The utility of the combined computational approach, as a valuable tool for in silico screening of COX-2-selective inhibitors, is further exemplified by identifying the physicochemical origins of the enantiospecific selective inhibition of COX-2 by alpha-substituted indomethacin ethanolamide inhibitors.

An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.

Structure–activity relationship for extracellular block of K+ channels by tetraalkylammonium ions

External tetraalkylammonium ion binding to potassium channels is studied using microscopic molecular modelling methods and the experimental structure of the KcsA channel. Relative binding free energies of the KcsA complexes with Me4N+, Et4N+, and n-Pr4N+ are calculated with the molecular dynamics free energy perturbation approach together with automated ligand docking. The four-fold symmetry of the entrance cavity formed by the Tyr82 residues is found to provide stronger binding for the D2d than for the S4 conformation of the ligands. In agreement with experiment the Et4N+ blocker shows several kcal/mol better binding than the other tetraalkylammonium ions.

How dihydrofolate reductase facilitates protonation of dihydrofolate.

Dihydrofolate Reductase (DHFR) catalyzes the reduction of dihydrofolate (H2F) to tetrahydrofolate. On the basis of 10-12.5 ns molecular dynamics simulations of two conformations (closed and occluded) of the ternary DHFR/NADPH/H2F complex from Escherichia coli and a free energy perturbation approach, we have calculated the pKa value for the N5 atom in H2F. Our results suggest that the N5 atom in H2F is responsible for the pH dependency of the catalyzed reaction, meaning that DHFR facilitates protonation of H2F by approximately 4 pKa units. The mechanism behind this increase is due to favorable electrostatic interactions between the Asp27 residue and a proton at the N5 atom. The electrostatic interactions are enhanced by a hydrophobic active site, which to a large extent is made hydrophobic by the M20 loop in DHFR. Moreover, we find that the conformation imposed on H2F by DHFR to some extent also favors protonation of the N5 atom. Our results add support to previous findings and suggestions by Callender and co-workers [e.g., Deng, J.; Callender, R. J. Am. Chem. Soc. 1998, 120, 7730-7737] and explain why mutation of Asp27 may lead to severely reduced activity at neutral pH.

Structure⧸Function Correlations of Proteins using MM, QM⧸MM, and Related Approaches: Methods, Concepts, Pitfalls, and Current Progress

Publisher Summary This chapter focuses on the progress and the current state of structure/function correlations of biomolecules. One of the primary goals of modern molecular biology is to provide a reliable correlation between the structure and function of biomolecules. The chapter reviews the classical force field approaches and their use in macromolecular modeling—such as molecular mechanics, molecular dynamics, and evaluation of free energy using the free energy perturbation (FEP) and linear response approximation (LRA) approaches. It also emphasizes the need for proper long range treatments and the problem with using a single very long trajectory to explore biological functions. The chapter considers models for studies of electrostatic energies in proteins and their use in evaluation of pKa's, redox potentials, and binding free energies. The main part of the review deals with the modeling of enzyme catalysis. This part will consider several modeling approaches ranging from the empirical valence bond (EVB) to different molecular hybrid quantum mechanics/molecular mechanics (QM/MM) orbital approaches as well as approaches that attempt to treat the entire enzyme quantum mechanically. The potential and actual current performance of QM/MM and related approaches is presented. It is pointed out that while the EVB provides the most reliable results, alternative QM (MO)/MM methods are starting to provide reasonable results as well. There is no doubt about massive progress in the use of reliable QM/MM and related approaches. It is also certain that there will be a growing realization that consistent simulation studies provide the ultimate tool for understanding structure/function correlations in proteins.