Relative Free Energy Calculations Research Articles

2P006 自由エネルギー変分原理に基づく相対的結合自由エネルギー計算とそのFK506結合タンパク阻害剤への応用(蛋白質(構造・構造機能相関),口頭発表,第45回日本生物物理学会年会)

2P006 自由エネルギー変分原理に基づく相対的結合自由エネルギー計算とそのFK506結合タンパク阻害剤への応用(蛋白質(構造・構造機能相関),口頭発表,第45回日本生物物理学会年会)

Interaction Simulation of hERG K+ Channel with Its Specific BeKm-1 Peptide: Insights into the Selectivity of Molecular Recognition

Potassium channels show a huge variability in the affinity when recognizing enormous bioactive peptides, and the elucidation of their recognition mechanism remains a great challenge due to an undetermined peptide-channel complex structure. Here, we employed combined computation methods to study the specific binding of BeKm-1 peptide to the hERG potassium channel, which is an essential determinant of the long-QT syndrome. By the use of a segment-assembly homology modeling method, the closed-state hERG structure containing unusual longer S5P linker was successfully constructed. It has a "petunia" shape, while four "petals" of symmetrically distributed S5P segments always decentralize. Starting from the hERG and BeKm-1 structures, a considerably reasonable BeKm-1-hERG complex structure was then screened out and identified by protein-protein docking, molecular dynamics (MD) simulations, and calculation of relative binding free energies. The validity of this predicted complex was further assessed by computational alanine-scanning, with the results correlating reasonably well with experimental data. In the novel complex structure, four considerably flexible S5P linkers are far from the BeKm-1 peptide. The BeKm-1 mainly uses its helical region to associate the channel outer vestibule, except for the S5P linker region; however, structural analysis further implies this neutral pore region with wiggling S5P linker is highly beneficial to the binding of BeKm-1 with lower positive charges. The most critical Lys18 of BeKm-1 plugs its side chain into the channel selectivity filter, while the secondarily important Arg20 forms three hydrogen bonds with spatially neighboring residues in the hERG channel. Different from the classical peptide-K+ channel interaction mainly induced by electrostatic interaction, a synergetic effect of the electrostatic and van der Waals interactions was found to mediate the molecular recognition between BeKm-1 and the hERG channel. And this specific binding process is revealed to be a dynamic change of reduction of binding free energy and conformational rearrangement mainly in the interface of both BeKm-1 and the hERG channel. All these structural and energy features yield deep insights on the high selective binding mechanism of hERG-specific peptides, present a diversity of peptide-K+ channel interactions, and also provide important clues to further study structure-function relationships of the hERG channel.

Theoretical Study of the Enol Imine ↔ Enaminone Tautomeric Equilibrium in Organic Solvents

The tautomeric enol imine <--> enaminone (phenol <--> quinone) equilibrium of the 1-hydroxy-2-naphthaldehyde Schiff base (2-phenyliminomethyl-naphthalen-1-ol) was investigated by density functional theory (B3LYP) and ab initio (MP2) methods in the IEF-PCM polarizable continuum dielectric solvent approximation and by a combined ab initio + FEP/MC study by considering an explicit solvent model. Special emphasis was put on the effect of solvation on this equilibrium by using an apolar (CCl4), polar aprotic (CH3CN), and polar protic (CH3OH) solvent. Compared with experimental tautomerization Gibbs free energies, the IEF-PCM/B3LYP calculations apparently overestimate the stability of the quinone form both when the 6-31G(d,p) and the 6-311++G(d,p) basis sets are applied. IEF-PCM/MP2 studies with the above basis sets predict the preference of the aromatic phenol tautomer, in contrast to the experiment in methanol and acetonitrile solvent. Calculation of the total relative free energy as DeltaG(tot) = DeltaE(int)(IEF-PCM/QCISD(T)/6-31G(d)) + DeltaG(solv, FEP/MC) + DeltaG(thermal) provided agreement with the experimental values up to 0.6 kcal/mol in the three solvents, and the predominant tautomer was always correctly predicted. In-solution relevant atomic charges, derived by a fit to the molecular electrostatic potential generated by the IEF-PCM/B3LYP/6-31G(d,p) wave function, show strong dependence on the fitting procedure (CHELPG or RESP) and are fairly insensitive to the chemical nature of the actual solvent. Use of the CHELPG charges in FEP/MC simulations revealed to be superior in comparison with the use of the RESP charge set.

