Estimating binding affinities by docking/scoring methods using variable protonation states

To investigate the effects of multiple protonation states on protein-ligand recognition, we generated alternative protonation states for selected titratable groups of ligands and receptors. The selection of states was based on the predicted pK(a) of the unbound receptor and ligand and the proximity of titratable groups of the receptor to the binding site. Various ligand tautomer states were also considered. An independent docking calculation was run for each state. Several protocols were examined: using an ensemble of all generated states of ligand and receptor, using only the most probable state of the unbound ligand/receptor, and using only the state giving the most favorable docking score. The accuracies of these approaches were compared, using a set of 176 protein-ligand complexes (15 receptors) for which crystal structures and measured binding affinities are available. The best agreement with experiment was obtained when ligand poses from experimental crystal structures were used. For 9 of 15 receptors, using an ensemble of all generated protonation states of the ligand and receptor gave the best correlation between calculated and measured affinities.     

Challenges and advances in computational docking: 2009 in review

Docking is a computational technique that places a small molecule (ligand) in the binding site of its macromolecular target (receptor) and estimates its binding affinity. This review addresses methodological developments that have occurred in the docking field in 2009, with a particular focus on the more difficult, and sometimes controversial, aspects of this promising computational discipline. These developments aim to address the main challenges of docking: receptor representation (such aspects as structural waters, side chain protonation, and, most of all, flexibility (from side chain rotation to domain movement)), ligand representation (protonation, tautomerism and stereoisomerism, and the effect of input conformation), as well as accounting for solvation and entropy of binding. This review is strongly focused on docking advances in the context of drug design, specifically in virtual screening and fragment-based drug design.     

The molecular basis for ligand specificity in a mouse olfactory receptor: a network of functionally important residues

Sequence differences between members of the mouse olfac-tory receptor MOR42 subfamily (MOR42-3 and MOR42-1) are likely to be the basis for variation in ligand binding preference among these receptors. We investigated the specificity of MOR42-3 for a variety of dicarboxylic acids. We used site-directed mutagenesis, guided by homology modeling and ligand docking studies, to locate functionally important residues. Receptors were expressed in Xenopus oocytes and assayed using high throughput electrophysiology. The importance of the Val-113 residue, located deep within the receptor, was analyzed in the context of interhelical interactions. We also screened additional residues predicted to be involved in ligand binding site, based on comparison of ortholog/paralog pairs from the mouse and human olfactory receptor genomes (Man, O., Gilad, Y., and Lancet, D. (2004) Protein Sci. 13, 240-254). A network of 8 residues in transmembrane domains III, V, and VI was identified. These residues form part of the ligand binding pocket of MOR42-3. C12 dicarboxylic acid did not activate the receptor in our functional assay, yet our docking simulations predicted its binding site in MOR42-3. Binding without activation implied that C12 dicarboxylic acid might act as an antagonist. In our functional assay, C12 dicarboxylic acid did indeed act as an antagonist of MOR42-3, in agreement with molecular docking studies. Our results demonstrate a powerful approach based on the synergy between computational predictions and physiological assays.     

Structural and Energetic Effects of A2A Adenosine Receptor Mutations on Agonist and Antagonist Binding

To predict structural and energetic effects of point mutations on ligand binding is of considerable interest in biochemistry and pharmacology. This is not only useful in connection with site-directed mutagenesis experiments, but could also allow interpretation and prediction of individual responses to drug treatment. For G-protein coupled receptors systematic mutagenesis has provided the major part of functional data as structural information until recently has been very limited. For the pharmacologically important A_2A adenosine receptor, extensive site-directed mutagenesis data on agonist and antagonist binding is available and crystal structures of both types of complexes have been determined. Here, we employ a computational strategy, based on molecular dynamics free energy simulations, to rationalize and interpret available alanine-scanning experiments for both agonist and antagonist binding to this receptor. These computer simulations show excellent agreement with the experimental data and, most importantly, reveal the molecular details behind the observed effects which are often not immediately evident from the crystal structures. The work further provides a distinct validation of the computational strategy used to assess effects of point-mutations on ligand binding. It also highlights the importance of considering not only protein-ligand interactions but also those mediated by solvent water molecules, in ligand design projects.     

Long-range coupling between separate docking sites in interleukin-1beta

