Computation of Absolute Hydration and Binding Free Energy with Free Energy Perturbation Distributed Replica-Exchange Molecular Dynamics (FEP/REMD)

Distributed Replica (REPDSTR) is a powerful parallelization technique enabling simulations of a group of replicas in a parallel/parallel fashion, where each replica is distributed to different nodes of a large cluster [Theor. Chem. Acc. 109: 140 (2003)]. Here, we use the framework provided by REPDSTR to combine a staged free energy perturbation protocol with different values of the thermodynamic coupling parameters with replica-exchange molecular dynamics (FEP/REMD. The structure of REPDSTR, which allows multiple parallel input/output (I/O), facilitates the treatment of replica-exchange to couple the N window simulations. As a result, each of the N synchronous window simulations benefit from the sampling carried out by the N-1 others. As illustrative examples of the FEP/REMD strategy, calculations of the absolute hydration and binding free energy of small molecules were performed using the biomolecular simulation program CHARMM adapted for the IBM Blue Gene/P platform. The computations show that a FEP/REMD strategy significantly improves the sampling and accelerate the convergence of absolute free energy computations.     

X‐Ray Crystallography and Free Energy Calculations Reveal the Binding Mechanism of A2A Adenosine Receptor Antagonists

We present a robust protocol based on iterations of free energy perturbation (FEP) calculations, chemical synthesis, biophysical mapping and X‐ray crystallography to reveal the binding mode of an antagonist series to the A_2A adenosine receptor (AR). Eight A_2AAR binding site mutations from biophysical mapping experiments were initially analyzed with sidechain FEP simulations, performed on alternate binding modes. The results distinctively supported one binding mode, which was subsequently used to design new chromone derivatives. Their affinities for the A_2AAR were experimentally determined and investigated through a cycle of ligand‐FEP calculations, validating the binding orientation of the different chemical substituents proposed. Subsequent X‐ray crystallography of the A_2AAR with a low and a high affinity chromone derivative confirmed the predicted binding orientation. The new molecules and structures here reported were driven by free energy calculations, and provide new insights on antagonist binding to the A_2AAR, an emerging target in immuno‐oncology.     

X-Ray Crystallography and Free Energy Calculations Reveal the Binding Mechanism of A2A Adenosine Receptor Antagonists

We present a robust protocol based on iterations of free energy perturbation (FEP) calculations, chemical synthesis, biophysical mapping and X-ray crystallography to reveal the binding mode of an antagonist series to the A_2A adenosine receptor (AR). Eight A_2AAR binding site mutations from biophysical mapping experiments were initially analyzed with sidechain FEP simulations, performed on alternate binding modes. The results distinctively supported one binding mode, which was subsequently used to design new chromone derivatives. Their affinities for the A_2AAR were experimentally determined and investigated through a cycle of ligand-FEP calculations, validating the binding orientation of the different chemical substituents proposed. Subsequent X-ray crystallography of the A_2AAR with a low and a high affinity chromone derivative confirmed the predicted binding orientation. The new molecules and structures here reported were driven by free energy calculations, and provide new insights on antagonist binding to the A_2AAR, an emerging target in immuno-oncology.     

Free energy perturbation calculations on parallel computers: Demonstrations of scalable linear speedup

A coarse-grain parallel implementation of the free energy perturbation (FEP) module of the AMBER molecular dynamics program is described and then demonstrated using five different molecular systems. The difference in the free energy of (aqueous) solvation is calculated for two monovalent cations ΔΔG_aq(Li⁺ Δ Cs⁺), and for the zero-sum ethane-to-ethane′ perturbation ΔΔG_aq(CH₃? methyl? X → X? methyl? CH₃), where X is a ghost methyl. The difference in binding free energy for a docked HIV-1 protease inhibitor into its ethylene mimetic is examined by mutating its fifth peptide bond, ΔG(CO? NH → CH?CH). A potassium ion (K⁺) is driven outward from the center of mass of ionophore salinomycin (SAL^?) in a potential of mean force calculation ΔG_MeOH(SAL^? · K⁺) carried out in methanol solvent. Parallel speedup obtained is linearly proportional to the number of parallel processors applied. Finally, the difference in free energy of solvation of phenol versus benzene, ΔΔG_oct(phenol → benzene), is determined in water-saturated octanol and then expressed in terms of relative partition coefficients, Δ log(P_o/w). Because no interprocessor communication is required, this approach is scalable and applicable in general for any parallel architecture or network of machines. FEP calculations run on the nCUBE/2 using 50 or 100 parallel processors were completed in clock times equivalent to or twice as fast as a Cray Y-MP. The difficulty of ensuring adequate system equilibrium when agradual configurational reorientation follows the mutation of the Hamiltonian is discussed and analyzed. The results of a successful protocol for overcoming this equilibration problem are presented. The types of molecular perturbations for which this method is expected to perform most efficiently are described. © 1994 by John Wiley & Sons, Inc.     

Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations

Free Energy Perturbation with Replica Exchange Molecular Dynamics (FEP/REMD) offers a powerful strategy to improve the convergence of free energy computations. In particular, it has been shown previously that a FEP/REMD scheme allowing random moves within an extended replica ensemble of thermodynamic coupling parameters "lambda" can improve the statistical convergence in calculations of absolute binding free energy of ligands to proteins [J. Chem. Theory Comput. 2009, 5, 2583]. In the present study, FEP/REMD is extended and combined with an accelerated MD simulations method based on Hamiltonian replica-exchange MD (H-REMD) to overcome the additional problems arising from the existence of kinetically trapped conformations within the protein receptor. In the combined strategy, each system with a given thermodynamic coupling factor lambda in the extended ensemble is further coupled with a set of replicas evolving on a biased energy surface with boosting potentials used to accelerate the inter-conversion among different rotameric states of the side chains in the neighborhood of the binding site. Exchanges are allowed to occur alternatively along the axes corresponding to the thermodynamic coupling parameter lambda and the boosting potential, in an extended dual array of coupled lambda- and H-REMD simulations. The method is implemented on the basis of new extensions to the REPDSTR module of the biomolecular simulation program CHARMM. As an illustrative example, the absolute binding free energy of p-xylene to the nonpolar cavity of the L99A mutant of T4 lysozyme was calculated. The tests demonstrate that the dual lambda-REMD and H-REMD simulation scheme greatly accelerates the configurational sampling of the rotameric states of the side chains around the binding pocket, thereby improving the convergence of the FEP computations.     

Use of quantum mechanics/molecular mechanics-based FEP method for calculating relative binding affinities of FBPase inhibitors for type-2 diabetes

A quantum mechanics (QM)/molecular mechanics (MM)-based free energy perturbation (FEP) method, developed recently, provides most accurate estimation of binding affinities. The validity of the method was evaluated for a large set of diverse inhibitors for fructose 1,6-bisphosphatase (FBPase), a target enzyme for type-II diabetes mellitus. The validation set comprises of 22 important structurally different mutations. The calculated relative binding free energies using the QM/MM-based FEP method reproduce the experimental values with exceptional precision of less than ±0.5 kcal/mol. The CPU requirements for QM/MM-based FEP are about fivefold greater than conventional FEP methods, but it is superior in accuracy of predictions. In addition, the QM/MM-based FEP method eliminates the need for time-consuming development of MM force field parameters, which are frequently required for novel inhibitors described by MM. Future automation of the method and parallelization of the code for 128/256/512 cluster computers is expected to enhance the speed and increase its use for drug design and lead optimization. The present application of QM/MM-based FEP method for structurally diverse set of analogs serves to enhance the scope of FEP method and demonstrate the utility of QM/MM-based FEP method for its potential in drug discovery.     

Binding free energies in the SAMPL6 octa-acid host–guest challenge calculated with MM and QM methods

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10–20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4–5.2 kJ/mol and a correlation coefficient (R²) of 0.77–0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD?=?3–8 kJ/mol and R²?=?0.1–0.9), but the results were improved compared to previous studies of this system with similar methods.     

