Modeling loop backbone flexibility in receptor‐ligand docking simulations

The relevance of receptor conformational change during ligand binding is well documented for many pharmaceutically relevant receptors, but is still not fully accounted for in in silico docking methods. While there has been significant progress in treatment of receptor side chain flexibility sampling of backbone flexibility remains challenging because the conformational space expands dramatically and the scoring function must balance protein–protein and protein–ligand contributions. Here, we investigate an efficient multistage backbone reconstruction algorithm for large loop regions in the receptor and demonstrate that treatment of backbone receptor flexibility significantly improves binding mode prediction starting from apo structures and in cross docking simulations. For three different kinase receptors in which large flexible loops reconstruct upon ligand binding, we demonstrate that treatment of backbone flexibility results in accurate models of the complexes in simulations starting from the apo structure. At the example of the DFG‐motif in the p38 kinase, we also show how loop reconstruction can be used to model allosteric binding. Our approach thus paves the way to treat the complex process of receptor reconstruction upon ligand binding in docking simulations and may help to design new ligands with high specificity by exploitation of allosteric mechanisms. © 2012 Wiley Periodicals, Inc.     

Intervals and the deduction of drug binding site models

In the search for new drugs, it often occurs that the binding affinities of several compounds to a common receptor macromolecule are known experimentally, but the structure of the receptor is not known. This article describes an extraordinarily objective computer algorithm for deducing the important geometric and energetic features of the common binding site, starting only from the chemical structures of the ligands and their observed binding. The user does not have to propose a pharmacophore, guess the bioactive conformations of the ligands, or suggest ways to superimpose the active compounds. The method takes into account conformational flexibility of the ligands, stereospecific binding, diverse or unrelated chemical structures, inaccurate or qualitative binding data, and the possibility that chemically similar ligands may or may not bind to the receptor in similar orientations. The resulting model can be viewed graphically and interpreted in terms of one or more binding regions of the receptor, each preferring to be occupied by various sorts of chemical groups. The model always fits the given data completely and can predict the binding of any other ligand, regardless of chemical structure. The method is an outgrowth of distance geometry and Voronoi polyhedra site modeling but incorporates several novel features. The geometry of the ligand molecules and the site is described in terms of intervals of internal distances. Determining the site model consists of reducing the uncertainty in the interregion distance intervals, and this uncertainty is described as intervals of intervals. Similarly, the given binding affinities and their experimental uncertainties are treated as intervals in the affinity scale. The final site model specifies an entire region of interaction energy parameters that satisfy the training set rather than a single set of parameters. Predicted binding for test compounds results in an interval which, when compared to the experimental interval, may be correct, incorrect, or vague. There is a pervasive ternary logic involved in the assessment of predictions, in the search for a satisfactory model, and in judging whether a given molecule may bind in a particular orientation: true, false, or maybe. The approach is illustrated on an extremely simple artificial example and on a real data set of cocaine analogues binding to a nerve membrane receptor in vitro. © 1995 by John Wiley & Sons, Inc.     

Fully automated molecular mechanics based induced fit protein-ligand docking method

We describe a method for docking a ligand into a protein receptor while allowing flexibility of the protein binding site. The method employs a multistep procedure that begins with the generation of protein and ligand conformations. An initial placement of the ligand is then performed by computing binding site hotspots. This initial placement is followed by a protein side-chain refinement stage that models protein flexibility. The final step of the process is an energy minimization of the ligand pose in the presence of the rigid receptor. Thus the algorithm models flexibility of the protein at two stages, before and after ligand placement. We validated this method by performing docking and cross docking studies of eight protein systems for which crystal structures were available for at least two bound ligands. The resulting rmsd values of the 21 docked protein-ligand complexes showed values of 2 A or less for all but one of the systems examined. The method has two critical benefits for high throughput virtual screening studies. First, no user intervention is required in the docking once the initial binding site selection has been made in the protein. Second, the initial protein conformation generation needs to be performed only once for a given binding region. Also, the method may be customized in various ways depending on the particular scenario in which dockings are being performed. Each of the individual steps of the method is fully independent making it straightforward to explore different variants of the high level workflow to further improve accuracy and performance.     

