Relative Binding Free Energies Research Articles

Why Mercury Prefers Soft Ligands

Mercury (Hg) is a major global pollutant arising from both natural and anthropogenic sources. Defining the factors that determine the relative affinities of different ligands for the mercuric ion, Hg2+, is critical to understanding its speciation, transformation, and bioaccumulation in the environment. Here, we use quantum chemistry to dissect the relative binding free energies for a series of inorganic anion complexes of Hg2+. Comparison of Hg2+–ligand interactions in the gaseous and aqueous phases shows that differences in interactions with a few, local water molecules led to a clear periodic trend within the chalcogenide and halide groups and resulted in the well-known experimentally observed preference of Hg2+ for soft ligands such as thiols. Our approach establishes a basis for understanding Hg speciation in the biosphere.

Water Network Perturbation in Ligand Binding: Adenosine A2AAntagonists as a Case Study

Recent efforts in the computational evaluation of the thermodynamic properties of water molecules have resulted in the development of promising new in silico methods to evaluate the role of water in ligand binding. These methods include WaterMap, SZMAP, GRID/CRY probe, and Grand Canonical Monte Carlo simulations. They allow the prediction of the position and relative free energy of the water molecule in the protein active site and the analysis of the perturbation of an explicit water network (WNP) as a consequence of ligand binding. We have for the first time extended these approaches toward the prediction of kinetics for small molecules and of relative free energy of binding with a focus on the perturbation of the water network and application to large diverse data sets. Our results support a qualitative correlation between the residence time of 12 related triazine adenosine A(2A) receptor antagonists and the number and position of high energy trapped solvent molecules. From a quantitative viewpoint, we successfully applied these computational techniques as an implicit solvent alternative, in linear combination with a molecular mechanics force field, to predict the relative ligand free energy of binding (WNP-MMSA). The applicability of this linear method, based on the thermodynamics additivity principle, did not extend to 375 diverse A(2A) receptor antagonists. However, a fast but effective method could be enabled by replacing the linear approach with a machine learning technique using probabilistic classification trees, which classified the binding affinity correctly for 90% of the ligands in the training set and 67% in the test set.

Estimation of relative binding free energy based on a free energy variational principle for the FKBP-ligand system

Predicting an accurate binding free energy between a target protein and a ligand can be one of the most important steps in a drug discovery process. Often, many molecules must be screened to find probable high potency ones. Thus, a computational technique with low cost is highly desirable for the estimation of binding free energies of many molecules. Several techniques have thus far been developed for estimating binding free energies. Some techniques provide accurate predictions of binding free energies but high large computational cost. Other methods give good predictions but require tuning of some parameters to predict them with high accuracy. In this study, we propose a method to predict relative binding free energies with accuracy comparable to the results of prior methods but with lower computational cost and with no parameter needing to be carefully tuned. Our technique is based on the free energy variational principle. FK506 binding protein (FKBP) with 18 ligands is taken as a test system. Our results are compared to those from other widely used techniques. Our method provides a correlation coefficient (r²) of 0.80 between experimental and calculated relative binding free energies and yields an average absolute error of 0.70 kcal/mol compared to experimental values. These results are comparable to or better than results from other techniques. We also discuss the possibility to improve our method further.

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.

Quantum polarized ligand docking investigation to understand the significance of protonation states in histone deacetylase inhibitors

The effects of different protonation states of the hydroxamic acid (HA) inhibitors against the class I histone deacetylase enzymes (HDACs) have been studied using the state of the art quantum polarized ligand docking (QPLD) and molecular mechanics-generalized Born surface area (MM-GBSA) approaches. The binding modes of the inhibitors and their inter-molecular interactions with class I HDACs, in response to the protonation states of the inhibitors, are explored. Our results indicate that the different protonation states of the inhibitors exhibit significant differences in their interactions with the catalytic zinc metal ion and the other active site residues in the HDAC enzymes, which in turn affect the ‘Histidine-Aspartate’ charge relay mechanism. The QPLD calculations show that the protonated states of the inhibitors display higher scores in all the class I HDACs in this study, while the deprotonated forms present lower scores. The molecular electrostatic potentials and the other physico-chemical descriptors support the results. The MM-GBSA approach employed in the present work has been able to accurately calculate the relative binding free energies of the neutral and the protonated HA inhibitors; those were close to the experimental values. However, the MM-GBSA approach breaks down while calculating the binding free energies of the deprotonated inhibitors, which resulted in unrealistic values. Large energetic differences were found in the polar electrostatic solvation energy terms and the coulombic contributions in the deprotonated inhibitors. Thus improvements in the present solvation models and force fields become inevitable for the inclusions of charged states of inhibitors in computational drug discovery.

