Free Energy Perturbation Simulations Research Articles

Exploring the C-X…π Halogen Bonding Motif: An Infrared and Raman Study of the Complexes of CF3X (X = Cl, Br and I) with the Aromatic Model Compounds Benzene and Toluene

The formation of halogen bonded complexes formed between the trifluorohalomethanes CF3Cl, CF3Br and CF3I and the Lewis bases benzene and toluene at temperatures below 150K was investigated using FTIR and Raman spectroscopy. Experiments using liquid krypton as solvent show that for both CF3Br and CF3I substantial fractions of the monomers can be involved in 1:1 complexes. In addition, weak absorptions illustrating the formation of 2:1 complexes between CF3I and benzene are observed. Using spectra recorded at temperatures between 120 and 140 K, observed information on the relative stability was obtained for all complexes by determining the complexation enthalpies in solution. The resulting values for CF3Br·benzene, CF3I·benzene and (CF3I)2·benzene are −6.5(3), −7.6(2) and −14.5(9) kJ mol−1. The values for CF3Br·toluene and CF3I·toluene are −6.2(5) and −7.4(5) kJ mol−1. The experimental complexation enthalpies are compared with theoretical data obtained by combining results from MP2/aug-cc-pVDZ(-PP) and MP2/aug-cc-pVTZ(-PP) ab initio calculations, from statistical thermodynamical calculations and from Monte Carlo Free Energy Perturbation simulations. The data are also compared with results derived for other C-X···π halogen bonded complexes involving unsaturated Lewis bases such as ethene and ethyne.

Insights into the Human Glutamate Transporters from the Bacterial Homolog Glt(Ph)

Human glutamate/aspartate transporters cotransport three Na+ and one H+ ions with the substrate and countertransport one K+ ion. The archeal homolog GltPh, which has a 36% sequence identity with the human transporters, transports aspartate and three sodium ions across the lipid membrane. We have previously performed Free Energy Perturbation (FEP) simulations of ligand binding to GltPh in the extracellular(EC)-facing state, showing that the substrate binds after two Na+ ions, with a strong water-mediated coupling between the second sodium ion to bind (Na1) and the substrate. We have also performed simulations and mutagenesis experiments to locate the binding site of the third sodium ion (Na3), which could not be resolved in the crystal structure. More recently, we have turned our focus to the intracellular(IC)-facing conformation of GltPh. By running long unbiased simulations with different ligands bound to the transporter, we were able to describe in details the gating mechanism and point out the differences to the EC-facing state. The IC-facing FEP results demonstrate that the unbinding of the ligands happens in a symmetrical manner, with the ligands being released in the reverse order of their binding. The Na3 ion is, therefore, the first ligand to bind and the last to unbind. By means of free energy calculations, Kramers’ theory and experiments, we are able to point out the unbinding of the Na3 ion as the rate-limiting step in the transport mechanism. With all this information, we have built a homology model for a human glutamate transporter (EAAT1) in different conformations, which can reproduce some of the differences observed in the transport mechanism of the human and the bacterial transporters.

Improving the Prediction of Absolute Solvation Free Energies Using the Next Generation OPLS Force Field.

Explicit solvent molecular dynamics free energy perturbation simulations were performed to predict absolute solvation free energies of 239 diverse small molecules. We use OPLS2.0, the next generation OPLS force field, and compare the results with popular small molecule force fields-OPLS_2005, GAFF, and CHARMm-MSI. OPLS2.0 produces the best correlation with experimental data (R(2) = 0.95, slope = 0.96) and the lowest average unsigned errors (0.7 kcal/mol). Important classes of compounds that performed suboptimally with OPLS_2005 show significant improvements.

On achieving high accuracy and reliability in the calculation of relative protein–ligand binding affinities

We apply a free energy perturbation simulation method, free energy perturbation/replica exchange with solute tempering, to two modifications of protein-ligand complexes that lead to significant conformational changes, the first in the protein and the second in the ligand. The approach is shown to facilitate sampling in these challenging cases where high free energy barriers separate the initial and final conformations and leads to superior convergence of the free energy as demonstrated both by consistency of the results (independence from the starting conformation) and agreement with experimental binding affinity data. The second case, consisting of two neutral thrombin ligands that are taken from a recent medicinal chemistry program for this interesting pharmaceutical target, is of particular significance in that it demonstrates that good results can be obtained for large, complex ligands, as opposed to relatively simple model systems. To achieve quantitative agreement with experiment in the thrombin case, a next generation force field, Optimized Potentials for Liquid Simulations 2.0, is required, which provides superior charges and torsional parameters as compared to earlier alternatives.