HIV蛋白酶抑制剂——利迪链菌素的分子对接研究

用分子对接的方法, 对利迪链菌素的抗HIV蛋白酶活性进行了研究. 为了更准确地反映利迪链菌素分子与酶蛋白结合的情况, 充分考虑受体活性部位的柔性, 采用了FlexX(初步对接)和Flexidock(精确对接)分两步将配体与受体进行对接. 在初步对接中, 设计了不同的受体活性部位来考察是否有结合水分子参与抑制剂与酶的结合. 对一种作用方式已知的非肽类HIV蛋白酶抑制剂Aha006进行的对接研究显示, 分子模拟的结果与实际情况吻合得较好, 证明了本文所采用的方法的可靠性. 利迪链菌素与蛋白酶活性部位的对接结果显示, 配体分子与受体之间的结合没有结合水分子的参与, 两者通过5对氢键作用结合成为稳定的复合物. 利迪链菌素占据结合腔, 覆盖了蛋白酶的活性三联体Asp25-Thr26-Gly27, 从而起到抑制其生物活性的作用.     

Pharmacophore and receptor models for neurokinin receptors

Three neurokinin (NK) antagonist pharmacophore models (Models 1-3) accounting for hydrogen bonding groups in the 'head' and 'tail' of NK receptor ligands have been developed by use of a new procedure for treatment of hydrogen bonds during superimposition. Instead of modelling the hydrogen bond acceptor vector in the strict direction of the lone pair, an angle is allowed between the hydrogen bond acceptor direction and the ideal lone pair direction. This approach adds flexibility to hydrogen bond directions and produces more realistic RMS values. By using this approach, two novel pharmacophore models were derived (Models 2 and 3) and a hydrogen bond acceptor was added to a previously published NK2 pharmacophore model [Poulsen et al., J. Comput.-Aided Mol. Design, 16 (2002) 273] (Model 1). Model 2 as well as Model 3 are described by seven pharmacophore elements: three hydrophobic groups, three hydrogen bond acceptors and a hydrogen bond donor. Model 1 contains the same hydrophobic groups and hydrogen bond donor as Models 2 and 3, but only one hydrogen bond acceptor. The hydrogen bond acceptors and donor are represented as vectors. Two of the hydrophobic groups are always aromatic rings whereas the other hydrophobic group can be either aromatic or aliphatic. In Model 1 the antagonists bind in an extended conformation with two aromatic rings in a parallel displaced and tilted conformation. Model 2 has the same two aromatic rings in a parallel displaced conformation whereas Model 3 has the rings in an edge to face conformation. The pharmacophore models were evaluated using both a structure (NK receptor homology models) and a ligand based approach. By use of exhaustive conformational analysis (MMFFs force field and the GB/SA hydration model) and least-squares molecular superimposition studies, 21 non-peptide antagonists from several structurally diverse classes were fitted to the pharmacophore models. More antagonists could be fitted to Model 2 with a low RMS and a low conformational energy penalty than to Models 1 and 3. Pharmacophore Model 2 was also able to explain the NK1, NK2 and NK3 subtype selectivity of the compounds fitted to the model. Three NK 7TM receptor models were constructed, one for each receptor subtype. The location of the antagonist binding site in the three NK receptor models is identical. Compounds fitted to pharmacophore Model 2 could be docked into the NK1, NK2 and NK3 receptor models after adjustment of the conformation of the flexible linker connecting the head and tail. Models I and 3 are not compatible with the receptor models.     

An automated method for predicting the positions of hydrogen-bonding atoms in binding sites

