Free Energy Perturbation Method Research Articles

Computation of Standard Binding Free Energies of Polar and Charged Ligands to the Glutamate Receptor GluA2

Accurate calculation of the binding affinity of small molecules to proteins has the potential to become an important tool in rational drug design. In this study, we use the free energy perturbation (FEP) method with restraints to calculate the standard binding free energy of five ligands (ACPA, AMPA, CNQX, DNQX, and glutamate) to the glutamate receptor GluA2, which plays an essential role in synaptic transmission. To deal with the convergence problem in FEP calculations with charged ligands, we use a protocol where the ligand is coupled in the binding site while it is decoupled in bulk solution simultaneously. The contributions from the conformational, rotational, and translational entropies to the standard binding free energy are determined by applying/releasing respective restraints to the ligand in bulk/binding site. We also employ the confine-and-release approach, which helps to resolve convergence problems in FEP calculations. Our results are in good agreement with the experimental values for all five ligands, including the charged ones which are often problematic in FEP calculations. We also analyze the different contributions to the binding free energy of each ligand to GluA2 and discuss the nature of these interactions.

A Potent and Selective Peptide Blocker of the Kv1.3 Channel: Prediction from Free-Energy Simulations and Experimental Confirmation

The voltage-gated potassium channel Kv1.3 is a well-established target for treatment of autoimmune diseases. ShK peptide from a sea anemone is one of the most potent blockers of Kv1.3 but its application as a therapeutic agent for autoimmune diseases is limited by its lack of selectivity against other Kv channels, in particular Kv1.1. Accurate models of Kv1.x-ShK complexes suggest that specific charge mutations on ShK could considerably enhance its specificity for Kv1.3. Here we evaluate the K18A mutation on ShK, and calculate the change in binding free energy associated with this mutation using the path-independent free energy perturbation and thermodynamic integration methods, with a novel implementation that avoids convergence problems. To check the accuracy of the results, the binding free energy differences were also determined from path-dependent potential of mean force calculations. The two methods yield consistent results for the K18A mutation in ShK and predict a 2 kcal/mol gain in Kv1.3/Kv1.1 selectivity free energy relative to wild-type peptide. Functional assays confirm the predicted selectivity gain for ShK[K18A] and suggest that it will be a valuable lead in the development of therapeutics for autoimmune diseases.

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree-Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self-consistent field method, Møller-Plesset perturbation theory method, and time-dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self-consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented.

Quantum mechanics/molecular mechanics dual Hamiltonian free energy perturbation

The dual Hamiltonian free energy perturbation (DH-FEP) method is designed for accurate and efficient evaluation of the free energy profile of chemical reactions in quantum mechanical/molecular mechanical (QM/MM) calculations. In contrast to existing QM/MM FEP variants, the QM region is not kept frozen during sampling, but all degrees of freedom except for the reaction coordinate are sampled. In the DH-FEP scheme, the sampling is done by semiempirical QM/MM molecular dynamics (MD), while the perturbation energy differences are evaluated from high-level QM/MM single-point calculations at regular intervals, skipping a pre-defined number of MD sampling steps. After validating our method using an analytic model potential with an exactly known solution, we report a QM/MM DH-FEP study of the enzymatic reaction catalyzed by chorismate mutase. We suggest guidelines for QM/MM DH-FEP calculations and default values for the required computational parameters. In the case of chorismate mutase, we apply the DH-FEP approach in combination with a single one-dimensional reaction coordinate and with a two-dimensional collective coordinate (two individual distances), with superior results for the latter choice.

Clarification on the Decarboxylation Mechanism in KasA Based on the Protonation State of Key Residues in the Acyl-Enzyme State