Grid-based steered thermodynamic integration accelerates the calculation of binding free energies

The calculation of binding free energies is important in many condensed matter problems. Although formally exact computational methods have the potential to complement, add to, and even compete with experimental approaches, they are difficult to use and extremely time consuming. We describe a Grid-based approach for the calculation of relative binding free energies, which we call Steered Thermodynamic Integration calculations using Molecular Dynamics (STIMD), and its application to Src homology 2 (SH2) protein cell signalling domains. We show that the time taken to compute free energy differences using thermodynamic integration can be significantly reduced: potentially from weeks or months to days of wall-clock time. To be able to perform such accelerated calculations requires the ability to both run concurrently and control in realtime several parallel simulations on a computational Grid. We describe how the RealityGrid computational steering system, in conjunction with a scalable classical MD code, can be used to dramatically reduce the time to achieve a result. This is necessary to improve the adoption of this technique and further allows more detailed investigations into the accuracy and precision of thermodynamic integration. Initial results for the Src SH2 system are presented and compared to a reported experimental value. Finally, we discuss the significance of our approach.

Free energies of ligand binding for structurally diverse compounds.

The one-step perturbation approach is an efficient means to calculate many relative free energies from a common reference compound. Combining lessons learned in previous studies, an application of the method is presented that allows for the calculation of relative binding free energies for structurally rather diverse compounds from only a few simulations. Based on the well known statistical-mechanical perturbation formula, the results do not require any empirical parameters, or training sets, only limited knowledge of the binding characteristics of the ligands suffices to design appropriate reference compounds. Depending on the choice of reference compound, relative free energies of binding rigid ligands to the ligand-binding domain of the estrogen receptor can be obtained that show good agreement with the experimental values. The approach presented here can easily be applied to many rigid ligands, and it should be relatively easy to extend the method to account for ligand flexibility. The free-energy calculations can be straightforwardly parallelized, allowing for an efficient means to understand and predict relative binding free energies.

Simulation of the Interaction Between ScyTx and Small Conductance Calcium-Activated Potassium Channel by Docking and MM-PBSA

Computational methods are employed to simulate interaction of scorpion toxin ScyTx in complex with the small conductance calcium-activated potassium channel rsk2. All of available 25 structures of ScyTx in the Protein Data Bank determined by NMR were considered for improving performance of rigid protein docking of ZDOCK. Four main binding modes were found among a large number of predicted complexes by using clustering analysis, screening with expert knowledge, energy minimization, and molecular dynamics simulations. The quality and validity of the resulting complexes were further evaluated by molecular dynamics simulations with the generalized Born solvation model and by calculation of relative binding free energies with the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) in the AMBER 7 suit of programs. The complex formed by the 22nd structure of the ScyTx and rsk2 channel was identified as the most favorable complex by using a combination of computational methods, which contain further introduction of flexibility without restraining residue side chain. From the resulted spatial structure of the ScyTx and rsk2 channel, ScyTx associates the mouth of the rsk2 channel with α-helix rather than β-sheet. Structural analysis first revealed that Arg 13 played a novel and vital role of blocking the pore of the rsk2 channel, whose role is remarkably different from that of highly homologous scorpion toxin P05. Between the interfaces in the ScyTx-rsk2 complex, strong electrostatic interaction and hydrogen bonds exist between Arg 13 of ScyTx and Gly-Tyr-Gly-Asp sequential residues located in the four symmetrical chains of the pore region. Simultaneously, five hydrogen bonds between Arg 6 of ScyTx and Asp 341(C), Val 366(C), and Pro 367(C), and electrostatic interaction between Arg 6 of ScyTx and Asp 364(B) and Asp 341(C) are also found by structural analysis. In addition, His 31 located at the C-terminal of ScyTx is surrounded by Val 342(A), Asp 364(A), Met 365(A), Pro 367(B), and Asn 366(B) within a contact distance of 4.0 Å. These simulation results are in good agreement with experimental data and can effectively explain the binding phenomena between ScyTx and the potassium channel at the level of molecular spatial structure. The consistency between results of molecular modeling and experimental data strongly suggests that our spatial structure model of the ScyTx-rsk2 complex is reasonable. Therefore, molecular docking combined with molecular dynamics simulations followed by molecular mechanics Poisson-Boltzmann surface area analysis is an attractive approach for modeling scorpion toxin-potassium channel complexes a priori for further biological studies.