Hydrogen bonds are the most specific, and therefore predictable of the intermolecular interactions involved in ligand–protein binding. Given the structure of a molecule, it is possible to estimate the positions at which complementary hydrogen-bonding atoms could be found. Crystal-survey data are used in the design of a program, HBMAP, that generates a hydrogen-bond map for any given ligand, which contains all the feasible positions at which a complementary atom could be found. On superposition of ligands, the overlapping regions of their maps represent positions of receptor atoms to which each molecule can bind. The certainty of these positions is increased by the incorporation of a larger number and diversity of molecules. In this work, superposition is achieved using the program HBMATCH, which uses simulated annealing to generate the correspondence between points from the hydrogen-bonding maps of the two molecules. Equivalent matches are distinguished on the basis of their steric similarity. The strategy is tested on a number of ligands for which ligand–protein complexes have been solved crystallographically, which allows validation of the techniques. The receptor atom positions of thermolysin are successfully predicted when the correct superposition is obtained.     

Computational methods for the structural alignment of molecules

In drug design, often enough, no structural information on a particular receptor protein is available. However, frequently a considerable number of different ligands is known together with their measured binding affinities towards a receptor under consideration. In such a situation, a set of plausible relative superpositions of different ligands, hopefully approximating their putative binding geometry, is usually the method of choice for preparing data for the subsequent application of 3D methods that analyze the similarity or diversity of the ligands. Examples are 3D-QSAR studies, pharmacophore elucidation, and receptor modeling. An aggravating fact is that ligands are usually quite flexible and a rigorous analysis has to incorporate molecular flexibility. We review the past six years of scientific publishing on molecular superposition. Our focus lies on automatic procedures to be performed on arbitrary molecular structures. Methodical aspects are our main concern here. Accordingly, plain application studies with few methodical elements are omitted in this presentation. While this review cannot mention every contribution to this actively developing field, we intend to provide pointers to the recent literature providing important contributions to computational methods for the structural alignment of molecules. Finally we provide a perspective on how superposition methods can effectively be used for the purpose of virtual database screening. In our opinion it is the ultimate goal to detect analogues in structure databases of nontrivial size in order to narrow down the search space for subsequent experiments.     

Improving binding mode and binding affinity predictions of docking by ligand-based search of protein conformations: evaluation in D3R grand challenge 2015

The growing number of protein–ligand complex structures, particularly the structures of proteins co-bound with different ligands, in the Protein Data Bank helps us tackle two major challenges in molecular docking studies: the protein flexibility and the scoring function. Here, we introduced a systematic strategy by using the information embedded in the known protein–ligand complex structures to improve both binding mode and binding affinity predictions. Specifically, a ligand similarity calculation method was employed to search a receptor structure with a bound ligand sharing high similarity with the query ligand for the docking use. The strategy was applied to the two datasets (HSP90 and MAP4K4) in recent D3R Grand Challenge 2015. In addition, for the HSP90 dataset, a system-specific scoring function (ITScore2_hsp90) was generated by recalibrating our statistical potential-based scoring function (ITScore2) using the known protein–ligand complex structures and the statistical mechanics-based iterative method. For the HSP90 dataset, better performances were achieved for both binding mode and binding affinity predictions comparing with the original ITScore2 and with ensemble docking. For the MAP4K4 dataset, although there were only eight known protein–ligand complex structures, our docking strategy achieved a comparable performance with ensemble docking. Our method for receptor conformational selection and iterative method for the development of system-specific statistical potential-based scoring functions can be easily applied to other protein targets that have a number of protein–ligand complex structures available to improve predictions on binding.     

Time-efficient flexible superposition of medium-sized molecules

We present an efficient algorithm for the structural alignment of medium-sized organic molecules. The algorithm has been developed for applications in 3D QSAR and in receptor modeling. The method assumes one of the molecules, the reference ligand, to be presented in the conformation that it adopts inside the receptor pocket. The second molecule, the test ligand, is considered to be flexible, and is assumed to be given in an arbitrary low-energy conformation. Ligand flexibility is modeled by decomposing the test ligand into molecular fragments, such that ring systems are completely contained in a single fragment. Conformations of fragments and torsional angles of single bonds are taken from a small finite set, which depends on the fragment and bond, respectively. The algorithm superimposes a distinguished base fragment of the test ligand onto a suitable region of the reference ligand and then attaches the remaining fragments of the test ligand in a step-by-step fashion. During this process, a scoring function is optimized that encompasses bonding terms and terms accounting for steric overlap as well as for similarity of chemical properties of both ligands. The algorithm has been implemented in the FLEXS system. To validate the quality of the produced results, we have selected a number of examples for which the mutual superposition of two ligands is experimentally given by the comparison of the binding geometries known from the crystal structures of their corresponding protein–ligand complexes. On more than two-thirds of the test examples the algorithm produces rms deviations of the predicted versus the observed conformation of the test ligand below 1.5 Å. The run time of the algorithm on a single problem instance is a few minutes on a common-day workstation. The overall goal of this research is to drastically reduce run times, while limiting the inaccuracies of the model and the computation to a tolerable level.     