Toward quantitative estimates of binding affinities for protein–ligand systems involving large inhibitor compounds: A steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein-ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein-ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine-based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein-ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis-Menten mechanism with a competitive inhibitor.

Kinetics and Thermodynamics of Small Molecule Binding to Pincer-PCP Rhodium(I) Complexes

The kinetics and thermodynamics of the binding of several small molecules, L (L = N2, H2, D2, and C2H4), to the coordinatively unsaturated pincer-PCP rhodium(I) complexes Rh[(t)Bu2PCH2(C6H3)CH2P(t)Bu2] (1) and Rh[(t)Bu2P(CH2)2(CH)(CH2)2P(t)Bu2] (2) in organic solvents (n-heptane, toluene, THF, and cyclohexane-d12) have been investigated by a combination of kinetic flash photolysis methods, NMR equilibrium studies, and density functional theory (DFT) calculations. Using various gas mixtures and monitoring by NMR until equilibrium was established, the relative free energies of binding of N2, H2, and C2H4 in cyclohexane-d12 were found to increase in the order C2H4 < N2 < H2. Time-resolved infrared (TRIR) and UV-vis transient absorption spectroscopy revealed that 355 nm excitation of 1-L and 2-L results in the photoejection of ligand L. The subsequent mechanism of binding of L to 1 and 2 to regenerate 1-L and 2-L is determined by the structure of the PCP ligand framework and the nature of the solvent. In both cases, the primary transient is a long-lived, unsolvated species (τ = 50-800 ns, depending on L and its concentration in solution). For 2, this so-called less-reactive form (LRF) is in equilibrium with a more-reactive form (MRF), which reacts with L at diffusion-controlled rates to regenerate 2-L. These two intermediates are proposed to be different conformers of the three-coordinate (PCP)Rh fragment. For 1, a similar mechanism is proposed to occur, but the LRF to MRF step is irreversible. In addition, a parallel reaction pathway was observed that involves the direct reaction of the LRF of 1 with L, with second-order rate constants that vary by almost 3 orders of magnitude, depending on the nature of L (in n-heptane, k = 6.7 × 10(5) M(-1) s(-1) for L = C2H4; 4.0 × 10(6) M(-1) s(-1) for L = N2; 5.5 × 10(8) M(-1) s(-1) for L = H2). Experiments in the more coordinating solvent, THF, revealed the binding of THF to 1 to generate 1-THF, and its subsequent reaction with L, as a competing pathway.

Separated topologies—A method for relative binding free energy calculations using orientational restraints

Orientational restraints can improve the efficiency of alchemical free energy calculations, but they are not typically applied in relative binding calculations, which compute the affinity difference been two ligands. Here, we describe a new "separated topologies" method, which computes relative binding free energies using orientational restraints and which has several advantages over existing methods. While standard approaches maintain the initial and final ligand in a shared orientation, the separated topologies approach allows the initial and final ligands to have distinct orientations. This avoids a slowly converging reorientation step in the calculation. The separated topologies approach can also be applied to determine the relative free energies of multiple orientations of the same ligand. We illustrate the approach by calculating the relative binding free energies of two compounds to an engineered site in Cytochrome C Peroxidase.

Modifying the catalytic preference of tributyrin in Bacillus thermocatenulatus lipase through in-silico modeling of enzyme-substrate complex