Development of a Quantum Mechanics-Based Free-Energy Perturbation Method: Use in the Calculation of Relative Solvation Free Energies

Free-energy perturbation (FEP) is considered the most accurate computational method for calculating relative solvation and binding free-energy differences. Despite some success in applying FEP methods to both drug design and lead optimization, FEP calculations are rarely used in the pharmaceutical industry. One factor limiting the use of FEP is its low throughput, which is attributed in part to the dependence of conventional methods on the user's ability to develop accurate molecular mechanics (MM) force field parameters for individual drug candidates and the time required to complete the process. In an attempt to find an FEP method that could eventually be automated, we developed a method that uses quantum mechanics (QM) for treating the solute, MM for treating the solute surroundings, and the FEP method for computing free-energy differences. The thread technique was used in all transformations and proved to be essential for the successful completion of the calculations. Relative solvation free energies for 10 structurally diverse molecular pairs were calculated, and the results were in close agreement with both the calculated results generated by conventional FEP methods and the experimentally derived values. While considerably more CPU demanding than conventional FEP methods, this method (QM/MM-based FEP) alleviates the need for development of molecule-specific MM force field parameters and therefore may enable future automation of FEP-based calculations. Moreover, calculation accuracy should be improved over conventional methods, especially for calculations reliant on MM parameters derived in the absence of experimental data.

Tools for channels: moving towards molecular calculations of gating and permeation in ion channel biophysics

Recent X-ray structures of voltage gated potassium channels provide an exciting opportunity to connect molecular structures with measured biological function. Two of the most important connections for these channels are: first, to the molecular basis behind selectivity and the associated free energy profile underlying ionic current flow and, second, to a true molecular understanding of the large-scale conformational transitions that underlie voltage dependent gating. But, existing computational tools need to be further developed to reach these goals. In this contribution to the symposia on sampling methods we outline our dynamic importance sampling method for sampling large-scale conformational transitions as well as our studies with non-equilibrium work events and equilibrium overlap sampling (OS) methods for sampling events related to the calculation of relative free energies.

The Development of Replica-Exchange-Based Free-Energy Methods

The calculation of relative free energies that involve large reorganizations of the environment is one of the great challenges of condensed-phase simulation. Such calculations are of particular importance in protein−ligand free-energy calculations. To meet this challenge, we have developed new free-energy techniques that combine the advantages of the replica-exchange method with free-energy perturbation (FEP) and finite-difference thermodynamic integration (FDTI). These new techniques are tested and compared with FEP, FDTI, and the adaptive umbrella weighted histogram analysis method (AdUmWHAM) on the challenging calculation of the relative hydration free energy of methane and water. This calculation involves a large solvent configurational change. Through the use of replica-exchange moves along the λ-coordinate, the configurations sampled along λ are allowed to mix, which leads to dramatic improvements in solvent configurational sampling, an efficient reduction of random sampling error, and a reduction o...

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

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

Staging Is More Important than Perturbation Method for Computation of Enthalpy and Entropy Changes in Complex Systems

Entropy and enthalpy compensation are frequently observed in thermodynamic measurements. To understand the molecular details of these compensating changes, molecular dynamics or Monte Carlo methods may, in principle, be used. Yet, in practice, computational methods to evaluate changes in entropy and enthalpy usually have much greater (e.g. order of magnitude worse) error than relative free energy calculations. In this paper, we examine ways to improve the computational ability to determine changes in both enthalpy and entropy. Toward that end we consider five different perturbation schemes for enthalpy determination and also consider the relative importance of staging the computation (i.e., dividing the entire perturbation calculation into a series of consecutive smaller ones) along the reaction coordinate. Two model systems are used for evaluation: a system of Lennard-Jones spheres, where the computation is directed to the addition/deletion of a sphere, and the alchemical transformation of an anion in e...

Free Energy Calculations of Conformational Equilibrium of Chorismate in Water. The Role of Solute Polarization

The conformational equilibrium of chorismate molecule in water has been studied with Monte Carlo free energy perturbation simulations. Relative solvation free energies of chorismate conformations have been calculated using three types of perturbation paths depending on the intermolecular solute−solvent energy terms considered. First, smooth transformation between solute conformations are performed without consideration of chorismate site charges. Second, the solvation free energy contribution because of the inclusion of electrostatic solute−solvent interactions are calculated. Finally, the role of polarization forces are studied introducing progressively atom polarizabilities in the chorismate molecule. The calculation of relative free energies through these three perturbation paths allows the comparison of the different energetic forces, and improves the convergence of the results. It has been found that the inclusion of solute polarization is crucial to obtain diequatorial−diaxial conformational proportions compatible with the experimental values. Different procedures to calculate solute polarization are analyzed depending on the way the intramolecular interactions among the polarization sites of chorismate are treated. Two types of water models, polarizable and nonpolarizable, have been considered. A hydrogen bond analysis and a study of the role of the charge variations among chorismate conformations are also done.