Geometrically feasible binding modes of a flexible ligand molecule at the receptor site

A very efficient algorithm for determining the geometrically feasible binding modes of a flexible ligand molecule at the receptor site is presented. It is based on distance geometry but maintains the requirements of three dimensions. The distance geometry manipulation can superimpose two bodies without explicitly calculating the necessary rigid rotation and translation. The whole conformation space of a flexible molecule can be efficiently examined by considering only a finite number of conformational points. The method is suitable only when the criterion for superposition is some minimum distance limit. It cannot, however, give the exact distance between two points in two different bodies.     

Improved FlexX docking using FlexS-determined base fragment placement

We report on a novel hybrid FlexX/FlexS docking approach, whereby the base fragment of the test ligand is chosen by FlexS superposition onto a cocrystallized template ligand and then fed into FlexX for the incremental construction of the final solution. The new approach is tested on the diverse 200 protein-ligand complex dataset that has been previously described for FlexX validation. In total, 62.9% of the complexes can be reproduced at rank 1 by our approach, which compares favorably with 46.9% when using FlexX alone. In addition, we report "cross-docking" experiments in which several receptor structures of complexes with identical proteins have been used for docking all cocrystallized ligands of these complexes. The results show that, in almost all cases, the hybrid approach can acceptably dock a ligand into a foreign receptor structure using a different ligand template, can give solutions where FlexX alone fails, and tends to give solutions that are more accurately positioned.     

Molecular modelling studies on the ORL1-receptor and ORL1-agonists

The ORLI (opioid receptor like 1)- receptor is a member of the family of rhodopsin-like G protein-coupled receptors (GPCR) and represents an interesting new therapeutical target since it is involved in a variety of biomedical important processes, such as anxiety, nociception, feeding, and memory. In order to shed light on the molecular basis of the interactions of the GPCR with its ligands, the receptor protein and a dataset of specific agonists were examined using molecular modelling methods. For that purpose, the conformational space of a very potent non-peptide ORL1-receptor agonist (Ro 64-6198) with a small number of rotatable bonds was analysed in order to derive a pharmacophoric arrangement. The conformational analyses yielded a conformation that served as template for the superposition of a set of related analogues. Structural superposition was achieved by employing the program FlexS. Using the experimental binding data and the superposition of the ligands, a 3D-QSAR analysis applying the GRID/GOLPE method was carried out. After the ligand-based modelling approach, a 3D model of the ORL1-receptor has been constructed using homology modelling methods based on the crystal structure of bovine rhodopsin. A representative structure of the model taken from molecular dynamics simulations was used for a manual docking procedure. Asp-130 and Thr-305 within the ORL1-receptor model served as important hydrophilic interaction partners. Furthermore, a hydrophobic cavity was identified stabilizing the agonists within their binding site. The manual docking results were supported using FlexX, which identified the same protein-ligand interaction points.     

Modelling of the binding site of the human m1 muscarinic receptor: Experimental validation and refinement

Our model of the human m1 muscarinic receptor has been refined on the basis of the recently published projection map of bovine rhodopsin. The refined model has a slightly different helix arrangement, which reveals the presence of an extra hydrophobic pocket located between helices 3, 4 and 5. The interaction of series of agonists and antagonists with the m1 muscarinic receptor has been studied experimentally by site-directed mutagenesis. In order to account for the observed results, three-dimensional models of m1 ligands docked in the target receptor are proposed. Qualitatively, the obtained models are in good agreement with the experimental observations. Agonists and partial agonists have a relatively small size. They can bind to the same region of the receptor using, however, different anchoring receptor residues. Antagonists are usually larger molecules, filling almost completely the same pocket as agonists. They can usually produce much stronger interactions with aromatic residues. Experimental data combined with molecular modelling studies highlight how subtle and diverse receptor–ligand interactions could be.     