The human cytokine interleukin-1beta (IL-1beta) interacts with the interleukin type I receptor using two large docking surfaces designated A and B. Crystallographic studies reveal that a single histidine residue (His30) in IL-1beta makes critical electrostatic interactions at the receptor/ligand interface. To study the function of this residue at site A, four mutant forms of IL-1beta (H30A, H30D, H30F and H30R) were investigated. The mutation that introduces charge repulsion at His30 destabilizes the protein, but paradoxically causes the least effect on receptor binding (H30D). Mutations that enhance hydrophobic or electrostatic interactions have little effect on protein stability yet markedly impair receptor binding (H30F, H30R). All mutations can transmit effects from site A to site B, as evidenced by changes in the binding of a single-chain antibody highly specific for site B. Dihedral scalar coupling constants for the wild-type IL-1beta and the four His mutant proteins showed changes in backbone angles in residues located around site B, some approximately 30 angstroms away from His30 in site A. A comparison of native solvent exchange in wild-type and mutated IL-1beta shows transmission of local destabilization along the hydrogen bond network of the beta-sheet. Taken together, the data indicate that a single residue in site A of IL-1beta can impact stability and function through perturbations in both local and long-range contacts.     

Steric interactions and the activity of fentanyl analogs at the mu-opioid receptor

Fentanyl is a highly potent and clinically widely used narcotic analgesic. The synthesis of its analogs remains a challenge in the attempt to develop highly selective mu-opioid receptor agonists with specific pharmacological properties. In this paper, the use of flexible molecular docking in a study of the formation of complexes between a series of active fentanyl analogs and the mu-opioid receptor is described. The optimal position and orientation of fourteen fentanyl analogs in the binding pocket of the mu-receptor were determined. The major receptor amino acids and the ligand functional groups participating in the complex formation were identified. Stereochemical effects on the potency and binding are explained. The proposed model of ligand-receptor binding is in agreement with point mutation experiments explaining the role of the amino acids: Asp147, Tyr148, Asn230, His297, Trp318, His319, Cys321, and Tyr326 in the complex formation. In addition, the following amino acids were identified as being important for ligand binding or receptor activation: Ile322, Gly325, Val300, Met203, Leu200, Val143, and Ile144.     

Flexible docking using tabu search and an empirical estimate of binding affinity

This article describes the implementation of a new docking approach. The method uses a Tabu search methodology to dock flexibly ligand molecules into rigid receptor structures. It uses an empirical objective function with a small number of physically based terms derived from fitting experimental binding affinities for crystallographic complexes. This means that docking energies produced by the searching algorithm provide direct estimates of the binding affinities of the ligands. The method has been tested on 50 ligand-receptor complexes for which the experimental binding affinity and binding geometry are known. All water molecules are removed from the structures and ligand molecules are minimized in vacuo before docking. The lowest energy geometry produced by the docking protocol is within 1.5 Å root-mean square of the experimental binding mode for 86% of the complexes. The lowest energies produced by the docking are in fair agreement with the known free energies of binding for the ligands. Proteins 33:367–382, 1998. © 1998 Wiley-Liss, Inc.     

Site-specific basicities regulate molecular recognition in receptor binding: in silico docking of thyroid hormones

Interactions between thyroid hormone α and β receptors and the eight protonation microspecies of each of the main thyroid hormones (thyroxine, liothyronine, and reverse liothyronine) were investigated and quantitated by molecular modeling. Flexible docking of the various protonation forms of thyroid hormones and high-affinity thyromimetics to the two thyroid receptors was carried out. In this method the role of the ionization state of each basic site could be studied in the composite process of molecular recognition. Our results quantitate at the molecular level how the ionization state and the charge distribution influence the protein binding. The anionic form of the carboxyl group (i.e., carboxylate site) is essential for protein binding, whereas the protonated form of amino group worsens the binding. The protonation state of the phenolate plays a less important role in the receptor affinity; its protonation, however, alters the electron density and the concomitant stacking propensity of the aromatic rings, resulting in a different binding score. The combined results of docking and microspeciation studies show that microspecies with the highest concentration at the pH of blood are not the strongest binding ones. The calculated binding free energy values can be well interpreted in terms of the interactions between the actual sites of the microspecies and the receptor amino acids. Our docking results were validated and compared with biological data from the literature. Since the thyroid hormone receptors influence several physiologic functions, such as metabolic rate, cholesterol and triglyceride levels, and heart frequency, our binding results provide a molecular basis for drug design and development in related therapeutic indications.     

The pH-dependent structural variation of complementarity-determining region H3 in the crystal structures of the Fv fragment from an anti-dansyl monoclonal antibody.