The β-ketoacyl ACP synthase I (KasA) is a promising drug target because it is essential for the survival of Mycobacterium tuberculosis , a causative agent of tuberculosis. It catalyzes a condensation reaction that comprises three steps. The resulting elongated acyl chains are subsequently needed for the cell wall construction. While the mechanism of the first step (acylation of Cys171 in the active site) is straightforward already, the second step (decarboxylation of malonyl substrate) has been controversial due to the difficulty in determining the correct protonation states of the involved residues (His311, His345, Lys340, Glu354). Available experimental data suggest three possible mechanisms which differ considerably. They are not consistent with each other because these studies could not be performed for KasA at the beginning of decarboxylation step (acyl-enzyme state of KasA). Instead, different mutants had to be used which are expected to resemble this situation. In this first computational study about this topic, we use the free energy perturbation (FEP) method to compute the relevant pKa values in the acyl-enzyme state of KasA and use molecular dynamics (MD) simulations to rationalize the results. Subsequent density functional theory (DFT)-based quantum mechanical/molecular mechanical (QM/MM) MD simulations and umbrella samplings have been used to disentangle the close relationships between the protonation states of the involved residues. By these simulations, we can address the preferred protonation states and roles of the residues involved in decarboxylation reaction, thereby suggesting the possible mechanism for the decarboxylation step.

Advances in Binding Free Energies Calculations: QM/MM-Based Free Energy Perturbation Method for Drug Design

Multiple approaches have been devised and evaluated to computationally estimate binding free energies. Results using a recently developed Quantum Mechanics (QM)/Molecular Mechanics (MM) based Free Energy Perturbation (FEP) method suggest that this method has the potential to provide the most accurate estimation of binding affinities to date. The method treats ligands/inhibitors using QM while using MM for the rest of the system. The method has been applied and validated for a structurally diverse set of fructose 1,6- bisphosphatase (FBPase) inhibitors suggesting that the approach has the potential to be used as an integral part of drug discovery for both lead identification lead optimization, where there is a structure available. In addition, this QM/MM-based FEP method was shown to accurately replicate the anomalous hydration behavior exhibited by simple amines and amides suggesting that the method may also prove useful in predicting physical properties of molecules. While the method is about 5-fold more computationally demanding than conventional FEP, it has the potential to be less demanding on the end user since it avoids development of MM force field parameters for novel ligands and thereby eliminates this time-consuming step that often contributes significantly to the inaccuracy of binding affinity predictions using conventional FEP methods. The QM/MM-based FEP method has been extensively tested with respect to important considerations such as the length of the simulation required to obtain satisfactory convergence in the calculated relative solvation and binding free energies for both small and large structural changes between ligands. Future automation of the method and parallelization of the code is expected to enhance the speed and increase its use for drug design and lead optimization.

Are the Intramolecular O–H···F and O–H···Cl Hydrogen Bonds Maintained in Solution? A Theoretical Study

The present case study aims at calculating the equilibrium conformer compositions for 2X-ethanol and 2X-phenol (X = F, Cl) in solution, and exploring the effect of the applied theoretical method and basis set on the obtained results, as well as considering the usefulness of the continuum solvent approach in comparison with the explicit solvent Monte Carlo model utilizing the free energy perturbation method. Gas-phase optimizations at the DFT/B97D/aug-cc-pvtz and ab initio MP2/aug-cc-pvtz levels predicted structures in good agreement with the available experimental data for three test molecules. Because in-solution geometries change only slightly according to the IEF-PCM continuum solvent calculations in carbon tetrachloride and water, the two theoretical levels were applied further on, and complete basis set (CBS) relative internal free energies were estimated for the conformers under study. The predicted OCCF gauche/trans ratio for 2F-ethanol was well reproduced in comparison with available experimental compositions. The predominant gauche structure maintains an intramolecular hydrogen bond in carbon tetrachloride (HB structure), whereas HB and NoHB gauche conformers appear in nearly the same fraction in aqueous solution. The internally hydrogen-bonded conformer is predominant also for 2X-phenol species, as calculated on the basis of relative CBS internal free energies and IEF-PCM and FEP/MC solvation free energies. Use of a trihydrate supermolecule model for 2F-ethanol conformers leads to the prediction of the aqueous-solution composition in contrast to the experiment. Solution structure modeling predicts weak hydrogen-bond formation capacity for both the covalently bound F and Cl atoms, even in conformations where they are fully exposed to hydration.