Free Energy Peturbation and Molecular Dynamics Simulation Studies on the Enantiomeric Discrimination of Amines by Dimethyldiketopyridino-18-Crown-6

Abstract

Discrimination of chiral amines by dimethyldiketopyridino-18-crown-6 (1) is studied by free energy peturbation (FEP) and molecular dynamics (MD) methods. 1 has two (S)-chiral centers and discriminates chiral amines through host-guest interactions. The optically active amines in this study are α-(1-naphthyl)ethylamine, methylbenzylamine, cyclohexylethylamine, and sec-butylamine. The trends in binding free energy differences obtained from FEP calculations were in excellent agreement with experimental results obtained in the gas phase. In order to explain the enantioselectivity of the host in terms of the host-guest interactions at the molecular level, we analyzed the structures generated by 10-ns MD simulations of host-guest complexes. The suggested chiral discrimination mechanism, the π-π interaction and the steric repulsion between the guest and the host, was verified by our MD simulation analysis.     

Large scale relative protein ligand binding affinities using non-equilibrium alchemy

Ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design. However, their broad uptake and impact is held back by the notoriously complex setup of the calculations. Only a few tools other than the free energy perturbation approach by Schrödinger Inc. (referred to as FEP+) currently enable end-to-end application. Here, we present for the first time an approach based on the open-source software pmx that allows to easily set up and run alchemical calculations for diverse sets of small molecules using the GROMACS MD engine. The method relies on theoretically rigorous non-equilibrium thermodynamic integration (TI) foundations, and its flexibility allows calculations with multiple force fields. In this study, results from the Amber and Charmm force fields were combined to yield a consensus outcome performing on par with the commercial FEP+ approach. A large dataset of 482 perturbations from 13 different protein–ligand datasets led to an average unsigned error (AUE) of 3.64 ± 0.14 kJ mol⁻¹, equivalent to Schrödinger''s FEP+ AUE of 3.66 ± 0.14 kJ mol⁻¹. For the first time, a setup is presented for overall high precision and high accuracy relative protein–ligand alchemical free energy calculations based on open-source software.