The Fv fragment from an anti-dansyl antibody was optimally crystallized into two crystal forms having slightly different lattice dimensions at pH 5.25 and 6.75. The two crystal structures were determined and refined at high resolution at 112 K (at 1.45 A for the crystal at pH 5.25 and at 1.55 A for that at pH 6.75). In the two crystal structures, marked differences were identified in the first half of CDRH3 s having an amino acid sequence of Ile95H-Tyr96H-Tyr97H-His98H-Tyr99H-Pro1 00H-Trp100aH-Phe100bH-Ala101H- Tyr102H. NMR pH titration experiments revealed the p Kavalues of four histidine residues (His27dL, His93L, His55H and His98H) exposed to solvent. Only His98H (p Ka=6.3) completely changed its protonation state between the two crystallization conditions. In addition, the environmental structures including hydration water molecules around the four histidine residues were carefully compared. While the hydration structures around His27dL, His93L and His55H were almost invariant between the two crystal structures, those around His98Hs showed great difference in spite of the small conformational difference of His98H between the two crystal structures. These spectroscopic and crystallographic findings suggested that the change in the protonation state in His98H was responsible for the structural differences between pH 5.25 and 6.75. In addition, the most plausible binding site of the dansyl group was mapped into the present structural models with our previous NMR experimental results. The complementarity-determining regions H1, H3 and the N-terminal region in the VH domain formed the site. The side-chain of Tyr96H occupied the site and interacted with Phe27H of H1, giving a clue for the binding mode of the dansyl group in the site.     

Molten globule formation in apomyoglobin monitored by the fluorescent probe Nile Red

The interaction of nile red (NR) with apomyoglobin (ApoMb) in the native (pH 7) and molten globule (pH 4) states was investigated using experimental and computational methods. NR binds to hydrophobic locations in ApoMb with higher affinity (K(d) = 25 +/- 5 microM) in the native state than in the molten globule state (K(d) = 52 +/- 5 microM). In the molten globule state, NR is located in a more hydrophobic environment. The dye does not bind to the holoprotein, suggesting that the binding site is located at the heme pocket. In addition to monitoring steady-state properties, the fluorescence emission of NR is capable of tracking submillisecond, time-resolved structural rearrangements of the protein, induced by a nanosecond pH jump. Molecular dynamics simulations were run on ApoMb at neutral pH and at pH 4. The structure obtained for the molten globule state is consistent with the experimentally available structural data. The docking of NR with the crystal structure shows that the ligand binds into the binding pocket of the heme group, with an orientation bringing the planar ring system of NR to overlap with the position of two of the heme porphyrin rings in Mb. The docking of NR with the ApoMb structure at pH 4 shows that the dye binds to the heme pocket with a slightly less favorable binding energy, in keeping with the experimental K(d) value. Under these conditions, NR is positioned in a different orientation, reaching a more hydrophobic environment in agreement with the spectroscopic data.     

Efficient docking of peptides to proteins without prior knowledge of the binding site

Reliability in docking of ligand molecules to proteins or other targets is an important challenge for molecular modeling. Applications of the docking technique include not only prediction of the binding mode of novel drugs, but also other problems like the study of protein-protein interactions. Here we present a study on the reliability of the results obtained with the popular AutoDock program. We have performed systematical studies to test the ability of AutoDock to reproduce eight different protein/ligand complexes for which the structure was known, without prior knowledge of the binding site. More specifically, we look at factors influencing the accuracy of the final structure, such as the number of torsional degrees of freedom in the ligand. We conclude that the Autodock program package is able to select the correct complexes based on the energy without prior knowledge of the binding site. We named this application blind docking, as the docking algorithm is not able to "see" the binding site but can still find it. The success of blind docking represents an important finding in the era of structural genomics.     

Receptor–ligand molecular docking

Docking methodology aims to predict the experimental binding modes and affinities of small molecules within the binding site of particular receptor targets and is currently used as a standard computational tool in drug design for lead compound optimisation and in virtual screening studies to find novel biologically active molecules. The basic tools of a docking methodology include a search algorithm and an energy scoring function for generating and evaluating ligand poses. In this review, we present the search algorithms and scoring functions most commonly used in current molecular docking methods that focus on protein–ligand applications. We summarise the main topics and recent computational and methodological advances in protein–ligand docking. Protein flexibility, multiple ligand binding modes and the free-energy landscape profile for binding affinity prediction are important and interconnected challenges to be overcome by further methodological developments in the docking field.     

Protein-Protein Docking with Dynamic Residue Protonation States

