Learn more: PMC Disclaimer | PMC Copyright Notice
Ligand Docking Methods to Recognize Allosteric Inhibitors for G-Protein-Coupled Receptors
Associated Data
Abstract
G-protein-coupled receptors (GPCRs) are membrane proteins which play an important role in many cellular processes and are excellent drug targets. Despite the existence of several US Food and Drug Administration (FDA)-approved GPCR-targeting drugs, there is a continuing challenge of side effects owing to the nonspecific nature of drug binding. We have investigated the diversity of the ligand binding site for this class of proteins against their cognate ligands using computational docking, even if their structures are known already in the ligand-complexed form. The cognate ligand of some of these receptors dock at allosteric binding site with better score than the binding at the conservative site. Interestingly, amino acid residues at such allosteric binding site are not conserved across GPCR subfamilies. Such a computational approach can assist in the prediction of specific allosteric binders for GPCRs.
Introduction
Membrane proteins constitute nearly 40% of the human genome.1,2 G-protein-coupled receptors (GPCRs), being important drug targets, are one of the most well studied3-5 among the membrane proteins. There are close to 800 known GPCRs in the human genome, and they are grouped into 5 families depending on their substrates (such as peptides and amines 6 ). GPCRs act as an important environmental sensor and play a vital role in diverse signalling processes. As a result, there is enormous diversity in their position in different biological pathways.
The hallmark of GPCRs is the presence of 7 transmembrane helices (TMHs), where extracellular loop regions and parts of TMHs recognize the cognate substrates and indentations in the intracellular loop regions provide the capacity to participate in diverse biological pathways. 7 In the past 2 decades, there have been substantial insights into the structural features of few important GPCRs, despite the inherent challenges in the structure determination of membrane proteins. Consequently, we now have a great deal of information about the conformational changes that occur subsequent to ligand binding, distinct sites for agonist and antagonist binding, and so on. Characteristic functional motifs in the intracellular regions, the presence of conserved charged residues in the intracellular face and the presence of prolyl residues are known to contribute to signalling and conformational changes. 8 We now have close to 389 Protein Data Bank (PDB) entries (also organized in specialized GPCR-EXP database https://zhanglab.dcmb.med.umich.edu/GPCR-EXP/) 9 which pertain to GPCRs, either in the apo-form or ligand-bound form (agonist or antagonist-bound states)). There has been successful design of GPCR drugs, such as haloperidol, 10 over the years. As of November 2017, the US Food and Drug Administration (FDA) has approved drugs against 134 GPCRs. 11
Despite the structural insights available, drug design for GPCRs remains highly challenging due to the inherent characteristics of substrate promiscuity, structural similarities of drug molecules and sequence similarity within subfamilies of GPCRs. There have been continuing demands to provide drug solutions to address this prevailing feature of GPCRs. Mutations in either GPCRs or their interacting proteins have been widely implicated in neurodegenerative and other diseases.12,13 Some of the important GPCR drug targets are in the area of neurodegenerative diseases, where there is appalling amounts of side effects observed in patients who are treated with GPCR drugs. Recent efforts have, therefore, focused on the design of allosteric inhibitors.14-16
In this article, we report systematic computational ligand docking experiments using selected GPCRs with known information on ligand binding, to present which factors enable best capture of near-native ligand binding and how this can be employed subsequently to identify allosteric inhibitors. Our choice of ligand docking algorithm is AutoDock which is a well-known software that has gone through rigorous analyses by other groups as well. 17 We first describe the analysis of binding poses, where different levels of known information on structural or evolutionary conservation can guide the docking process. Objective measures, such as Tanimoto co-efficient, have been employed to assess the comparison between docking poses and the structural data for ligand binding. This approach enables the identification of structurally reasonable docking poses which are not close to the native pose, namely allosteric binding. Finally, to test this approach, we have performed blind docking on one of the GPCRs, chemokine receptors, using specialized set of ligands known to bind GPCRs. 18 Few allosteric binding sites could indeed be recognized using this novel computational approach. Interactions between ligand and chemokine receptor at the allosteric site were found to be stable, when subject to molecular dynamics (MD) simulations.
Methods
Selection of GPCR structures and ligands for analysis
All the available GPCR structures were downloaded from PDB (April 2017). 19 Among 71 such structures, several identical GPCR structures with different resolution and agonist/antagonist bound forms were observed. Thus, for every GPCR with bound ligand, one best structure with highest resolution was selected to avoid redundancy. A set of 27 GPCRs were selected for final analysis (Table 1). G-protein-coupled receptor alignments were obtained from GPCRdb for the subfamily of GPCR receptors to study evolutionary conservation. 20 Ligands were retrieved from PDB, PubChem, 21 or GLASS database. 18 GLASS database is a resource to retrieve ligands that are known to bind to GPCRs and were used to explore allosteric binding. Out of 873 ligands (Supplementary Table 1), 46 were chosen based on Lipinski’s drug likeness and XlogP (less than 2). The 46 ligands were clustered using ChemMine tool 22 and finally 29 compounds were used for docking (Supplementary File). These compounds were converted to PDB format using Open babel tool. 23
Table 1.
List of GPCRs used in the study.
| Serial no. | GPCR name | PDB ID |
|---|---|---|
| GPCR-1 | 5-HT1B receptor | 4IAQ |
| GPCR-2 | 5-HT2B receptor | 4IB4 |
| GCPR-3 | A2A receptor | 2YDO |
| GPCR-4 | M2 receptor | 3UON |
| GPCR-5 | M3 receptor | 4DAJ |
| GPCR-6 | β1-adrenoceptor | 2VT4 |
| GPCR-7 | β2-adrenoceptor | 2RH1 |
| GPCR-8 | AT1 receptor | 4YAY |
| GPCR-9 | CCR5 | 4MBS |
| GPCR-10 | CRF1 receptor | 4K5Y |
| GPCR-11 | CXCR4 | 3ODU |
| GPCR-12 | D3 receptor | 3PBL |
| GPCR-13 | FFA1 receptor | 4PHU |
| GPCR-14 | mGlu5 receptor | 4OO9 |
| GPCR-15 | H1 receptor | 3RZE |
| GPCR-16 | LPA1 receptor | 4Z34 |
| GPCR-17 | δ receptor | 4BUO |
| GPCR-18 | κ receptor | 4N6H |
| GPCR-19 | μ receptor | 4DJH |
| GPCR-20 | NOP receptor | 4DKL |
| GPCR-21 | Rhodopsin | 4EA3 |
| GPCR-22 | OX2 receptor | 4S0V |
| GPCR-23 | P2Y1 receptor | 4XNV |
| GPCR-24 | P2Y12 receptor | 4NTJ |
| GPCR-25 | PAR1 | 3VW7 |
| GPCR-26 | S1P1 receptor | 3V2Y |
| GPCR-27 | SMO | 4JKV |
Abbreviations: GPCR, G-protein-coupled receptors; PDB, Protein Data Bank.
GPCR proteins were selected based on the highest-resolution structure available in PDB. All class-A GPCRs structures were selected from PDB and a nonredundant data set was obtained.
Grid setting for GPCR-ligand docking
The protein and ligand coordinates were separated from PDB file and used for docking. In 2 out of 3 experiments, ligand coordinates were derived from PubChem (please see below). Both the agonist and the antagonist bind in the extracellular domains for most of the GPCRs in distinct sites. Hence, the grid region of choice for docking was chosen in the extracellular region. However, different grid settings were considered to identify the parameters in grid setting that will help in correct prediction of ligand binding to GPCRs. These are as follows:
Level 1: Ligand coordinates were taken as such from PDB, but ligand coordinates were kept flexible for docking. Only the extracellular half of the receptor region was selected as grid box (semiblind docking) (Figure 1A).
Different grids used for docking of ligands to GPCRs. (A) For Level 1, extracellular half of GPCRs was considered for grid box generation and ligand coordinates were taken from PDB. (B) For Level 2, smaller grid box covering only the known ligand-binding site was defined, and ligand coordinates were retrieved from PubChem. (C) For Level 3, the ligands coordinates were retrieved from PubChem while the grid box was defined by considering the extracellular half region of GPCRs.
GPCR indicates G-protein-coupled receptors.
Level 2: Ligand coordinates were taken from PubChem and ligand was kept flexible. A smaller grid that covers just the known ligand binding site was selected as grid box (guided docking) (Figure 1B).
Level 3: Ligand coordinates were taken from PubChem and ligand was kept flexible. Only the extracellular half of the receptor region was selected as grid box (semiblind docking) (Figure 1C).
Molecular docking
AutoDock 4.2.6 was used for docking studies. 17 Protein and ligand structures were prepared using ADT (AutoDock Tool). The hydrogen atoms were added and water molecules were removed followed by addition of gasteiger charges. The grid parameter file was generated with default distance (0.375 Å) between grid points. Grid space was defined around the extracellular site of GPCRs. The search parameter was set to 100 Genetic Algorithm (GA) runs and Lamarckian genetic algorithm was used. Autogrid was run followed by AutoDock. Similar parameters were used for all ligands docking to a single protein. The docking results were analysed, and protein-ligand complex was made using ADT.
Calculation of Tanimoto co-efficient and identification of allosteric sites
To identify how similar the protein-ligand docking results are as observed in the X-ray crystallographic structure of complexes, we used the Tanimoto co-efficient score. Any atom within 4 Å distance was considered to be a contact with the ligand.
Let ‘a’ be the contacts in PDB structure and ‘b’ be the contacts in the AutoDock complex then the Tanimoto co-efficient (TC) is calculated as
Furthermore, for each ligand, all 100 poses were checked for their binding site compared to native ligand. The minimum distance between all the atoms of native ligand and all the atoms of docked ligand was calculated and if the minimum distance is more than 5 Å for any pose, then that pose was considered as allosteric site-binding pose.
MD simulation
The protein-ligand complex structure was subjected to MD simulation using Desmond module of Schrodinger to assess the interaction stability of complex. 24 Initially, the complex structure was processed using protein preparation wizard of maestro (Schrodinger Release 2019-4: Maestro, Schrodinger, LLC, New York, NY, 2019) which assigns bond orders, sets protonation state, optimizes H-bonds and then minimizes the structure. 25 The structure (CXCR4-ligand complex) was protonated at pH 5.5 before minimization in protein preparation wizard. The membrane information was retrieved from OPM database. 26 After processing, the structure was solvated in the presence of membrane (POPC) using TIP4P water model. Orthorhombic box, with buffer distance of 10 Å, was used followed by minimization of box size. The system was neutralized and 150 mM salt (NaCl) was added. System builder was run and output of system builder (solvated system) was used for MD simulation. The system was subjected to default relaxation protocol of Desmond followed by production MD run for 100 ns. Simulation results were analysed using Simulation Interaction Diagram and Simulation Event analysis modules of Schrodinger. The analysis was performed for entire range of simulation time. Root-mean-square deviation (RMSD) is calculated for each frame by aligning the complex to protein backbone of the reference frame. Significantly higher values of "Lig fit Prot” than protein RMSD signifies the diffusion of ligand away from its initial binding site. “Lig fit lig RMSD” is calculated by aligning the ligand on the reference ligand conformation and it indicates the internal fluctuation of ligand.
Results
Prior knowledge of the binding mode helps to predict correct protein-ligand interaction
Twenty-seven GPCRs from PDB were selected for analysis. Each of the protein-ligand complexes was separated, and docking was performed as mentioned in Methods. The 3 different grid conditions were used as mentioned in Methods (Level 1, 2, and 3) (Figure 1A-C). We observed that in 50% of the instances, when the ligand coordinates were taken from PDB and semiguided docking was performed, we could get the protein-ligand complex with high TC score with the original PDB structure and AutoDock energy. In Level 2 of docking (runs with use of ligand coordinates from PubChem and guided docking), complexes with high TC score and low energy could be obtained in 30% of instances. Level 3 (ligand from PubChem and semiguided docking), however, did not result in complexes similar to the original structure in most of the cases (only 3 out of 27 GPCRs had high TC scores) (Figure 2: TC score Vs Energy correlation plot for Level 1). This trend can be explained owing to the increased levels of difficulty in docking and lower clues we provide from the native pose of the ligand under question. In level 1, the position of the ligand is quite close to the native pose of the ligand (please see Figure 1).
Identifying novel binding modes using the TC score and energy values for GPCR-ligand complexes
We selected complexes which had low TC score (that indicates a different binding mode compared to original ligand binding site [hereafter referred as ‘OLBS’]) and low-and-negative docking energy (that indicates stable binding) on the basis of Level 1 results. First 3 highest scoring complexes were selected for detailed analysis of the docking poses. These include beta-2-adrenergic receptor, chemokine receptor, and 5-HT2B receptor. It was observed that for beta-2-adrenergic receptor and 5-HT2B receptors, the ligand binds in the same pocket as in original complex, but is displaced along Y-axis of the protein with a large overlap with OLBS (Figure 3A and andB).B). This results in low TC score for these complexes. For chemokine receptor, however, we could find novel binding sites for the cognate ligand with minimum overlap with OLBS (Figure 3C). Thus, we used chemokine receptor-ligand complex for further analysis on allosteric binding in GPCRs.
Novel ligand binding poses for high scoring complexes. (A) Ligand poses for beta-2 adrenergic receptor indicating the overlap between novel docked ligand pose (grey) with native ligand-binding pose (red) and best scoring ligand pose is shown in blue. (B) Ligand poses for 5HT2b receptor. The novel binding sites (grey) are close to the PDB-complex ligand binding site (red) and the best scoring pose (blue). (C) Ligand poses for chemokine receptor CXCR4. Different structures of CXCR4 are bound to different ligands in PDB. Native ligand-binding pose (red) and best scoring ligand pose is shown in blue. Other binding modes observed in PDB are represented by yellow, grey, and magenta colour ligand poses.
PDB indicates Protein Data Bank.
Residues at the alternate novel binding sites are not conserved
Twenty-three homologous sequences of chemokine receptors were used to identify conservation of residues at the endogenous (orthosteric) binding site and allosteric binding site. It was observed that 7 out of 12 residues at OLBS were conserved, while only 3 out of 12 residues were evolutionarily conserved at the allosteric site, even within the same subfamily of receptors. This clearly suggests that allosteric binding site is novel (Figure 4).
Identifying novel allosteric binders from GLASS database for chemokine receptor
In order to explore allosteric binders for CXCR4, ligands were retrieved from GLASS database and blind docking was initiated. The docking energy of best docking pose for ligands varied from −6.25 kcal/mol (for ligand ‘25178561’) to −11.27 kcal/mol (for ligand ‘483559’). Out of 29 ligands, the overall best pose (in terms of docking energy) for 6 ligands also belongs to the largest cluster. Similar analysis was performed for all the ligands. Interestingly, it was observed that all 100 poses are allosteric poses for ligand ‘76381’ (Figure 5). The best pose of ligand ‘76381’ in terms of best docking energy is also the best pose among the largest cluster. Therefore, this best pose for ‘76381’ ligand to CXCR4 was subjected to MD simulation to assess the interaction stability of the complex.
Bar graph indicating the percentage of each ligand conformation in either native site or allosteric site.
MD indicates molecular dynamics; RMSD, root-mean-square deviation.
The CXCR4-76381 complex was subjected to 100 ns simulation and the RMSD was calculated. Root-mean-square deviation plot indicates that the ligand moves from its original binding position during simulations (Figure 6A). Major interacting residues of CXCR4 are Tyr256 and Gln200 (Figure 6B). Most of the interactions involve H-bonding and atom-wise interactions of 76381 are shown (Figure 6C). All contacts with the CXCR4 are lost at around 40 ns but regains subsequently, as indicated by total contact plot (Figure 6D). H-bond analysis between 76381 and CXCR4 (including H-bond with water) was also performed, and it was found that at multiple instants during simulation, the H-bond contact was lost and re-formed (Figure 7). However, at no point during simulation, there was any overlap between the native ligand binding mode (ITD) and 76381 binding mode (Figure 8; Supplementary Video file https://drive.google.com/file/d/1GvO7XiVRjFz6S4H6e5sQwtCxP4VV3jBO/view?usp = sharing) suggesting that these are true allosteric binding sites.
MD simulation results of CXCR4-76381 complex. (A) RMSD plot of CXCR4 and 76381. (B) Interacting residues and interaction types of CXCR4 with 76381 over the course of simulation time. Normalized stacked bars indicate the fraction of simulation time for which a particular type of interaction was maintained. Values more than 1.0 suggests that the residue forms multiple interactions of same subtype with ligand (C) interactions of 76381 atoms with residues of CXCR4 along with type and duration of interactions. Interactions that persist for more than 10% of simulation time have been shown. If some residues form multiple interaction of same type with the same atom of ligand, then interaction value can be more than 100% (D) total number of contacts (H-bonds, Water bridges, Hydrophobic, Ionic) between CXCR4 and ligand 76381 throughout the simulation.
MD indicates molecular dynamics; RMSD, root-mean-square deviation.
H-bond interactions between 76381 ligand and CXCR4 (and water) throughout the simulation period. Simulation time is in picoseconds.
76381 poses with CXCR4 compared to native ligand (ITD) of CXCR4 at different time points (in nanoseconds) during simulation. Native ligand pose (ITD) is shown as green-blue sticks. Cyan: 0 ns; Magenta: 33 ns; Yellow: 55 ns; Pink: 61 ns; White: 77 ns; Violet: 85 ns; Orange: 90 ns; Green: 100 ns.
ITD, isothiourea derivative.
Discussion
G-protein-coupled receptors are important drug targets for several approved human drugs. Since all the GPCRs have similar 7-transmembrane structure, it becomes challenging to specifically design drugs to target particular GPCRs. Thus, designing specific allosteric modulators is an important and challenging task. 14 In this study, we first understand how well different protocols (Levels 1-3) play a role in modelling GPCR-ligand interaction. We found that accurate guidance of docking (smaller grid) is a crucial parameter and any biochemical knowledge on ligand binding sites (as in Level 1 and 2) greatly improves the chances of reproducing native ligand pose using computational docking methods. We selected structural entries of GPCRs with alternate binding sites and a specialized database of ligands known to bind GPCRs, as documented in GLASS database, for an in-depth study. Structural analysis with smaller set of ligands against different proteins has been reported in the literature as well.27-32 We performed docking to check the binding of ligands from GLASS database with GPCRs. We were able to obtain ligands that are reported as targets for particular GPCRs that bind better to the allosteric site as compared to the orthosteric binding site.
Conclusion
Identification or computational design of highly specific drugs for GPCRs has been a problem for many decades. This is due to similar overall shape of ligand binding site, existence of subfamilies, and broad sequence variation among GPCRs. Nonavailability of drugs, which is highly specific for one GPCR subfamily member, gives rise to serious consequences and side effects to patients. Thus, design strategies to invest on allosteric sites acquire prime importance. Our computational approach employs reliable ligand binding software to recognize nonorthosteric ligand-binding sites which appear energetically favourable. Furthermore, our sequence analysis of such GPCRs, showed that residues at the allosteric site are not as conserved as in the orthosteric binding site, which provides hope for achieving the specificity of a particular GPCR. Furthermore, for certain GPCRs, some of their known cognate ligands were predicted to have better binding preference towards the allosteric site than orthosteric site. We conclude that few of the known GPCR-specific ligands may bind to the allosteric site and alter the function of the GPCRs. Therefore, this computational approach for ligand binding can be used for prediction of allosteric binders for GPCRs in general.
Supplemental Material
Supplemental material, sj-pdf-1-bbi-10.1177_11779322211037769 for Ligand Docking Methods to Recognize Allosteric Inhibitors for G-Protein-Coupled Receptors by K Harini, S Jayashree, Vikas Tiwari, Sneha Vishwanath and Ramanathan Sowdhamini in Bioinformatics and Biology Insights
Acknowledgments
The authors would like to acknowledge Prof. N. Srinivasan of Indian Institute of Science for useful discussions.
Footnotes
Funding: The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: R.S. would like to acknowledge her JC Bose Fellowship (JC Bose fellowship (SB/S2/JC-071/2015) from Science and Engineering Research Board, Bioinformatics Centre Grant funded by Department of Biotechnology, India (BT/PR40187/BTIS/137/9/2021 and NCBS (TIFR) for infrastructural facilities.
Declaration of Conflicting Interests: The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Author Contributions: K.H., S.J., and V.T. carried out all the work and analyses. R.S. conceptualized the work and both R.S. and S.V. participated in discussions. K.H. and V.T. wrote first draft of the manuscript and S.J., S.V., and R.S. improved it.
ORCID iD: Ramanathan Sowdhamini 
https://orcid.org/0000-0002-6642-2367
Supplemental Material: Supplemental material for this article is available online.