Binding of Isoflurane to a Bacterial Voltage-Gated Sodium Channel: Structure and Accessibility of Distinct Interaction Sites

The most likely targets for volatile anesthetics (VA) are ion channels. The mechanism of volatile anesthetic action is not completely understood. Identifying the molecular pathway for drug binding is crucial to understand the effect of VA on voltage gated sodium channels. We use Molecular Dynamics simulations to identify the binding sites for a hydrophobic general anesthetic isoflurane, on voltage gated bacterial sodium channel NaChBac. Apart from the voltage sensing domain (S1-S4), linker and the pore domain (S5-S6), bacterial sodium channels also have fenestrations, which provide a hydrophobic tunnel through the lipid-embedded portion of the channel to the central cavity, where the known local anesthetic site is located. Unbiased “flooding” simulations were performed on the activated open confirmation of NaChBac. We performed a cluster analysis to identify all the possible binding sites of isoflurane. The three most important ones among them are: a region near the selectivity filter, called the extracellular site, a region near the S4-S5 linker, called the linker site, and a region within the cavity, called the cavity site. The most important observation is that isoflurane enters the central cavity through the fenestrations. Free energy perturbation method was employed to calculate the binding affinities of isoflurane for each of these sites. We also studied the interactions between isoflurane and the amino acids in these three binding sites.

Simulation Studies on Nucleotide Selection in T7 RNA Polymerase Transcription

Transcription by RNA polymerases (RNAPs) is essential in gene expression. During transcription elongation, an incoming RNA nucleotide binds to the catalytic center of the polymerase, and is added to the end of an existing transcript when the nucleotide correctly pairs with a template DNA nucleotide. Accurate selection of rNTPs is key to maintaining high transcription fidelity. Detailed mechanisms on how the nucleotides are selected are still lack of. using T7 RNAP as a model system, we calculate binding affinities of various NTP molecules with the RNAP elongation complex, in both its insertion and pre-insertion states. Stable complex structures are obtained by using molecular dynamics simulation. Techniques such as MM-GB/SA, QM/MM-GB/SA are applied to compare the rNTP binding affinities. In addition, free energy perturbation and thermodynamic integration methods are also utilized to calculate the rNTP binding free energies. Based on the simulation results, we would be able to decide which checkpoints (pre-insertion and/or insertion state) the RNAP utilizes essentially to select against wrong nucleotides. We would also like to know how each selection proceeds (rejection and/or inhibition), as well as how strong each selection is in terms of the binding free energy differences between right and wrong nucleotides.

Catalytic Effects of Mutations of Distant Protein Residues in Human DNA Polymerase β: Theory and Experiment