Relative ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design.     

Calculation of relative binding affinities of fructose 1,6-bisphosphatase mutants with adenosine monophosphate using free energy perturbation method

The free energy perturbation (FEP) methodology is the most accurate means of estimating relative binding affinities between inhibitors and protein variants. In this article, the importance of hydrophobic and hydrophilic residues to the binding of adenosine monophosphate (AMP) to the fructose 1,6-bisphosphatase (FBPase), a target enzyme for type-II diabetes, was examined by FEP method. Five mutations were made to the FBPase enzyme with AMP inhibitor bound: 113Tyr --> 113Phe, 31Thr --> 31Ala, 31Thr --> 31Ser, 177Met --> 177Ala, and 30Leu --> 30Phe. These mutations test the strength of hydrogen bonds and van der Waals interactions between the ligand and enzyme. The calculated relative free energies indicated that: 113Tyr and 31Thr play an important role, each via two hydrogen bonds affecting the binding affinity of inhibitor AMP to FBPase, and any changes in these hydrogen bonds due to mutations on the protein will have significant effect on the binding affinity of AMP to FBPase, consistent to experimental results. Also, the free energy calculations clearly show that the hydrophilic interactions are more important than the hydrophobic interactions of the binding pocket of FBPase.     

Computing Alchemical Free Energy Differences with Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) Simulations

Alchemical free energy calculations play a very important role in the field of molecular modeling. Efforts have been made to improve the accuracy and precision of those calculations. One of the efforts is to employ a Hamiltonian replica exchange molecular dynamics (H-REMD) method to enhance conformational sampling. In this paper, we demonstrated that HREMD method not only improves convergence in alchemical free energy calculations but also can be used to compute free energy differences directly via the Free Energy Perturbation (FEP)algorithm. We show a direct mapping between the H-REMD and the usual FEP equations, which are then used directly to compute free energies. The H-REMD alchemical free energy calculation (Replica exchange Free Energy Perturbation, REFEP) was tested on predicting the pK(a) value of the buried Asp26 in thioredoxin. We compare the results of REFEP with TI and regular FEP simulations. REFEP calculations converged faster than those from TI and regular FEP simulations. The final predicted pK(a) value from the H-REMD simulation was also very accurate, only 0.4 pK(a) unit above the experimental value. Utilizing the REFEP algorithm significantly improves conformational sampling, and this in turn improves the convergence of alchemical free energy simulations.     

Free energy perturbation calculations on models of active sites: Applications to adenosine deaminase inhibitors

Free energy perturbation calculations were performed to determine the free energy of binding associated with the presence of perhaps an unusual hydroxyl group in the transition state analog of nebularine, an inhibitor of the enzyme adenosine deaminase. The presence of a single hydroxyl group in this inhibitor has been found to contribute ?9.8 kcal/mol to the free energy of binding, with a 10⁸-fold increase in the binding affinity by the enzyme. In this work, we calculate the difference in solvation free energy for the 1,6-dihydropurine complex versus that of the 6-hydroxyl-1,6-dihydropurine complex to determine if this marked increase in binding affinity is attributed to an unusually hydrophobic hydroxyl group. The calculated ΔG associated for the solvation free energy is ?11.8 kcal/mol. This large change in the solvation free energy suggests that this hydroxyl is instead unusually hydrophilic and that the difference in free energy of interaction for the two inhibitors to the enzyme must be at least ca. 20 kcal/mol. Although the crystal structure for adenosine deaminase is currently not known, we attempt to mimic the nature of the active site by constructing models which simulate the enzyme-inhibitor complex. We present a first attempt at determining the change in free energy of binding for a system in which structural data for the enzyme is incomplete. To do this, we construct what we believe is a minimal model of the binding between adenosine deaminase and an inhibitor. The active site is simulated as a single charged carboxyl group which can form a hydrogen bond with the hydroxyl group of the analog. Two different carboxyl anion models are used. In the first model, the association is modeled between an acetic acid anion and the modified inhibitor. The second model consists of a hydrophobic amino acid pocket with an interior Glu residue in the active site. From these models we calculate the change in free energy of association and the overall change in free energy of binding. We calculate the free energies of interaction both in the absence and presence of water. We conclude from this that the presence of a single suitably placed-CO^?₂ group probably cannot explain the binding effect of the-OH group and that additional interactions will be found in the adenosine deaminase active site.     