Protein-protein interactions depend on a host of environmental factors. Local pH conditions influence the interactions through the protonation states of the ionizable residues that can change upon binding. In this work, we present a pH-sensitive docking approach, pHDock, that can sample side-chain protonation states of five ionizable residues (Asp, Glu, His, Tyr, Lys) on-the-fly during the docking simulation. pHDock produces successful local docking funnels in approximately half (79/161) the protein complexes, including 19 cases where standard RosettaDock fails. pHDock also performs better than the two control cases comprising docking at pH 7.0 or using fixed, predetermined protonation states. On average, the top-ranked pHDock structures have lower interface RMSDs and recover more native interface residue-residue contacts and hydrogen bonds compared to RosettaDock. Addition of backbone flexibility using a computationally-generated conformational ensemble further improves native contact and hydrogen bond recovery in the top-ranked structures. Although pHDock is designed to improve docking, it also successfully predicts a large pH-dependent binding affinity change in the Fc–FcRn complex, suggesting that it can be exploited to improve affinity predictions. The approaches in the study contribute to the goal of structural simulations of whole-cell protein-protein interactions including all the environmental factors, and they can be further expanded for pH-sensitive protein design.     

Mechanism of auxin interaction with Auxin Binding Protein (ABP1): a molecular dynamics simulation study

Auxin Binding Protein 1 (ABP1) is ubiquitous in green plants. It binds the phytohormone auxin with high specificity and affinity, but its role in auxin-induced processes is unknown. To understand the proposed receptor function of ABP1 we carried out a detailed molecular modeling study. Molecular dynamics simulations showed that ABP1 can adopt two conformations differing primarily in the position of the C-terminus and that one of them is stabilized by auxin binding. This is in agreement with experimental evidence that auxin induces changes at the ABP1 C-terminus. Simulations of ligand egress from ABP1 revealed three main routes by which an auxin molecule can enter or leave the ABP1 binding site. Assuming the previously proposed orientation of ABP1 to plant cell membranes, one of the routes leads to the membrane and the other two to ABP1's aqueous surroundings. A network of hydrogen-bonded water molecules leading from the bulk water to the zinc-coordinated ligands in the ABP1 binding site was formed in all simulations. Water entrance into the zinc coordination sphere occurred simultaneously with auxin egress. These results suggest that the hydrogen-bonded water molecules may assist in protonation and deprotonation of auxin molecules and their egress from the ABP1 binding site.     

Some insights into the binding mechanism of the GABAA receptor: a combined docking and MM-GBSA study

Gamma-aminobutyric type A receptor (GABA_AR) is a member of the Cys-loop family of pentameric ligand gated ion channels (pLGICs). It has been identified as a key target for many clinical drugs. In the present study, we construct the structure of human 2α₁2β₂γ₂ GABA_AR using a homology modeling method. The structures of ten benzodiazepine type drugs and two non-benzodiazepine type drugs were then docked into the potential benzodiazepine binding site on the GABA_AR. By analyzing the docking results, the critical residues His102 (α₁), Phe77 (γ₂) and Phe100 (α₁) were identified in the binding site. To gain insight into the binding affinity, molecular dynamics (MD) simulations were performed for all the receptor–ligand complexes. We also examined single mutant GABA_AR (His102A) in complexes with the three drugs (flurazepam, eszopiclone and zolpidem) to elucidate receptor–ligand interactions. For each receptor–ligand complex (with flurazepam, eszopiclone and zolpidem), we calculated the average distance between the C_α of the mutant residue His102A (α₁) to the center of mass of the ligands. The results reveal that the distance between the C_α of the mutant residue His102A (α₁) to the center of flurazepam is larger than that between His102 (α₁) to flurazepam in the WT type complex. Molecular mechanic-generalized Born surface area (MM-GBSA)-based binding free energy calculations were performed. The binding free energy was decomposed into ligand-residue pairs to create a ligand-residue interaction spectrum. The predicted binding free energies correlated well (R ²?=?0.87) with the experimental binding free energies. Overall, the major interaction comes from a few groups around His102 (α₁), Phe77 (γ₂) and Phe100 (α₁). These groups of interaction consist of at least of 12 residues in total with a binding energy of more than 1 kcal mol^?1. The simulation study disclosed herein provides a meaningful insight into GABA_AR–ligand interactions and helps to arrive at a binding mode hypothesis with implications for drug design.     

Investigating the Putative Binding-mode of GABA and Diazepam within GABA<Subscript>A</Subscript> Receptor Using Molecular Modeling

The three-dimensional structure of the GABA A receptor that included the ligand/agonist binding site was constructed and validated by using molecular modeling technology. Moreover, the putative binding-mode of GABA and diazepam with GABAA receptor were investigated by means of docking studies. Based on an rmsd-tolerance of 1.0 angstroms, the docking of GABA to alpha1/beta2 interface resulted in three multi-member conformational clusters and model 2 was supported by homologous sequence alignment data and experimental evidence. On the other hand, the docking of diazepam to alpha1/gamma2 interface revealed five multi-member conformational clusters in the binding site and model 1 seemed to represent the correct orientation of diazepam in the binding site.     