We carried out free-energy calculations and transient kinetic experiments for the insertion of the right (dC) and wrong (dA) nucleotides by wild-type (WT) and six mutant variants of human DNA polymerase β (Pol β). Since the mutated residues in the point mutants, I174S, I260Q, M282L, H285D, E288K, and K289M, were not located in the Pol β catalytic site, we assumed that the WT and its point mutants share the same dianionic phosphorane transition-state structure of the triphosphate moiety of deoxyribonucleotide 5'-triphosphate (dNTP) substrate. On the basis of this assumption, we have formulated a thermodynamic cycle for calculating relative dNTP insertion efficiencies, Ω = (k(pol)/K(D))(mut)/(k(pol)/K(D))(WT) using free-energy perturbation (FEP) and linear interaction energy (LIE) methods. Kinetic studies on five of the mutants have been published previously using different experimental conditions, e.g., primer-template sequences. We have performed a presteady kinetic analysis for the six mutants for comparison with wild-type Pol β using the same conditions, including the same primer/template DNA sequence proximal to the dNTP insertion site used for X-ray crystallographic studies. This consistent set of kinetic and structural data allowed us to eliminate the DNA sequence from the list of factors that can adversely affect calculated Ω values. The calculations using the FEP free energies scaled by 0.5 yielded 0.9 and 1.1 standard deviations from the experimental log Ω values for the insertion of the right and wrong dNTP, respectively. We examined a hybrid FEP/LIE method in which the FEP van der Waals term for the interaction of the mutated amino acid residue with its surrounding environment was replaced by the corresponding van der Waals term calculated using the LIE method, resulting in improved 0.4 and 1.0 standard deviations from the experimental log Ω values. These scaled FEP and FEP/LIE methods were also used to predict log Ω for R283A and R283L Pol β mutants.

Chemically transferable coarse-grained potentials from conditional reversible work calculations

The representability and transferability of effective pair potentials used in multiscale simulations of soft matter systems is ill understood. In this paper, we study liquid state systems composed of n-alkanes, the coarse-grained (CG) potential of which may be assumed pairwise additive and has been obtained using the conditional reversible work (CRW) method. The CRW method is a free-energy-based coarse-graining procedure, which, by means of performing the coarse graining at pair level, rigorously provides a pair potential that describes the interaction free energy between two mapped atom groups (beads) embedded in their respective chemical environments. The pairwise nature of the interactions combined with their dependence on the chemically bonded environment makes CRW potentials ideally suited in studies of chemical transferability. We report CRW potentials for hexane using a mapping scheme that merges two heavy atoms in one CG bead. It is shown that the model is chemically and thermodynamically transferable to alkanes of different chain lengths in the liquid phase at temperatures between the melting and the boiling point under atmospheric (1 atm) pressure conditions. It is further shown that CRW-CG potentials may be readily obtained from a single simulation of the liquid state using the free energy perturbation method, thereby providing a fast and versatile molecular coarse graining method for aliphatic molecules.

ChemInform Abstract: Computation of Kinetic Isotope Effects for Enzymatic Reactions

AbstractReview: 72 refs.

Free-energy molecular simulations of the inclusion complex of Ne with fullerene C60 in water

Binding free energy of Ne in the endohedral fullerene Ne@C60 in water was calculated using closed thermodynamic cycle that involves annihilation in silico of Ne in the free and bound states. The underlying microscopic force field molecular dynamics simulations were performed at 300 K in the constant-volume 27.9 Å periodic box of the TIP3P water molecules. The reference binding free energy of a dummy atom to C60 was calculated from the ideal gas theory equations. The free energies of atom annihilation were computed by means of the free-energy perturbation and Bennett acceptance ratio methods. Introduction of intermediate particle of small size in the mutation scheme in C60 substantially improved free-energy convergence. The predicted binding free energies of Ne and Ne cavity to C60 at the standard state concentration are about − 5 kJ mol− 1 and +22 kJ mol− 1, respectively. The free energy of reversible annihilation of Ne is controlled by the insertion step for simulations both in water and in C60.

Free-Energy Simulations Reveal that Both Hydrophobic and Polar Interactions Are Important for Influenza Hemagglutinin Antibody Binding