In this study, rational design for Bacillus thermocatenulatus lipase (BTL2) was carried out to lower the activation barrier for hydrolysis of short-chain substrates. In this design, we used computational models for the enzyme-substrate (ES) complexes of tributyrin (C4) and tricaprylin (C8), which were generated through docking and molecular dynamics (MD) simulations. These ES complexes were employed in steered MD (SMD) simulations with Jarzynski's equality to estimate their relative binding free energies. Potential mutation sites for modifying the chain-length selectivity of BTL2 were found by inspecting the SMD trajectories and fine-tuning the volume and hydrophobicity of the cleft. Seven mutations (F17A, L57F, V175A, V175F, I320A, I320F and L360F) were performed to cover three binding pockets for sn-1, sn-2 and sn-3 acyl chains. The relative binding free energies of the mutant ES complexes formed by C4 and C8 ligands were calculated similarly. The experimental routines of protein engineering including site-directed mutagenesis, heterologous protein expression and purification were performed for all lipases. Steady-state specific activities towards C4 and C8 were determined for wild-type and mutant lipases, which gave an estimate of the relative change in the binding free energy of transition state complex (ES(‡)). The chain-length selectivity of mutants was determined from the relative changes in the activation barrier of hydrolysis of C4 and C8 triacylglycerol with respect to wild-type using computational and experimental findings. The most promising mutant for C4 over C8 preference was found to be L360F. We suggest that L360F may be at a critical position to lower the activation barrier for C4 and elevate it for C8 hydrolysis.

Study of base pair mutations in proline-rich homeodomain (PRH)–DNA complexes using molecular dynamics

Proline-rich homeodomain (PRH) is a regulatory protein controlling transcription and gene expression processes by binding to the specific sequence of DNA, especially to the sequence 5'-TAATNN-3'. The impact of base pair mutations on the binding between the PRH protein and DNA is investigated using molecular dynamics and free energy simulations to identify DNA sequences that form stable complexes with PRH. Three 20-ns molecular dynamics simulations (PRH-TAATTG, PRH-TAATTA and PRH-TAATGG complexes) in explicit solvent water were performed to investigate three complexes structurally. Structural analysis shows that the native TAATTG sequence forms a complex that is more stable than complexes with base pair mutations. It is also observed that upon mutation, the number and occupancy of the direct and water-mediated hydrogen bonds decrease. Free energy calculations performed with the thermodynamic integration method predict relative binding free energies of 0.64 and 2kcal/mol for GC to AT and TA to GC mutations, respectively, suggesting that among the three DNA sequences, the PRH-TAATTG complex is more stable than the two mutated complexes. In addition, it is demonstrated that the stability of the PRH-TAATTA complex is greater than that of the PRH-TAATGG complex.

Alchemical FEP Calculations of Ligand Conformer Focusing in Explicit Solvent.

Slow rotational degrees of freedom in ligands can make alchemical FEP simulations unreliable due to inadequate sampling. We addressed this problem by introducing a FEP-based protocol of ligand conformer focusing in explicit solvent. Our method involves FEP transformations between conformers using equilibrium dihedral angle as a reaction coordinate and provides the cost of "focusing" on one specific conformational state that binds to a protein. The calculated conformer focusing term made a considerable difference of 5-10 kJ/mol in computed relative binding free energies of studied Syk inhibitors and significantly improved the resulting accuracy of predictions.

Protein Kinase G Ligand Binding Specificity via Confine-And-Release Thermodynamic Integration

Protein Kinase G (PKG) is a second messenger dependent protein implicated in smooth muscle relaxation. Active PKG signals for blood vessel dilation, making it an attractive target in cardiovascular drug discovery. The PKG regulatory domain consists of two, potentially druggable, cyclic nucleotide binding domains (CNBDs), which activate PKG when simultaneously occupied by cGMP. Here we contribute to the characterization of each CNBD by examining its ability to discriminate between cGMP and a highly similar ligand, cAMP. Relative binding free energies were calculated using molecular dynamic simulations and thermodynamic integration (TI). The contribution of individual residues was assessed by repetition over a series of mutant structures. In some cases, multiple ligand binding conformations (syn vs. anti about the glycosyl bond) were found to meaningfully contribute to the binding free energies. As the barrier to rotation between these conformations prevented adequate sampling on computationally feasible TI timescales, the confine-and-release method was used to account for the contribution of each conformation. This resulted in improved convergence and agreement with experiment.

Reaction Control in Fungal Fatty Acid Synthase