Applying efficient implicit nongeometric constraints in alchemical free energy simulations.

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.

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.

A Direct Comparison of the MM-GB/SA Scoring Procedure and Free-Energy Perturbation Calculations Using Carbonic Anhydrase as a Test Case: Strengths and Pitfalls of Each Approach.

MM-GB/SA scoring and free energy perturbation (FEP) calculations have emerged as reliable methodologies to understand structural and energetic relationships to binding. In spite of successful applications to elucidate the structure-activity relationships for few pairs of ligands, the reality is that the performance of FEP calculations has rarely been tested for more than a handful of compounds. In this work, a series of 13 benzene sulfonamide inhibitors of carbonic anhydrase with binding free energies determined by isothermal titration calorimetry was selected as a test case. R(2) values of 0.70, 0.71, and 0.49 with the experiment were obtained with MM-GB/SA and FEP simulations run with MCPRO+ and Desmond, respectively. All methods work well, but the results obtained with Desmond are inferior to MM-GB/SA and MCPRO+. The main contrast between the methods is the level of sampling, ranging from full to restricted flexibility to single conformation for the complexes in Desmond, MCPRO+, and MM-GB/SA, respectively. The current and historical results obtained with MM-GB/SA qualify this approach as a more attractive alternative for rank-ordering; it can achieve equivalent or superior predictive accuracy and handle more structurally dissimilar ligands at a fraction of the computational cost of the rigorous free-energy methods. As for the large theoretical dynamic range for the binding energies, that seems to be a direct result of the degree of sampling in the simulations since MCPRO+ as well as MM-GB/SA are plagued by this. Van't Hoff analysis for selected pairs of ligands suggests that the wider scoring spread is not only affected by missing entropic contributions due to restricted sampling but also exaggerated enthalpic separation between the weak and potent compounds caused by diminished shielding of electrostatic interactions, thermal effects, and protein relaxation/strain.

QM/MM Modeling of Ba 2+ Blockades in Potassium ion Channels

The robust selectivity of potassium channels for K+ over Na+ ions is a major component of the regulation of intracellular K+ concentrations. The selectivity for K+ was quantified through experiments measuring the Na+ and K+ dependence on Ba2+-blockades,(1) indicating that K+ has greater permeability by at least 150 fold. In thermodynamic terms, the relative binding free energy of Na+ to the pore must be at least 3 kcal/mol less favorable than K+. Na+ vs K+ free energy perturbation (FEP) simulations are consistent with this, although no simulations to date have modeled the actual Ba2+ blockade experiment. We have used MD simulations to calculate the relative binding energies of Na+ and K+ in the KcsA ion channel when the S4 site is occupied by Ba2+. As Ba2+ is a strongly polarizing ion, we have used QM/MM FEP calculations using CHARMM interfaced to the deMon DFT code (2), as well as the polarizable Drude force field to correctly model the Ba2+-filter interactions, with the aim of better interpreting the original selectivity experiments.1. J. Neyton, C. Miller, J. Gen. Physiol. 92: 549-567.2. B. Lev et al., J. Comp. Chem., 31: 1015-1023.View Large Image | View Hi-Res Image | Download PowerPoint Slide

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.

Free energy perturbation simulation on transition states and high-activity mutants of human butyrylcholinesterase for (-)-cocaine hydrolysis.

A unified computational approach based on free energy perturbation (FEP) simulations of transition states has been employed to calculate the mutation-caused shifts of the free energy change from the free enzyme to the rate-determining transition state for (-)-cocaine hydrolysis catalyzed by the currently most promising series of mutants of human butyrylcholinesterase (BChE) that contain the A199S/A328W/Y332G mutations. The FEP simulations were followed by Michaelis-Menten kinetics analysis determining the individual k(cat) and K(M) values missing for the A199S/F227A/A328W/Y332G mutant in this series. The calculated mutation-caused shifts of the free energy change from the free enzyme to the rate-determining transition state are in good agreement with the experimental kinetic data, demonstrating that the unified computational approach based on the FEP simulations of the transition states may be valuable for future computational design of new BChE mutants with a further improved catalytic efficiency against (-)-cocaine.