Solvent models for protein–ligand binding: Comparison of implicit solvent poisson and surface generalized born models with explicit solvent simulations

Solvent effects play a crucial role in mediating the interactions between proteins and their ligands. Implicit solvent models offer some advantages for modeling these interactions, but they have not been parameterized on such complex problems, and therefore, it is not clear how reliable they are. We have studied the binding of an octapeptide ligand to the murine MHC class I protein using both explicit solvent and implicit solvent models. The solvation free energy calculations are more than 10³ faster using the Surface Generalized Born implicit solvent model compared to FEP simulations with explicit solvent. For some of the electrostatic calculations needed to estimate the binding free energy, there is near quantitative agreement between the explicit and implicit solvent model results; overall, the qualitative trends in the binding predicted by the explicit solvent FEP simulations are reproduced by the implicit solvent model. With an appropriate choice of reference system based on the binding of the discharged ligand, electrostatic interactions are found to enhance the binding affinity because the favorable Coulomb interaction energy between the ligand and protein more than compensates for the unfavorable free energy cost of partially desolvating the ligand upon binding. Some of the effects of protein flexibility and thermal motions on charging the peptide in the solvated complex are also considered. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 591–607, 2001     

Relative binding affinities of alkali metal cations to

The binding affinity and selectivity of a new ionophore, [1(8)]starand (1), toward alkali metal cations in methanol were examined through NMR titration experiments and free energy perturbation (FEP) and molecular dynamics simulations. The preference was determined to be K(+) > Rb(+) > Cs(+) > Na(+) > Li(+) in both FEP simulations and NMR experiments. The FEP simulation results were able to predict the relative binding free energies with errors less than 0.13 kcal/mol, except for the case between Li(+) and Na(+). The cation selectivity was rationalized by analyzing the radial distribution functions of the M-O and M-C distances of free metal cations in methanol and those of metal-ionophore complexes in methanol.     

On Restraints in End-Point Protein–Ligand Binding Free Energy Calculations

The impact of harmonic restraints on protein heavy atoms and ligand atoms on end-point free energy calculations is systematically characterized for 54 protein–ligand complexes. We observe that stronger restraints reduce the equilibration time and statistical inefficiency, suppress conformational sampling, influence correlation with experiment, and monotonically decrease the estimated loss of entropy upon binding, leading to stronger estimated binding free energies in most systems. A statistical estimator that reweights for the biasing potential and includes data prior to the estimated equilibration time has the highest correlation with experiment. A spring constant of 20 cal mol⁻¹ Å⁻² maintains a near-native energy landscape and suppresses artifactual energy minima while minimally limiting thermal fluctuations about the crystal structure. © 2019 Wiley Periodicals, Inc.     

Analysis of drug interactions with very low density lipoprotein by high-performance affinity chromatography

High-performance affinity chromatography (HPAC) was utilized to examine the binding of very low density lipoprotein (VLDL) with drugs, using R/S-propranolol as a model. These studies indicated that two mechanisms existed for the binding of R- and S-propranolol with VLDL. The first mechanism involved non-saturable partitioning of these drugs with VLDL, which probably occurred with the lipoprotein’s non-polar core. This partitioning was described by overall affinity constants of 1.2 (±0.3)?×?10⁶ M^?1 for R-propranolol and 2.4 (±0.6)?×?10⁶ M^?1 for S-propranolol at pH 7.4 and 37 °C. The second mechanism occurred through saturable binding by these drugs at fixed sites on VLDL, such as represented by apolipoproteins on the surface of the lipoprotein. The association equilibrium constants for this saturable binding at 37 °C were 7.0 (±2.3)?×?10⁴ M^?1 for R-propranolol and 9.6 (±2.2)?×?10⁴ M^?1 for S-propranolol. Comparable results were obtained at 20 and 27 °C for the propranolol enantiomers. This work provided fundamental information on the processes involved in the binding of R- and S-propranolol to VLDL, while also illustrating how HPAC can be used to evaluate relatively complex interactions between agents such as VLDL and drugs or other solutes.     