Antibodies binding to conserved epitopes can provide a broad range of neutralization to existing influenza subtypes and may also prevent the propagation of potential pandemic viruses by fighting against emerging strands. Here we propose a computational framework to study structural binding patterns and detailed molecular mechanisms of viral surface glycoprotein hemagglutinin (HA) binding with a broad spectrum of neutralizing monoclonal antibody fragments (Fab). We used rigorous free-energy perturbation (FEP) methods to calculate the antigen-antibody binding affinities, with an aggregate underlying molecular-dynamics simulation time of several microseconds (∼2 μs) using all-atom, explicit-solvent models. We achieved a high accuracy in the validation of our FEP protocol against a series of known binding affinities for this complex system, with <0.5 kcal/mol errors on average. We then introduced what to our knowledge are novel mutations into the interfacial region to further study the binding mechanism. We found that the stacking interaction between Trp-21 in HA2 and Phe-55 in the CDR-H2 of Fab is crucial to the antibody-antigen association. A single mutation of either W21A or F55A can cause a binding affinity decrease of ΔΔG > 4.0 kcal/mol (equivalent to an ∼1000-fold increase in the dissociation constant Kd). Moreover, for group 1 HA subtypes (which include both the H1N1 swine flu and the H5N1 bird flu), the relative binding affinities change only slightly (< ±1 kcal/mol) when nonpolar residues at the αA helix of HA mutate to conservative amino acids of similar size, which explains the broad neutralization capability of antibodies such as F10 and CR6261. Finally, we found that the hydrogen-bonding network between His-38 (in HA1) and Ser-30/Gln-64 (in Fab) is important for preserving the strong binding of Fab against group 1 HAs, whereas the lack of such hydrogen bonds with Asn-38 in most group 2 HAs may be responsible for the escape of antibody neutralization. These large-scale simulations may provide new insight into the antigen-antibody binding mechanism at the atomic level, which could be essential for designing more-effective vaccines for influenza.

Achieving fast convergence of ab initio free energy perturbation calculations with the adaptive force-matching method

This paper studies the possibility of improving the convergence of ab initio free energy perturbation (FEP) calculations by developing customized force fields with the adaptive force-matching (AFM) method. The ab initio FEP method relies on a molecular mechanics (MM) potential to sample configuration space. If the Boltzmann weight of the MM sampling is close to that of the ab initio method, the efficiency of ab initio FEP will be optimal. The difference in the Boltzmann weights can be quantified by the relative energy difference distribution (REDD). The force field developed through AFM significantly improves the REDD when compared with standard MM models, thus improving the convergence of the ab initio FEP calculation. The static dielectric constant es of ice-Ih was studied with PW-91 through ab initio FEP. With a customized force field developed through AFM, we were able to converge es to 80 ± 4 with 3,600 configurations. A similar ab initio FEP calculation with the TIP4P model would require 220 times more configurations to achieve the same accuracy. Our study indicates that the PW-91 functional underestimates ice-Ih es by about 20%.

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.

Computation of the binding affinities of catechol‐O‐methyltransferase inhibitors: Multisubstate relative free energy calculations

Alchemical free energy simulations are amongst the most accurate techniques for the computation of the free energy changes associated with noncovalent protein-ligand interactions. A procedure is presented to estimate the relative binding free energies of several ligands to the same protein target where multiple, low-energy configurational substates might coexist, as opposed to one unique structure. The contributions of all individual substates were estimated, explicitly, with the free energy perturbation method, and combined in a rigorous fashion to compute the overall relative binding free energies and dissociation constants. It is shown that, unless the most stable bound forms are known a priori, inaccurate results may be obtained if the contributions of multiple substates are ignored. The method was applied to study the complex formed between human catechol-O-methyltransferase and BIA 9-1067, a newly developed tight-binding inhibitor that is currently under clinical evaluation for the therapy of Parkinson's disease. Our results reveal an exceptionally high-binding affinity (K(d) in subpicomolar range) and provide insightful clues on the interactions and mechanism of inhibition. The inhibitor is, itself, a slowly reacting substrate of the target enzyme and is released from the complex in the form of O-methylated product. By comparing the experimental catalytic rate (k(cat)) and the estimated dissociation rate (k(off)) constants of the enzyme-inhibitor complex, one can conclude that the observed inhibition potency (K(i)) is primarily dependent on the catalytic rate constant of the inhibitor's O-methylation, rather than the rate constant of dissociation of the complex.