Molecular simulations provide new insights into the role of the accessory immunoglobulin-like domain of Cel9A

Cel9A from the thermoacidophilic bacterium Alicyclobacillus acidocaldarius belongs to the subfamily E1 of family 9 glycoside hydrolases, many members of which have an N-terminal Ig-like domain followed by the catalytic domain. The Ig-like domain is not directly involved in either carbohydrate binding or biocatalysis; however, deletion of the Ig-domain promotes loss of enzymatic activity. We have investigated the functional role of the Ig-like domain using molecular dynamics simulations. Our simulations indicate that residues within the Ig-like domain are dynamically correlated with residues in the carbohydrate-binding pocket and with key catalytic residues of Cel9A. Free energy perturbation simulations indicate that the Ig-like domain stabilizes the catalytic domain and may be responsible for the enhanced thermostability of Cel9A.

Calculation of the standard binding free energy of sparsomycin to the ribosomal peptidyl‐transferase P‐site using molecular dynamics simulations with restraining potentials

The standard (absolute) binding free energy of the antibiotic sparsomycin with the 50S bacteria ribosomal subunit is calculated using molecular dynamics (MD) free energy perturbation (FEP) simulations with restraining potentials developed by Wang et al. [Biophys. J. 91:2798-2814 (2006)]. In the simulation protocol, restraining potentials are activated for the orientational and translational movements of the ligand relative to the binding site when it is decoupled from the binding pocket, and then released once the ligand fully interacts with the rest of the system. A reduced system is simulated to decrease the computational cost of the FEP/MD calculations and the effects of the surrounding atoms are incorporated using the generalized solvent boundary potential (GSBP) method. The loss of conformational freedom of the ligand upon binding is characterized using the potential of mean force (PMF) as a function of the root-mean-square deviation (RMSD) relative to the bound conformation. The number of water molecules in the binding pocket is allowed to fluctuate dynamically in response to the ligand during the calculations by combining FEP/MD with grand canonical Monte Carlo (GCMC) simulations. The calculated binding free energy is about -6 kcal/mol, which is in reasonable agreement with the experimental value. The information gleaned from this study provides new insight on the recognition of ribosome by sparsomycin and highlights the challenges in calculations of absolute binding free energies in these systems.

The Role of the Protein-Conducting Channel in the Membrane Insertion of Transmembrane Segments

In all domains of life, the majority of membrane proteins are inserted into the membrane via a protein-conducting channel, also known as the SecY or Sec61 complex. In addition to a translocation pathway across the membrane, this channel features a unique lateral gate, which can open toward the membrane, permitting the sequential insertion of transmembrane segments (TMs). How this insertion occurs is still unclear, although a thermodynamic partioning between channel and membrane environments has been proposed. However, experiment- and simulation-based scales for the free-energy insertion cost of various amino acids differ, sometimes significantly as in the case of arginine (2-3 kcal/mol in experiment compared to 17 kcal/mol in simulation). Using free energy perturbation (FEP) simulations, we have calculated the insertion cost for an arginine located on a background poly-leucine helix, both in the center of a pure bilayer and in the center of a model of SecY featuring an open lateral gate. We find that the presence of SecY greatly reduces the membrane insertion cost for arginine, in agreement with prior simulations. We also consider the free energy cost for the insertion of the background helix from SecY to the membrane, which had been neglected previously.

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi-particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force-field (CHARMM27) or a polarizable force-field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual-thermostat method and the VV2 integrator to run molecular dynamics simulations.

The energetics of the acetylation switch in p53‐mediated transcriptional activation