Identification and analysis of stereoselective drug interactions with low-density lipoprotein by high-performance affinity chromatography

Columns containing immobilized low-density lipoprotein (LDL) were prepared for the analysis of drug interactions with this agent by high-performance affinity chromatography (HPAC). R/S-Propranolol was used as a model drug for this study. The LDL columns gave reproducible binding to propranolol over 60 h of continuous use in the presence of pH 7.4 0.067 M potassium phosphate buffer. Experiments conducted with this type of column through frontal analysis indicated that two types of interactions were occurring between R-propranolol and LDL, while only a single type of interaction was observed between S-propranolol and LDL. The first type of interaction, which was seen for both enantiomers, involved non-saturable binding; this interaction had an overall affinity (nK _a) of 1.9 (±0.1) × 10⁵ M⁻¹ for R-propranolol and 2.7 (±0.2) × 10⁵ M⁻¹ for S-propranolol at 37 °C. The second type of interaction was observed only for R-propranolol and involved saturable binding that had an association equilibrium constant (K _a) of 5.2 (±2.3) × 10⁵ M⁻¹ at 37 °C. Similar differences in binding behavior were found for the two enantiomers at 20 °C and 27 °C. This is the first known example of stereoselective binding of drugs by LDL or other lipoproteins. This work also illustrates the ability of HPAC to be used as a tool for characterizing mixed-mode interactions that involve LDL and related binding agents.     

烷基芳基磺酸盐的分子动力学模拟与自由能微扰计算

为研究不同结构的表面活性剂分子在溶液中胶束化能力的差异, 采用分子动力学方法模拟三种烷基芳基磺酸盐在真空和水溶液环境下的结构与相互作用. 利用自由能微扰(FEP)方法计算了水合自由能, 发现与用传统热力学表面张力法测定自制的烷基芳基磺酸盐结果一致. 研究表明: 烷基芳基磺酸盐在水溶液中的胶束化过程是自发进行的, 随着分子结构中芳环向长烷基链中间位置移动, 胶束化能力和胶束稳定性均下降; 疏水基周围水分子的“冰山结构”会影响胶束的稳定性, 而水分子中氢键的生存周期是反映冰山结构变化的重要指标; 同时, 亲水基与水分子间形成氢键的数目会增强或减弱分子脱离胶束体的趋势, 从而影响胶束结构的稳定性.     

Radical heterogeneous polymerization behavior of N‐methyl‐α‐fluoroacrylamide in benzene

The polymerization of N‐methyl‐α‐fluoroacrylamide (NMFAm) initiated with dimethyl 2,2′‐azobisisobutyrate (MAIB) in benzene was studied kinetically and with electron spin resonance. The polymerization proceeded heterogeneously with the highly efficient formation of long‐lived poly(NMFAm) radicals. The overall activation energy of the polymerization was 111 kJ/mol. The polymerization rate (R_p) at 50 °C is given by R_p = k[MAIB]^0.75±0.05 [NMFAm]^0.44±0.05. The concentration of the long‐lived polymer radical increased linearly with time. The formation rate (R_p?) of the long‐lived polymer radical at 50 °C is expressed by R_p? = k[MAIB]^1.0±0.1 [NMFAm]^0±0.1. The overall activation energy of the long‐lived radical formation was 128 kJ/mol, which agreed with the energy of initiation (129 kJ/mol), which was separately estimated. A comparison of R_p? with the initiation rate led to the conclusion that 1‐methoxycarbonyl‐1‐methylethyl radicals (primary radicals from MAIB), escaping from the solvent cage, were quantitatively converted into the long‐lived poly(NMFAm) radicals. Thus, this polymerization involves completely unimolecular termination due to polymer radical occlusion. ¹H NMR‐determined tacticities of resulting poly(NMFAm) were estimated to be rr = 0.34, mr = 0.48, and mm = 0.18. The copolymerization of NMFAm(M₁) and St(M₂) with MAIB at 50 °C in benzene gave monomer reactivity ratios of r₁ = 0.61 and r₂ = 1.79. © 2001 John Wiley & Sons, Inc. J Polym Sci Part A: Polym Chem 39: 2196–2205, 2001