Theoretical studies of HIV-1 reverse transcriptase inhibition

Computational methods for accurately calculating the binding affinity of a ligand for a protein play a pivotal role in rational drug design. We herein present a theoretical study of the binding of five different ligands to one of the proteins responsible for the human immunodeficiency virus type 1 (HIV-1) cycle replication; the HIV-1 reverse transcriptase (RT). Two types of approaches are used based on molecular dynamics (MD) simulations within hybrid QM/MM potentials: the alchemical free energy perturbation method, FEP, and the pathway method, in which the ligand is physically pulled away from the binding site, thus rendering a potential of mean force (PMF) for the binding process. Our comparative analysis stresses their advantages and disadvantages and, although the results are not in quantitative agreement, both methods are capable of distinguishing the most and the less potent inhibitors of HIV-1 RT activity on an RNase H site. The methods can then be used to select the proper scaffold to design new drugs. A deeper analysis of these inhibitors through molecular electrostatic potential (MEP) and calculation of the binding contribution of the individual residues shows that, in a rational design, apart from the strong interactions established with the two magnesium cations present in the RNase H site, it is important to take into account interactions with His539 and with those residues that are anchoring the metals; Asp443, Glu478, Asp498 and Asp549. The MEPs of the active site of the protein and the different ligands show a better complementarity in those inhibitors that present higher binding energies, but there are still possibilities of improving the favourable interactions and decreasing those that are repulsive in order to design compounds with higher inhibitory activity.

Contribution of both catalytic constant and Michaelis constant to CALB enantioselectivity: Use of FEP calculations for prediction studies

Candida antarctica lipase B (CALB) is characterized by its stability and ease of production and is widely used in the pharmaceutical industry. Here we report on the enantioselectivity of the enzyme using both experimental and computational methods. The apparent kinetic parameters were first experimentally determined for enantiopure butan-2-ol and pentan-2-ol substrates. We demonstrate that enantiopreference for the R form of butan-2-ol arises mainly from a lower apparent KM. This corresponds to a major contribution of ΔΔGES, the free energy difference between the ES complex formed with the R and S enantiomers, to ΔΔG‡, the free energy difference between both transition states, in comparison with ΔΔGkcat, the activation free energy difference. In the case of pentan-2-ol, we show that the enantiopreference for the R form comes from both a lower KM and a higher kcat. In addition, we used, for the first time, the free energy perturbation method to evaluate the free energy difference between tetrahedral intermediates formed with R and S alcohol enantiomers for a series of secondary alcohols. This is a valid model for ΔΔG‡. Computational results were found to be in qualitative agreement with experimental data, and enable the determination of substrate orientation in the active site with fair confidence.

Computation of kinetic isotope effects for enzymatic reactions

We describe a computational approach, incorporating quantum mechanics into enzyme kinetics modeling with a special emphasis on computation of kinetic isotope effects. Two aspects are highlighted: (1) the potential energy surface is represented by a combined quantum mechanical and molecular mechanical (QM/MM) potential in which the bond forming and breaking processes are modeled by electronic structure theory, and (2) a free energy perturbation method in path integral simulation is used to determine both kinetic isotope effects (KIEs). In this approach, which is called the PI-FEP/UM method, a light (heavy) isotope is mutated into a heavy (light) counterpart in centroid path integral simulations. The method is illustrated in the study of primary and secondary KIEs in two enzyme systems. In the case of nitroalkane oxidase, the enzymatic reaction exhibits enhanced quantum tunneling over that of the uncatalyzed process in water. In the dopa delarboxylase reaction, there appears to be distinguishable primary carbon-13 and secondary deuterium KIEs when the internal proton tautomerism is in the N-protonated or in the O-protonated positions. These examples show that the incorporation of quantum mechanical effects in enzyme kinetics modeling offers an opportunity to accurately and reliably model the mechanisms and free energies of enzymatic reactions.