Using the fast fourier transform in binding free energy calculations

According to implicit ligand theory, the standard binding free energy is an exponential average of the binding potential of mean force (BPMF), an exponential average of the interaction energy between the unbound ligand ensemble and a rigid receptor. Here, we use the fast Fourier transform (FFT) to efficiently evaluate BPMFs by calculating interaction energies when rigid ligand configurations from the unbound ensemble are discretely translated across rigid receptor conformations. Results for standard binding free energies between T4 lysozyme and 141 small organic molecules are in good agreement with previous alchemical calculations based on (1) a flexible complex ( for 24 systems) and (2) flexible ligand with multiple rigid receptor configurations ( for 141 systems). While the FFT is routinely used for molecular docking, to our knowledge this is the first time that the algorithm has been used for rigorous binding free energy calculations. © 2017 Wiley Periodicals, Inc.     

Two‐dimensional replica‐exchange method for predicting protein–ligand binding structures

We have developed a two‐dimensional replica‐exchange method for the prediction of protein–ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein‐ligand complex systems (MDM2 and HSP 90‐alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica‐exchange umbrella sampling method that we have previously developed. © 2013 Wiley Periodicals, Inc.     

Calculation of binding free energies of inhibitors to plasmepsin II

An understanding at the atomic level of the driving forces of inhibitor binding to the protein plasmepsin (PM) II would be of interest to the development of drugs against malaria. To this end, three state of the art computational techniques to compute relative free energies-thermodynamic integration (TI), Hamiltonian replica-exchange (H-RE) TI, and comparison of bound versus unbound ligand energy and entropy-were applied to a protein-ligand system of PM II and several exo-3-amino-7-azabicyclo[2.2.1]heptanes and the resulting relative free energies were compared with values derived from experimental IC(50) values. For this large and flexible protein-ligand system, the simulations could not properly sample the relevant parts of the conformational space of the bound ligand, resulting in failure to reproduce the experimental data. Yet, the use of Hamiltonian replica exchange in conjunction with thermodynamic integration resulted in enhanced convergence and computational efficiency compared to standard thermodynamic integration calculations. The more approximate method of calculating only energetic and entropic contributions of the ligand in its bound and unbound states from conventional molecular dynamics (MD) simulations reproduced the major trends in the experimental binding free energies, which could be rationalized in terms of energetic and entropic characteristics of the different structural and physico-chemical properties of the protein and ligands.     

Improvement of Ligand Affinity and Thermodynamic Properties by NMR‐Based Evaluation of Local Dynamics and Surface Complementarity in the Receptor‐Bound State

The thermodynamic properties of a ligand in the bound state affect its binding specificity. Strict binding specificity can be achieved by introducing multiple spatially defined interactions, such as hydrogen bonds and van der Waals interactions, into the ligand–receptor interface. These introduced interactions are characterized by restricted local dynamics and improved surface complementarity in the bound state. In this study, we experimentally evaluated the local dynamics and the surface complementarity of weak‐affinity ligands in the receptor‐bound state by forbidden coherence transfer analysis in free‐bound exchange systems (Ex‐FCT), using the interaction between a ligand, a myocyte‐enhancer factor 2A (MEF2A) docking peptide, and a receptor, p38α, as a model system. The Ex‐FCT analyses successfully provided information for the rational design of a ligand with higher affinity and preferable thermodynamic properties for p38α.     

Polymer-tethered ligand–receptor interactions between surfaces

We study the interaction between two parallel surfaces having a polymer-tethered ligand on one and a random distribution of receptors on the other. We examine the interplay between the specific ligand–receptor binding and the conformation degrees of freedom of the polymer tether, and address the difference between the cases of mobile (annealed) and immobile (quenched) receptors. The annealed case is solved exactly and the quenched case is treated by both Monte Carlo direct sampling and an analytical density expansion. The combination of the entropic repulsion due to chain confinement at small surface separations, and the attraction due to ligand–receptor binding which can take place at significant chain stretching, results in a minimum in the interaction free energy. For the same set of parameters, stronger binding is obtained for the annealed case than for the quenched case, reflecting the ability of the mobile receptors to migrate into the region of the ligand. In the quenched case, binding is limited by the availability of receptors within the reach of the ligand; for a given receptor density and binding energy, there exists an optimal chain length that yields the lowest minimum in the interaction free energy for the quenched case. © 2006 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 44: 2621–2637, 2006     

Classification of water molecules in protein binding sites

Water molecules play a crucial role in mediating the interaction between a ligand and a macromolecular receptor. An understanding of the nature and role of each water molecule in the active site of a protein could greatly increase the efficiency of rational drug design approaches: if the propensity of a water molecule for displacement can be determined, then synthetic effort may be most profitably applied to the design of specific ligands with the displacement of this water molecule in mind. In this paper, a thermodynamic analysis of water molecules in the binding sites of six proteins, each complexed with a number of inhibitors, is presented. Two classes of water molecules were identified: those conserved and not displaced by any of the ligands, and those that are displaced by some ligands. The absolute binding free energies of 54 water molecules were calculated using the double decoupling method, with replica exchange thermodynamic integration in Monte Carlo simulations. It was found that conserved water molecules are on average more tightly bound than displaced water molecules. In addition, Bayesian statistics is used to calculate the probability that a particular water molecule may be displaced by an appropriately designed ligand, given the calculated binding free energy of the water molecule. This approach therefore allows the numerical assessment of whether or not a given water molecule should be targeted for displacement as part of a rational drug design strategy.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Exploring Multi-Subsite Binding Pockets in Proteins: DEEP-STD NMR Fingerprinting and Molecular Dynamics Unveil a Cryptic Subsite at the GM1 Binding Pocket of Cholera Toxin B

Ligand-based NMR techniques to study protein–ligand interactions are potent tools in drug design. Saturation transfer difference (STD) NMR spectroscopy stands out as one of the most versatile techniques, allowing screening of fragments libraries and providing structural information on binding modes. Recently, it has been shown that a multi-frequency STD NMR approach, differential epitope mapping (DEEP)-STD NMR, can provide additional information on the orientation of small ligands within the binding pocket. Here, the approach is extended to a so-called DEEP-STD NMR fingerprinting technique to explore the binding subsites of cholera toxin subunit B (CTB). To that aim, the synthesis of a set of new ligands is presented, which have been subject to a thorough study of their interactions with CTB by weak affinity chromatography (WAC) and NMR spectroscopy. Remarkably, the combination of DEEP-STD NMR fingerprinting and Hamiltonian replica exchange molecular dynamics has proved to be an excellent approach to explore the geometry, flexibility, and ligand occupancy of multi-subsite binding pockets. In the particular case of CTB, it allowed the existence of a hitherto unknown binding subsite adjacent to the GM1 binding pocket to be revealed, paving the way to the design of novel leads for inhibition of this relevant toxin.     

Protein Allostery at Atomic Resolution

Protein allostery is a phenomenon involving the long range coupling between two distal sites in a protein. In order to elucidate allostery at atomic resoluion on the ligand-binding WW domain of the enzyme Pin1, multistate structures were calculated from exact nuclear Overhauser effect (eNOE). In its free form, the protein undergoes a microsecond exchange between two states, one of which is predisposed to interact with its parent catalytic domain. In presence of the positive allosteric ligand, the equilibrium between the two states is shifted towards domain–domain interaction, suggesting a population shift model. In contrast, the allostery-suppressing ligand decouples the side-chain arrangement at the inter-domain interface thereby reducing the inter-domain interaction. As such, this mechanism is an example of dynamic allostery. The presented distinct modes of action highlight the power of the interplay between dynamics and function in the biological activity of proteins.     

Fast,metadynamics‐based method for prediction of the stereochemistry‐dependent relative free energies of ligand–receptor interactions

The computational approach applicable for the molecular dynamics (MD)‐based techniques is proposed to predict the ligand–protein binding affinities dependent on the ligand stereochemistry. All possible stereoconfigurations are expressed in terms of one set of force‐field parameters [stereoconfiguration‐independent potential (SIP)], which allows for calculating all relative free energies by only single simulation. SIP can be used for studying diverse, stereoconfiguration‐dependent phenomena by means of various computational techniques of enhanced sampling. The method has been successfully tested on the β₂‐adrenergic receptor (β₂‐AR) binding the four fenoterol stereoisomers by both metadynamics simulations and replica‐exchange MD. Both the methods gave very similar results, fully confirming the presence of stereoselective effects in the fenoterol‐β₂‐AR interactions. However, the metadynamics‐based approach offered much better efficiency of sampling which allows for significant reduction of the unphysical region in SIP. © 2014 Wiley Periodicals, Inc.     

Identifying ligand binding sites and poses using GPU-accelerated Hamiltonian replica exchange molecular dynamics

We present a method to identify small molecule ligand binding sites and poses within a given protein crystal structure using GPU-accelerated Hamiltonian replica exchange molecular dynamics simulations. The Hamiltonians used vary from the physical end state of protein interacting with the ligand to an unphysical end state where the ligand does not interact with the protein. As replicas explore the space of Hamiltonians interpolating between these states, the ligand can rapidly escape local minima and explore potential binding sites. Geometric restraints keep the ligands from leaving the vicinity of the protein and an alchemical pathway designed to increase phase space overlap between intermediates ensures good mixing. Because of the rigorous statistical mechanical nature of the Hamiltonian exchange framework, we can also extract binding free energy estimates for all putative binding sites. We present results of this methodology applied to the T4 lysozyme L99A model system for three known ligands and one non-binder as a control, using an implicit solvent. We find that our methodology identifies known crystallographic binding sites consistently and accurately for the small number of ligands considered here and gives free energies consistent with experiment. We are also able to analyze the contribution of individual binding sites to the overall binding affinity. Our methodology points to near term potential applications in early-stage structure-guided drug discovery.     

A Monte Carlo method for finding important ligand fragments from receptor data

A simulated annealing method for finding important ligand fragments is described. At a given temperature, ligand fragments are randomly selected and randomly placed within the given receptor cavity, often replacing or forming bonds with existing ligand fragments. For each new ligand fragment combination, the bonded, nonbonded, polarization and solvation energies of the new ligand–receptor system are compared to the previous configuration. Acceptance or rejection of the new system is decided using the Boltzmann distribution , where E is the energy difference between the old and new systems, k is the Boltzmann constant and T is the temperature. Thus, energetically unfavorable fragment switches are sometimes accepted, sacrificing immediate energy gains in the interest of finding a system with minimum energy. By lowering the temperature, the rate of unfavorable switches decreases and energetically favorable combinations become more difficult to change. The process is terminated when the frequency of switches becomes too small. As a test, the method predicted positions and types of important ligand fragments for neuraminidase that were in accord with the known ligand, sialic acid.     

Gaussian mapping of chemical fragments in ligand binding sites

We present a new approach to automatically define a quasi-optimal minimal set of pharmacophoric points mapping the interaction properties of a user-defined ligand binding site. The method is based on a fitting algorithm where a grid of sampled interaction energies of the target protein with small chemical fragments in the binding site is approximated by a linear expansion of Gaussian functions. A heuristic approximation selects from this expansion the smallest possible set of Gaussians required to describe the interaction properties of the binding site within a prespecified accuracy. We have evaluated the performance of the approach by comparing the computed Gaussians with the positions of aromatic sites found in experimental protein-ligand complexes. For a set of 53 complexes, good correspondence is found in general. At a 95% significance level, approximately 65% of the predicted interaction points have an aromatic binding site within 1.5 A. We then studied the utility of these points in docking using the program DOCK. Short docking times, with an average of approximately 0.18 s per conformer, are obtained, while retaining, both for rigid and flexible docking, the ability to sample native-like binding modes for the ligand. An average 4-5-fold speed-up in docking times and a similar success rate is estimated with respect to the standard DOCK protocol.     

Modified replica exchange simulation methods for local structure refinement

Parallel tempering, also known as replica exchange molecular dynamics (REMD), has recently been successfully used to study the structure and thermodynamic properties of biomolecules such as peptides and small proteins. For large systems, however, applying REMD can be costly since the number of replicas needed increases as the square root of the number of degrees of freedom in the system. Often, enhanced sampling is only needed for a subset of atoms, such as a loop region of a large protein or a small ligand binding to a receptor. In such applications, it is often reasonable to assume a weak dependence of the structure of the larger region on the instantaneous conformation of the smaller region of interest. For these cases, we derived two variant replica exchange methods, partial replica exchange molecular dynamics (PREMD) and local replica exchange molecular dynamics (LREMD). The Hamiltonian for the system is separated, with replica exchange carried out only for terms involving the subsystem of interest while the remainder of the system is maintained at a single temperature. The number of replicas required for efficient exchange thus depends on the number of degrees of freedom in the fragment needing refinement rather than on the size of the full system. The method can be applied to much larger systems than was previously practical. This also provides a means to preserve the integrity of the structure outside the refinement region without introduction of restraints. LREMD takes this weak coupling approximation a step further, employing only a single representation of the large fragment that simultaneously interacts with all of the replicas of the subsystem of interest. This is obtained by combining replica exchange with the locally enhanced sampling approximation (LES), reducing the computational expense of replica exchange simulations to near that of a single standard molecular dynamics (MD) simulation. Use of LREMD also permits the use of LES without requiring the specification of a single temperature, a known difficulty for standard LES simulations. We tested these two methods on the loop region of an RNA hairpin model system and find significant advantages over standard MD and REMD simulations.     

Multiple-site ligand binding to flexible macromolecules: Separation of global and local conformational change and an iterative mobile clustering approach

This article concerns the calculation of equilibria of ligand binding to multiple sites in macromolecules in the presence of conformational flexibility and conformation-dependent interaction among the sites. A formulation of this problem is presented in which global conformational changes are distinguished from conformational changes that are confined to “locally flexible regions.” The formalism is quite general in that ligands of different types, multivalent binding sites, tautomeric binding sites, and sites that bind more than one type of ligand can be accommodated. Strictly speaking, the separation of the conformational problem into global and local parts does not impose any loss of generality, although in practice it is necessary to restrict the number of global and local conformers. Because of the combinatorics of binding and conformational states, the computational complexity of a problem having only local conformational flexibility grows exponentially with the number of sites and the number of locally flexible regions. An iterative mobile clustering method for cutting off this exponential growth and obtaining approximate solutions with low computational cost is presented and tested. In this method, a binding site is selected, and a “cluster” of strongly interacting sites is set up around it; within the cluster, the binding and conformational states are fully enumerated, whereas the influences of sites outside the cluster on the sites inside are treated by a mean field approximation. The procedure then moves to the next site around which another (possibly overlapping) cluster is formed and the calculation is repeated. The procedure iterates through the list of sites in this way, using the results of previous iterations for the mean-field terms of current iterations until a convergence criterion is met. The method is tested on a large set of randomly generated problems of varying size, whose geometries are chosen to have protein-like statistical properties. It is found that the method is accurate and rapid with the computational cost scaling linearly to quadratically with the number of sites, except for a minority of cases in which large clusters occur by chance. The new method is more accurate than a Monte Carlo method, and may be faster or slower depending on the clustering criteria and details of the macromolecule. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1091–1111, 1999     

MCDOCK: A Monte Carlo simulation approach to the molecular docking problem

Prediction of the binding mode of a ligand (a drug molecule) to its macromolecular receptor, or molecular docking, is an important problem in rational drug design. We have developed a new docking method in which a non-conventional Monte Carlo (MC) simulation technique is employed. A computer program, MCDOCK, was developed to carry out the molecular docking operation automatically. The current version of the MCDOCK program (version 1.0) allows for the full flexibility of ligands in the docking calculations. The scoring function used in MCDOCK is the sum of the interaction energy between the ligand and its receptor, and the conformational energy of the ligand. To validate the MCDOCK method, 19 small ligands, the binding modes of which had been determined experimentally using X-ray diffraction, were docked into their receptor binding sites. To produce statistically significant results, 20 MCDOCK runs were performed for each protein–ligand complex. It was found that a significant percentage of these MCDOCK runs converge to the experimentally observed binding mode. The root-mean-square (rms) of all non-hydrogen atoms of the ligand between the predicted and experimental binding modes ranges from 0.25 to 1.84 Å for these 19 cases. The computational time for each run on an SGI Indigo2/R10000 varies from less than 1 min to 15 min, depending upon the size and the flexibility of the ligands. Thus MCDOCK may be used to predict the precise binding mode of ligands in lead optimization and to discover novel lead compounds through structure-based database searching.     

π-Hole Interactions Involving Nitro Aromatic Ligands in Protein Structures

Studying noncanonical intermolecular interactions between a ligand and a protein constitutes an emerging research field. Identifying synthetically accessible molecular fragments that can engage in intermolecular interactions is a key objective in this area. Here, it is shown that so-called “π-hole interactions” are present between the nitro moiety in nitro aromatic ligands and lone pairs within protein structures (water and protein carbonyls and sulfurs). Ample structural evidence was found in a PDB analysis and computations reveal interaction energies of about −5 kcal mol⁻¹ for ligand–protein π-hole interactions. Several examples are highlighted for which a π-hole interaction is implicated in the superior binding affinity or inhibition of a nitro aromatic ligand versus a similar non-nitro analogue. The discovery that π-hole interactions with nitro aromatics are significant within protein structures parallels the finding that halogen bonds are biologically relevant. This has implications for the interpretation of ligand–protein complexation phenomena, for example, involving the more than 50 approved drugs that contain a nitro aromatic moiety.     

Q-Dock: Low-resolution flexible ligand docking with pocket-specific threading restraints

The rapidly growing number of theoretically predicted protein structures requires robust methods that can utilize low-quality receptor structures as targets for ligand docking. Typically, docking accuracy falls off dramatically when apo or modeled receptors are used in docking experiments. Low-resolution ligand docking techniques have been developed to deal with structural inaccuracies in predicted receptor models. In this spirit, we describe the development and optimization of a knowledge-based potential implemented in Q-Dock, a low-resolution flexible ligand docking approach. Self-docking experiments using crystal structures reveals satisfactory accuracy, comparable with all-atom docking. All-atom models reconstructed from Q-Dock's low-resolution models can be further refined by even a simple all-atom energy minimization. In decoy-docking against distorted receptor models with a root-mean-square deviation, RMSD, from native of approximately 3 A, Q-Dock recovers on average 15-20% more specific contacts and 25-35% more binding residues than all-atom methods. To further improve docking accuracy against low-quality protein models, we propose a pocket-specific protein-ligand interaction potential derived from weakly homologous threading holo-templates. The success rate of Q-Dock employing a pocket-specific potential is 6.3 times higher than that previously reported for the Dolores method, another low-resolution docking approach.