Computational analysis of plasmepsin IV bound to an allophenylnorstatine inhibitor

The plasmepsin proteases from the malaria parasite Plasmodium falciparum are attracting attention as putative drug targets. A recently published crystal structure of Plasmodium malariae plasmepsin IV bound to an allophenylnorstatine inhibitor [Clemente, J.C. et al. (2006) Acta Crystallogr. D 62, 246-252] provides the first structural insights regarding interactions of this family of inhibitors with plasmepsins. The compounds in this class are potent inhibitors of HIV-1 protease, but also show nM binding affinities towards plasmepsin IV. Here, we utilize automated docking, molecular dynamics and binding free energy calculations with the linear interaction energy LIE method to investigate the binding of allophenylnorstatine inhibitors to plasmepsin IV from two different species. The calculations yield excellent agreement with experimental binding data and provide new information regarding protonation states of active site residues as well as conformational properties of the inhibitor complexes.     

Novel statistical-thermodynamic methods to predict protein-ligand binding positions using probability distribution functions

We present two novel methods to predict native protein-ligand binding positions. Both methods identify the native binding position as the most probable position corresponding to a maximum of a probability distribution function (PDF) of possible binding positions in a protein active site. Possible binding positions are the origins of clusters composed, on the basis of root-mean square deviations (RMSD), from the multiple ligand positions determined by a docking algorithm. The difference between the methods lies in the ways the PDF is derived. To validate the suggested methods, we compare the averaged RMSD of the predicted ligand docked positions relative to the experimentally determined positions for a set of 135 PDB protein-ligand complexes. We demonstrate that the suggested methods improve docking accuracy by as much as 21-24% in comparison with a method that simply identifies the binding position as the energy top-scored ligand position.     

Exploring the binding pocket for pyridopyrimidine ligands at the CCK1 receptor by molecular docking

Pyridopyrimidine-based analogues are among the most highly potent and selective antagonists of cholecystokinin receptor subtype-1 (CCK1R) described to date. To better understand the structural and chemical features responsible for the recognition mechanism, and to explore the binding pocket of these compounds, we performed automated molecular docking using GOLD2.2 software on some derivatives with structural diversity, and propose a putative binding conformation for each compound. The docking protocol was guided by the key role of the Asn333 residue, as revealed by site directed mutagenesis studies. The results suggest two putative binding modes located in the same pocket. Both are characterized by interaction with the main residues revealed by experiment, Asn333 and Arg336, and differ in the spatial position of the Boc-Trp moiety of these compounds. Hydrophobic contacts with residues Thr117, Phe107, Ile352 and Ile329 are also in agreement with experimental data. Despite the poor correlation obtained between the estimated binding energies and the experimental activity, the proposed models allow us to suggest a plausible explanation of the observed binding data in accordance with chemical characteristics of the compounds, and also to explain the observed diastereoselectivity of this family of antagonists towards CCK1R. The most reasonable selected binding conformations could be the starting point for future studies. Figure Superimposition of the two putative binding conformations revealed by molecular docking for pyridopyrimidine-based CCK1 antagonists     

Molecular basis of the recognition of arachidonic acid by cytochrome P450 2E1 along major access tunnel

Cytochrome P450 2E1 is widely known for its ability to oxidize both low molecular weight xenobiotics and endogenous fatty acids (e.g., arachidonic acid (AA)). In this study, we investigated the structural features of the AA‐bound CYP2E1 complex utilizing molecular dynamics (MD) and found that the distinct binding modes for both AA and fatty acid analog are conserved. Moreover, multiple random acceleration MD simulations and steered MD simulations uncovered the most possible tunnel for fatty acids. The main attractions are derived from three key residues, His107, Ala108, and His109, whose side chains reorient to keep ligands bound via hydrogen bonds during the initial unbinding process. More importantly, based on the calculated binding free energy results, we hypothesize that the hydrogen bonds between the receptor and the ligand are the most important contributors involved in the binding affinity. Thus, it is inferred that the hydrogen bonds between these three residues and the ligand may help offer insights into the structural basis of the different ligand egress mechanisms for fatty acids and small weight compounds. Our investigation provides detailed atomistic insights into the structural features of human CYP2E1–fatty acid complex structures. Furthermore, the ligand‐binding characteristics obtained in the present study are helpful for both experimental and computational studies of CYPs and may allow future researchers to achieve desirable changes in enzymatic activities. © 2014 Wiley Periodicals, Inc. Biopolymers 103: 53–66, 2015.