One-Step Perturbation Methods for Solvation Free Energies of Polar Solutes

The one-step perturbation approach to free energy calculations permits the calculation of relative free energies for a number of end states from a single simulation of a carefully chosen reference state. This approach has successfully been used for the calculation of relative and absolute solvation free energies of nonpolar solutes in water, as well as relative binding free energies for a family of nonpolar protein ligands. We have investigated several possible approaches to extend one-step perturbation methods to also obtain accurate solvation free energies of polar species. A “soft dipole” reference state permits the determination of accurate relative solvation free energies for polar solutes (average unsigned error < 2.0 kJ/mol) but does not yield accurate absolute solvation free energies. However, contrary to expectations, simulations of a neutral reference state even yield accurate (average unsigned error < 2.9 kJ/mol) absolute solvation free energies for polar solutes when translational and rotational sampling of the solute is included. In general, the choice of an appropriate reference state for one-step perturbation is a challenge because of complex many-body effects. Nonetheless, the one-step perturbation method combined with an unphysical reference state yields accurate free energy estimates with an efficiency several times greater than traditional approaches.

Simulations of the estrogen receptor ligand-binding domain: affinity of natural ligands and xenoestrogens.

We have carried out molecular dynamics (MD) simulations and free energy calculations on the alpha-subtype of the human estrogen receptor ligand-binding domain (ERalpha LBD) complexed with a number of known agonists and putative xenoestrogens. Our dynamical simulations of ligand-receptor complexes underscore the highly structured nature of the complex and offer some interesting insights into the structure-activity relationship (SAR) for these ligands. With traditional thermodynamic integration (TI) calculations, we calculate relative binding free energies for three known agonists, in good agreement with experimental values. The sheer number of possible xenoestrogenic compounds makes an approach using traditional free energy calculations unfeasible. Instead, we have made use of a single-step perturbation methodology that allows the calculation of relative free energies for a large number of related polyaromatic hydrocarbons (PAHs) from a single simulation. Our results show good (maximum deviation 3.3 kJ mol(-1)) agreement with experimental data, suggesting the possibility of large-scale xenoestrogen screening in silico to obtain strongly estrogenic compounds for subsequent experimental testing.

Calculation of Relative Binding Free Energies of Peptidic Inhibitors to HIV-1 Protease and Its I84V Mutant

A methodology is presented for calculating relative binding free energies of enzyme−inhibitor associations in aqueous solvent. The methodology uses synthesis of semiempirical quantum chemistry to determine the protonation state of important residues in the enzyme active site, molecular mechanics to determine the gas-phase energetic contributions to the relative binding free energy, and dielectric continuum solvation to calculate electrostatic hydration contributions. The methodology is then applied to the calculation of the relative binding free energy of the inhibitors KNI-272, Ro31-8959, L-735,524, and A-77003 to HIV-1 protease and its I84V mutant. The calculated relative binding free energy is sensitive to the active-site protonation state of the aspartic acid residues of HIV-1 protease. The protonation state is inhibitor dependent. Given a particular protonation state, it was found that quantitatively accurate relative binding free energies could only be achieved when solvent effects were included. Three categories of binding were found. In the first, the change in binding free energy due to mutation is mainly due to the change in enthalpic interactions within the inhibitor−enzyme complex (Ro31-8959). In the second (L-735,524 and A-77003), the change in affinity is caused both by a change in enthalpic interactions within the enzyme and by a change in the hydration energy of the enzyme and inhibitor−enzyme complexes. In the third case (KNI-272), the change in affinity is mainly a solvent effectit is due to changes in hydration of the enzyme only. In all cases, it was found that the I84V mutant enzyme was more stable than the wild-type enzyme. This alone (without consideration of the inhibitor−enzyme complexes) can qualitatively explain the reduction in binding affinity due to mutation.

Advancing beyond the atom-centered model in additive and nonadditive molecular mechanics