Targeted therapeutic intervention in receptor-ligand interactions of p53-mediated tumor suppression can impact progression of disease, aging, and variation in genetic expression. Here, we conducted a number of molecular simulations, based on structures of p53 in complex with its transcriptional coactivating CBP bromodomain, determined by NMR spectroscopy, to investigate the energetics of the binding complex. Building on the observation that acetylation of K382 in p53 serves as the essential triggering switch for a specific interaction with CBP, we assessed the differential effect of acetylation on binding from simulations of an octapeptide derived from p53 with acetylated and nonacetylated K382 (residues 379-386). Cluster analysis of the simulations shows that acetylation of the free peptide does not significantly change the population of the preferred conformation of the peptide in solution for binding to CBP. Conversion of the acetylated K382 to nonacetylated form with free energy perturbation (FEP) simulations of the p53 CBP complex and the free peptide showed that the relative contribution of the acetyl group to binding is 4.8 kcal/mol. An analysis of residue contributions to the binding energy using an MM-GBSA approach agrees with the FEP results and sheds additional light on the origin of selectivity in p53 binding to the CBP bromodomain.

Free-Energy Perturbation Simulation on Transition States and Redesign of Butyrylcholinesterase

It is recognized that an ideal anti-cocaine treatment is to accelerate cocaine metabolism by producing biologically inactive metabolites via a route similar to the primary cocaine-metabolizing pathway, i.e., butyrylcholinesterase (BChE)-catalyzed hydrolysis of cocaine. BChE mutants with a higher catalytic activity against (-)-cocaine are highly desired for use as an exogenous enzyme in humans. To develop a rational design for high-activity mutants, we carried out free-energy perturbation (FEP) simulations on various mutations of the transition-state structures in addition to the corresponding free-enzyme structures by using an extended FEP procedure. The FEP simulations on the mutations of both the free-enzyme and transition-state structures allowed us to calculate the mutation-caused shift of the free-energy change from the free enzyme (BChE) to the transition state, and thus to theoretically predict the mutation-caused shift of the catalytic efficiency (kcat/KM). The computational predictions are supported by the kinetic data obtained from the wet experiments, demonstrating that the FEP-based computational design approach is promising for rational design of high-activity mutants of an enzyme. One of the BChE mutants designed and discovered in this study has an ∼1800-fold improved catalytic efficiency against (-)-cocaine compared to wild-type BChE. The high-activity mutant may be therapeutically valuable.

Charges for Large Scale Binding Free Energy Calculations with the Linear Interaction Energy Method.

The linear interaction energy method (LIE), which combines force field based molecular dynamics (MD) simulations and linear response theory, has previously been shown to give fast and reliable estimates of ligand binding free energies, suggesting that this type of technique could be used also in a high-throughput fashion. However, a limiting step in such applications is the assignment of atomic charges for compounds that have not been parametrized within the given force field, in this case OPLS-AA. In order to reach an automatable solution to this problem, we have examined the performance of nine different ab initio and semiempirical charge methods, together with estimates of solvent induced polarization. A test set of ten HIV-1 reverse transcriptase inhibitors was selected, and LIE estimates of their relative binding free energies were calculated using the resulting 23 different charge variants. Over 800 ns of MD simulation show that the LIE method provides excellent estimates with several different charge methods and that the semiempirically derived CM1A charges, in particular, emerge as a fast and reliable alternative for fully automated LIE based virtual screens with the OPLS-AA force field. Our conclusions regarding different charge models are also expected to be valid for other types of force field based binding free energy calculations, such as free energy perturbation and thermodynamic integration simulations.

Single Mutation Induced H3N2 Hemagglutinin Antibody Neutralization: A Free Energy Perturbation Study