Long-chain fatty acids form the membranes of every living cell and are an important component of energy storage and metabolism for higher organisms. Biosynthesis of fatty acids occurs over multiple iterations of a cycle of seven successively repeated reaction steps, each extending the growing acyl chain by two carbon atoms. Chain length is precisely controlled, with eukaryotic metabolic systems producing almost exclusively palmitic (C16) and stearic (C18) acids.Here we consider the fungal fatty acid synthase (FAS), a 2.6MDa multi-enzyme complex consisting of six heterodimers, which catalyses nine separate reactions across eight distinct types of active site. The high local concentrations resulting from this arrangement, further facilitated by an integral ‘shuttle’ domain, enable more efficient processing than could be achieved with freely diffusing enzymes.To elucidate how the extremely tight control of chain length is achieved from the perspective of the underlying molecular recognition processes, we used atomistic molecular modelling, molecular dynamics simulation and calculation of relative enzyme-substrate binding free energies. We focus on two key active sites, respectively responsible for elongation and termination, and establish that a relatively simple model of their interplay can explain the observed product spectrum, and that our method correctly predicts changes in chain length arising from tested mutations.These findings will direct ongoing work in which the synthetic machinery of the FAS complex is manipulated for the production of non-native products. The developed methods will contribute to techniques available for rational enzyme engineering.

Prediction of the Binding Affinities of PSD95 PDZ Domain in Complex with the CRIPT Peptide

In this work, we have applied the state of art molecular dynamics simulations in combination with solvation free energy and conformational entropy calculations to predict the binding affinities of PSD95 PDZ domain in complex with the CRIPT peptide. Four diffident computational protocols were evaluated on reproducing the relative binding free energies of the wild type PDZ and its five mutants. The protocol of MM-GB/SA in combination with normal mode analysis (NMA), which has a correlation coefficient square of 0.84, apparently outperforms the others especially for the two MM-PB/SA-based protocols. Free energy decomposition was also performed in order to identify the hot spots that contribute significantly to the binding. 

Comparison of the Molecular Dynamics and Calculated Binding Free Energies for Nine FDA‐Approved HIV‐1 PR Drugs Against Subtype B and C‐SA HIV PR

We report the first account of a comparative analysis of the binding affinities of nine FDA-approved drugs against subtype B as well as the South African subtype C HIV PR (C-SA). A standardized protocol was used to generate the inhibitor/C-SA PR complexes with the relative positions of the inhibitors taken from the corresponding X-ray structures for subtype B complexes. The dynamics and stability of these complexes were investigated using molecular dynamics calculations. Average relative binding free energies for these inhibitors were calculated from the molecular dynamics simulation using the molecular mechanics generalized Born surface area method. The calculated energies followed a similar trend to the reported experimental binding free energies. Postdynamic hydrogen bonding and electrostatic interaction analysis of the inhibitors with both subtypes reveal similar interactions. Most inhibitors show slightly weaker binding affinities for C-SA PR. Molecular dynamics studies demonstrated increased flap movement for C-SA PR, which can perhaps explain the weaker affinities. This study serves as a standardized platform for optimizing the design of future more potent HIV C-SA PR inhibitors.

An Integrated Computational Approach to Rationalize the Activity of Non-Zinc-Binding MMP-2 Inhibitors

Matrix metalloproteinases are a family of Zn-proteases involved in tissue remodeling and in many pathological conditions. Among them MMP-2 is one of the most relevant target in anticancer therapy. Commonly, MMP inhibitors contain a functional group able to bind the zinc ion and responsible for undesired side effects. The discovery of potent and selective MMP inhibitors not bearing a zinc-binding group is arising for some MMP family members and represents a new opportunity to find selective and non toxic inhibitors.In this work we attempted to get more insight on the inhibition process of MMP-2 by two non-zinc-binding inhibitors, applying a general protocol that combines several computational tools (docking, Molecular Dynamics and Quantum Chemical calculations), that all together contribute to rationalize experimental inhibition data. Molecular Dynamics studies showed both structural and mechanical-dynamical effects produced by the ligands not disclosed by docking analysis. Thermodynamic Integration provided relative binding free energies consistent with experimentally observed activity data. Quantum Chemical calculations of the tautomeric equilibrium involving the most active ligand completed the picture of the binding process. Our study highlights the crucial role of the specificity loop and suggests that enthalpic effect predominates over the entropic one.

Molecular Basis of Drug Resistance in A/H1N1 Virus

New mutants of human influenza virus (A/H1N1) exhibit resistance to antiviral drugs. The mechanism whereby they develop insensitivity to these medications is, however, not yet completely understood. A crystallographic structure of A/H1N1 neuraminidase has been published recently. Using molecular dynamic simulations, it is now possible to characterize at the atomic level the mechanism that underlies the loss of binding affinity of the drugs. In this study, free-energy perturbation was used to evaluate the relative binding free energies of Tamiflu and Relenza with H274Y, N294S, and Y252H neuraminidase mutants. Our results demonstrate a remarkable correlation between theoretical and experimental data, which quantitatively confirms that the mutants are resistant to Tamiflu but are still strongly inhibited by Relenza. The simulations further reveal the key interactions that govern the affinity of the two drugs for each mutant. This information is envisioned to prove useful for the design of novel neuraminidase inhibitors and for the characterization of new potential mutants.