A computational approach to the inclusion of off-center charges in both additive and nonadditive molecular mechanics calculations is presented. The additional sites in the molecular skeleton are placed in the approximate locations of the chemically intuitive electron lone pair, and are treated as formal particles throughout the calculation. The increase in the number of charge sites results in overall improvement in the energy associated with the angular dependence of hydrogen bonds and improved statistical accuracy of the electrostatic potential derived charges. The addition of lone pairs also results in improved accuracy in relative solvation free energy calculation for the pyridine to benzene and methanol to methane mutations. Because the number of atoms that require lone pairs is small, the extra accuracy can be achieved with little computational overhead. © 1997 John Wiley & Sons, Inc. J Comput Chem18: 1632–1646, 1997

Molecular modeling studies of 2-substituted cephalosporin esters binding to human leukocyte elastase

Summary A series of 17 putative human leukocyte elastase (HLE) inhibitors belonging to the family of cephalosporin esters have been studied by means of molecular modeling techniques. The optimized conformation (molecular mechanics) of HLE was used as input for a series of simulated annealing calculations meant to locate a lower energy minimum than that identified by the previous minimization. Manual and automated docking experiments with all the 17 compounds have then been made with elastase in the lowest energy conformation found by the annealing phase. Further molecular dynamics studies of the complexes made by the enzyme with two different inhibitors, as well as a relative free-energy calculation of the two above-mentioned complexes, have been performed in order to get information about the recognition/binding process. The relationships between the steric, electrostatic and lipophilic descriptors (some of which were obtained with semiempirical MO calculations) of the inhibitors, their intermolecular non-bond interaction energies (INIs) and their IC50s have been described with a series of statistical equations. PLS and MLR-like models explaining such relationships have been generated.

Calculation of Relative Binding Free Energies and Configurational Entropies: A Structural and Thermodynamic Analysis of the Nature of Non-polar Binding of Thrombin Inhibitors Based on Hirudin 55–65

Free energy calculations were carried out on a series of exosite-binding inhibitors of thrombin. These inhibitors are based on the C-terminal fragment of hirudin and have the sequence Phe-Glu-Glu-Ile H59-Pro-Glu-Glu-Tyr-Leu, where the superscript over Ile indicates its relative position in the natural sequence of hirudin. In this study, the effect of replacing Ile H59with ten other non-polar amino acids was examined. Three preferred interaction sites for methyl/methylene groups for the various Xaa H59side-chains in the complex were identified from conformational search calculations. The corresponding thermodynamic changes were determined using combination of systematic search and energy minimization in a manner that locates the local minima in the system and in the process simultaneously builds up the partition function. The free energy, internal energy and entropic contributions are readily calculated from the partition function. Very good agreement in the resulting relative binding free energies was obtained between theory and experiment. The calculations allowed us to dissect out the enthalpic, entropic and solvation contributions to ΔΔ G. The contribution from desolvation was found to be relatively weak. The binding of these non-polar side-chains to thrombin is found to be driven mainly by favourable protein-ligand interactions rather than by the desire for non-polar groups to be desolvated. We also find that the configurational entropy contributes about 0.48 kcal/mol (0.81 kT) in average for each torsional angle “frozen” in binding.

Calculation of relative binding of free energy of vancomycin with diacetyl- l-lysyl- d-alanyl- d-alanine and diacetyl- l-lysyl- l-alanyl- l-alanine

By means of a free energy perturbation method implemented with molecular dynamics, the relative binding of free energies of vancomycin with Ac 2- l-Lys- d-Ala- d-Ala and with Ac 2- l-Lys- l-Ala- l-Ala peptide ligands were calculated. A transition state between the vancomycin complex with the D-form peptide and the complex with the L-form peptide was employed as the perturbation final stage in order to overcome the difficulties of calculating the perturbation between isomers. The calculated relative binding of free energy was 5.15 kcal/mol in water, which compared well with an experimental value of 8.4 kcal/mol. Structure analysis were performed on the number of possible hydrogen bondings, the hydrophobic interactions between peptide ligands and vancomycin, and the bonding angle between designated nuclei. The calculated structure results were all well supported by the NMR data. An investigation of the relative free energy calculation suggests that relative binding affinity in a host-guest system can be achieved by molecular dynamic simulation and the finite difference thermodynamic integration method.

A relative free energy calculation method with ion screening effect using molecular dynamics simulation

A relative free energy calculation method with ion screening effect using molecular dynamics simulation