Electrostatic complementarity between proteins and ligands. 3. Structural basis

Summary Electrostatic potential complementarity between ligands and their receptor sites is evaluated by the superposition of the electrostatic potential, generated by the receptor, onto the ligand potential over the ligand van der Waals surface. We would like to examine which structural factors generate this pattern of superposition. Example studies suggest that in many ligand-protein pairs, there exist principal formal charges on each molecule, largely responsible for the electrostatic potential complementarity observed. Electrostatic potential complementarity depends on the relative disposition of these principal charges and the ligand van der Waals surface. Simple mathematical models were constructed to predict the complementarity solely from structural considerations. The essential conditions for electrostatic potential complementarity were elucidated. These can be used in ligand design strategies to obtain an electrostatically optimal ligand.     

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.     

Harmonic modes as variables to approximately account for receptor flexibility in ligand–receptor docking simulations: Application to DNA minor groove ligand complex

An approach to approximately account for receptor flexibility in ligand–receptor docking simulations is described and applied to a DNA/Hoechst 33258 analogue complex. Harmonic modes corresponding to eigenvectors with small eigenvalues of the Hessian matrix of the potential energy function were used as independent variables to describe receptor flexibility. For the DNA minor groove ligand case most of the conformational difference between an energy minimized free DNA and ligand-bound structure could be assigned to 5–40 harmonic receptor modes with small eigenvalues. During docking, deformations of the DNA receptor structure in the subset of harmonic modes were limited using a simple penalty function that avoided the summation over all intrareceptor atom pairs. Significant improvement of the sterical fit between ligand and receptor was found upon relaxation of the DNA in the subset of harmonic modes after docking of the ligand at the position found in the known crystal structure. In addition, the harmonic mode relaxation resulted in DNA structures that were more similar to the energy minimized ligand-bound form. Although harmonic mode relaxation also leads to improved sterical fit for other ligand placements, the placement as observed in the crystal structure could still be identified as the site with the most favorable sterical interactions. Because relaxation in the harmonic modes is orders of magnitude faster than conventional energy minimization using all atom coordinates as independent variables, the approach might be useful as a preselection tool to recognize ligand binding sites accessible only upon small conformational changes of the receptor. The harmonic mode relaxed structures can only be considered as approximate structures because deformation of the receptor in the harmonic modes can lead to small perturbations of the stereochemical geometry of the molecule. Energy minimization of preselected ligand–DNA docking candidates in all atom coordinates is required to reduce these deviations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 287–300, 1999     

Conformational selection of protein kinase A revealed by flexible-ligand flexible-protein docking

Protein kinases have high structural plasticity: their structure can change significantly, depending on what ligands are bound to them. Rigid-protein docking methods are not capable of describing such effects. Here, we present a new flexible-ligand flexible-protein docking model in which the protein can adopt conformations between two extremes observed experimentally. The model utilized a molecular dynamics-based simulated annealing cycling protocol and a distance-dependent dielectric model to perform docking. By testing this model on docking four diverse ligands to protein kinase A, we found that the ligands were able to dock successfully to the protein with the proper conformations of the protein induced. By imposing relatively soft conformational restraints to the protein during docking, this model reduced computational costs yet permitted essential conformational changes that were essential for these inhibitors to dock properly to the protein. For example, without adequate movement of the glycine-rich loop, it was difficult for the ligands to move from the surface of the protein to the binding site. In addition, these simulations called for better ways to compare simulation results with experiment other than using the popular root-mean-square deviation between the structure of a ligand in a docking pose and that in experiment because the structure of the protein also changed. In this work, we also calculated the correlation coefficient between protein-ligand/protein-protein distances in the docking structure and those in the crystal structure to check how well a ligand docked into the binding site of the protein and whether the proper conformation of the protein was induced.