Adsorption Behavior of Hydrophobin Proteins on Polydimethylsiloxane Substrates

The design of a bioactive surface with appropriate wettability for effective protein immobilization has attracted much attention. Previous experiments showed that the adsorption of hydrophobic protein HFBI onto a polydimethylsiloxane (PDMS) substrate surface can reverse the inherent hydrophobicity of the surface, hence making it suitable for immobilization of a secondary protein. In this study, atomistic molecular dynamics simulations have been conducted to elucidate the adsorption mechanism of HFBI on the PDMS substrate in an aqueous environment. Nine independent simulations starting from three representative initial orientations of HFBI toward the solid surface were performed, resulting in different adsorption modes. The main secondary structures of the protein in each mode are found to be preserved in the entire course of adsorption due to the four disulfide bonds. The relative binding free energies of the different adsorption modes were calculated, showing that the mode, in which the binding residues of HFBI fully come from its hydrophobic patch, is most energetically favored. In this favorable binding mode, the hydrophilic region of HFBI is fully exposed to water, leading to a high hydrophilicity of the modified PDMS surface, consistent with experiments. Furthermore, a set of residues consisting of Leu12, Leu24, Leu26, Ile27, Ala66, and Leu68 were found to play an important role in the adsorption of HFBI on different hydrophobic substrates, irrespective of the structural features of the substrates.

Computation of relative binding free energy for an inhibitor and its analogs binding with Erk kinase using thermodynamic integration MD simulation

In the present study, we carried out thermodynamic integration molecular dynamics simulation for a pair of analogous inhibitors binding with Erk kinase to investigate how computation performs in reproducing the relative binding free energy. The computation with BCC-AM1 charges for ligands gave -1.1kcal/mol, deviated from experimental value of -2.3kcal/mol by 1.2kcal/mol, in good agreement with experimental result. The error of computed value was estimated to be 0.5kcal/mol. To obtain convergence, switching vdw interaction on and off required approximately 10 times more CPU time than switching charges. Residue-based contributions and hydrogen bonding were analyzed and discussed. Furthermore, subsequent simulation using RESP charge for ligand gave ΔΔG of -1.6kcal/mol. The computed results are better than the result of -5.6kcal/mol estimated using PBSA method in a previous study. Based on these results, we further carried out computations to predict ΔΔG for five new analogs, focusing on placing polar and nonpolar functional groups at the meta site of benzene ring shown in the Fig.1, to see if these ligands have better binding affinity than the above ligands. The computations resulted that a ligand with polar -OH group has better binding affinity than the previous examined ligand by ~2.0kcal/mol and two other ligands have better affinity by ~1.0kcal/mol. The predicted better inhibitors of this kind should be of interest to experimentalist for future experimental enzyme and/or cell assays.

Configurational Preferences of Arylamide α-Helix Mimetics via Alchemical Free Energy Calculations of Relative Binding Affinities

We use molecular docking and free energy calculations to estimate the relative free energy of binding of six arylamide compounds designed to inhibit the hDM2-p53 interaction. We show that using docking methods to predict or rank the binding affinity of a series of arylamide inhibitors of the hDM2-p53 interaction is problematic. However, using free energy calculations, we show that we can achieve levels of accuracy that can guide the development of novel arylamide compounds. We perform alchemical free energy calculations using the Desmond molecular dynamics package with the same arylamide inhibitors of the hDM2-p53 system and illustrate the challenges of performing accurate free energy calculations for realistic systems. To our knowledge, these are the first calculations for inhibitors of the hDM2 system that employ a full treatment of statistical mechanics including explicit water representation and full protein flexibility. We show that mutating three functional groups in a single transformation can be more efficient than mutating the groups one by one if proper intermediates are used. We also show that Hamiltonian exchanges can improve the efficiency of the calculation compared to standard alchemical methods, with a novel use of the phase space overlap to monitor sampling extent. We show that, despite sampling limitations, this approach can achieve levels of accuracy sufficient to bias further inhibitor modification toward binding, and identifies antiparallel configurations as stable or more stable than the parallel configurations that are typically considered.