The single mutation effect on the binding affinity of H3N2 viral protein hemagglutinin (HA) with the monoclonical antibody fragment (Fab) is studied in this paper using the free energy perturbation (FEP) simulations. An all-atom protein model with explicit solvents is used to perform an aggregate of several microsecond FEP molecular dynamics simulations. A recent experiment shows that a single mutation in H3N2 HA, T131I, increases the antibody-antigen dissociation constant Kd by a factor of approximately 4000 (equivalent to a binding affinity decrease of approximately 5 kcal/mol), thus introducing an escape of the antibody (Ab) neutralization. Our FEP result confirms this experimental finding by estimating the HA-Ab binding affinity decrease of 5.2 +/- 0.9 kcal/mol but with a somewhat different molecular mechanism from the experimental findings. Detailed analysis reveals that this large binding affinity decrease in the T131I mutant is mainly due to the displacement of two bridge water molecules otherwise present in the wild-type HA/Ab interface. The decomposition of the binding free energy supports this observation, as the major contribution to the binding affinity is from the electrostatic interactions. In addition, we find that the loss of the binding affinity is also related to the large conformational distortion of one loop (loop 155-161) in the unbound state of the mutant. We then simulate all other possible mutations for this specific mutation site T131, and predict a few more mutations with even larger decreases in the binding affinity (i.e., better candidates for antibody neutralization), such as T131W, T131Y, and T131F. As for further validation, we have also modeled another mutation, S157L, with experimental binding affinity available (Kd increasing approximately 500 times), and found a binding affinity decrease of 4.1 +/- 1.0 kcal/mol, which is again in excellent agreement with experiment. These large scale simulations might provide new insights into the detailed physical interaction, possible future escape mutation, and antibody-antigen coevolution relationship between influenza virus and human antibodies.

Optimized Parameters for Continuum Solvation Calculations with Carbohydrates

Continuum solvent models have been shown to be an efficient method for the calculation of the energetics of biomolecules in solution. However, for these methods to produce accurate results, an appropriate set of atomic radii or volumes is needed. While these have been developed for proteins and nucleic acids, the same is not true of carbohydrates. Here, a set of optimized parameters for continuum solvation calculations of carbohydrates in conjunction with the Carbohydrate Solution Force Field are presented. Explicit solvent free-energy perturbation simulations were performed on a set of hexapyranose sugars and used to fit atomic radii for Poisson-Boltzmann and generalized-Born calculations, and to fit atomic volumes for use with the analytical continuum electrostatics model. The solvation energetics computed with the optimized radii and a Poisson-Boltzmann model show remarkable agreement with explicit solvent simulation, with a root-mean-square error of 1.19 kcal/mol over a large test set of sugars in many conformations. The generalized-Born model gives slightly poorer agreement, but still correlates very strongly, with an error of 1.69 kcal/mol. The analytical continuum electrostatics model correlates well with the explicit solvent results, but gives a larger error of 4.71 kcal/mol. The remarkable agreement between the solvation free energies computed in explicit and implicit solvent provides strong motivation for the use of implicit solvent models in the simulation of carbohydrate-containing systems.

Free Energy Perturbation (FEP) Simulation on the Transition States of Cocaine Hydrolysis Catalyzed by Human Butyrylcholinesterase and Its Mutants

A novel computational protocol based on free energy perturbation (FEP) simulations on both the free enzyme and transition state structures has been developed and tested to predict the mutation-caused shift of the free energy change from the free enzyme to the rate-determining transition state for human butyrylcholinesterase (BChE)-catalyzed hydrolysis of (-)-cocaine. The calculated shift, denoted by DeltaDeltaG(1 --> 2), of such kind of free energy change determines the catalytic efficiency (kcat/KM) change caused by the simulated mutation transforming enzyme 1 to enzyme 2. By using the FEP-based computational protocol, the DeltaDeltaG(1 --> 2) values for the mutations A328W/Y332A --> A328W/Y332G and A328W/Y332G --> A328W/Y332G/A199S were calculated to be -0.22 and -1.94 kcal/mol, respectively. The calculated DeltaDeltaG(1 --> 2) values predict that the change from the A328W/Y332A mutant to the A328W/Y332G mutant should slightly improve the catalytic efficiency and that the change from the A328W/Y332G mutant to the A328W/Y332G/A199S mutant should significantly improve the catalytic efficiency of the enzyme for the (-)-cocaine hydrolysis. The predicted catalytic efficiency increases are supported by the experimental data showing that kcat/KM = 8.5 x 10(6), 1.4 x 10(7), and 7.2 x 10(7) min(-1) M(-1) for the A328W/Y332A, A328W/Y332G, and A328W/Y332G/A199S mutants, respectively. The qualitative agreement between the computational and experimental data suggests that the FEP simulations may provide a promising protocol for rational design of high-activity mutants of an enzyme. The general computational strategy of the FEP simulation on a transition state can be used to study the effects of a mutation on the activation free energy for any enzymatic reaction.