Next Article in Journal
Possible Adverse Effects of High-Dose Nicotinamide: Mechanisms and Safety Assessment
Next Article in Special Issue
New Insights into Key Determinants for Adenosine 1 Receptor Antagonists Selectivity Using Supervised Molecular Dynamics Simulations
Previous Article in Journal
A Fluorescence-Based Method to Measure ADP/ATP Exchange of Recombinant Adenine Nucleotide Translocase in Liposomes
Previous Article in Special Issue
Evolution of Angiotensin Peptides and Peptidomimetics as Angiotensin II Receptor Type 2 (AT2) Receptor Agonists
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Investigations on the Binding Mode of Ligands for the Cannabinoid-Activated G Protein-Coupled Receptor GPR18

1
PharmaCenter Bonn, Pharmaceutical Institute, Pharmaceutical Sciences Bonn (PSB), Pharmaceutical & Medicinal Chemistry, University of Bonn, 53121 Bonn, Germany
2
Research Training Group 1873, University of Bonn, 53127 Bonn, Germany
3
Department of Technology and Biotechnology of Drugs, Faculty of Pharmacy, Jagiellonian University Medical College, Medyczna 9, 30-688 Krakow, Poland
*
Author to whom correspondence should be addressed.
Submission received: 10 March 2020 / Revised: 22 April 2020 / Accepted: 23 April 2020 / Published: 29 April 2020

Abstract

:
GPR18 is an orphan G protein-coupled receptor (GPCR) expressed in cells of the immune system. It is activated by the cannabinoid receptor (CB) agonist ∆9-tetrahydrocannabinol (THC). Several further lipids have been proposed to act as GPR18 agonists, but these results still require unambiguous confirmation. In the present study, we constructed a homology model of the human GPR18 based on an ensemble of three GPCR crystal structures to investigate the binding modes of the agonist THC and the recently reported antagonists which feature an imidazothiazinone core to which a (substituted) phenyl ring is connected via a lipophilic linker. Docking and molecular dynamics simulation studies were performed. As a result, a hydrophobic binding pocket is predicted to accommodate the imidazothiazinone core, while the terminal phenyl ring projects towards an aromatic pocket. Hydrophobic interaction of Cys251 with substituents on the phenyl ring could explain the high potency of the most potent derivatives. Molecular dynamics simulation studies suggest that the binding of imidazothiazinone antagonists stabilizes transmembrane regions TM1, TM6 and TM7 of the receptor through a salt bridge between Asp118 and Lys133. The agonist THC is presumed to bind differently to GPR18 than to the distantly related CB receptors. This study provides insights into the binding mode of GPR18 agonists and antagonists which will facilitate future drug design for this promising potential drug target.

Graphical Abstract

1. Introduction

G protein-coupled receptors (GPCR) represent the largest family of membrane proteins in eukaryotes. They are structurally characterized by seven transmembrane (TM) regions connected by three extracellular (ECL1-3) and three intracellular loops (ICL1-3), an extracellular N-terminal and an intracellular C-terminal domain. Upon binding of the cognate agonist (e.g., biogenic amine neurotransmitter, nucleotide, lipid, amino acid, peptide, glycoprotein) conformational changes are induced. These result in coupling with G proteins, and thereby transducing information from the extracellular to the intracellular compartment and inducing or inhibiting downstream signaling pathways [1,2]. Despite persistent efforts, nearly 100 GPCRs remain orphan, with their endogenous ligands unidentified or unconfirmed [3]. The functionalities and roles of orphan GPCRs under (patho)physiological conditions are in most cases poorly understood. The identification of the endogenous ligands would be helpful for target validation studies and the design of novel therapeutic drugs for orphan GPCRs.
GPR18 is such an orphan GPCR of therapeutic interest, phylogenetically belonging to the δ-branch of class A, rhodopsin-like GPCRs. GPR18 was first described in 1997 and reported to be highly expressed in different tissues and cell lines of the immune system, including spleen, thymus, and leukocytes [4]. The role of GPR18 is still unclear and controversially debated. GPR18 has been proposed by independent groups to be involved in immunological [5,6,7,8] and neurodegenerative processes including Alzheimer’s disease and multiple sclerosis [9,10,11,12,13]. Based on the observation that the activation of GPR18 lowers the intraocular pressure in mice, GPR18 agonists have been proposed for the treatment of glaucoma [14,15]. Antagonists targeting GPR18 may be effective as anticancer drugs [16,17,18], since the receptor was found to be abundantly overexpressed in melanoma metastases and reported to contribute to tumor cell survival through inhibition of apoptosis [17].
In recent years, several studies aimed at the deorphanization of GPR18 have been published. Due to the lack of selective agonists, the moderately potent cannabinoid (CB) receptor agonist ∆9-tetrahydrocannabinol (THC, 1) has been used in pharmacological studies to activate human GPR18, which led to the suggestion to classify GPR18 as a cannabinoid receptor subtype besides CB1 and CB2 [19,20,21,22]. N-Arachidonoylglycine (NAGly, 2) and resolvin D2 (RvD2, 3) were proposed as endogenous agonists of GPR18 [23,24]. However, independent confirmation for both lipids is still lacking, as other groups, including ours, have not been able to confirm their activation of GPR18 [25,26]. We recently described the first GPR18 antagonists based on an imidazothiazinone core structure [21,27]. These were discovered by screening a compound library at the human receptor in a β-arrestin recruitment assay using enzyme complementation technology and THC as an agonist. Based on the screening results, a library of imidazothiazinones was synthesized and their structure–activity relationships (SARs) were investigated. PSB-CB-27 (4) and PSB-CB-5 (5; for structures, see Figure 1) were reported as the first potent and selective GPR18 antagonists [21].
In the present study, we constructed a homology model of the human GPR18 to elucidate the binding mode of the only confirmed agonist so far, the natural product THC, and of selected antagonists by docking and molecular dynamics (MD) simulation studies. Insights into the binding interactions of agonists and antagonists will provide a basis for the rational design of more potent ligands and may eventually contribute to the deorphanization of GPR18.

2. Material and Methods

2.1. Homology Modeling

The crystal structures of the murine μ-opioid receptor in complex with the agonist BU72 (PDB-ID: 5C1M), the human P2Y1 receptor in complex with the allosteric antagonist BPTU (PDB-ID: 4XNV) and the zebrafish lysophosphatidic acid receptor LPA6 in complex with 1-oleoyl-R-glycerol (PDB-ID: 5XSZ) were obtained from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) [28,29,30]. The crystal structures of all three receptors were used as templates for generating homology models of the human GPR18 sequence (accession number: Q14330) retrieved from the UniProt sequence database (http://www.uniprot.org) [31]. The sequences of the murine μ-opioid receptor, the P2Y1R, and the zebrafish lysophosphatidic acid receptor LPA6 were aligned with that of the human GPR18 using Clustal Omega [32]. We generated 500 models for the human GPR18 based on the triple template approach using the standard comparative modeling automodel class implemented in MODELLER (version 9.16, University of California, San Francisco, CA, USA). To ensure correct tertiary protein structure prediction, we introduced a disulfide bridge between Cys94 and Cys172. The best model was selected on the basis of Discrete Optimized Protein Energy (DOPE) scores calculated for the models [33,34]. The generated models were analyzed, and the best models for the human GPR18 were used to perform molecular docking studies, based on the DOPE and GA341 score, PROSA II Z score, and Ramachandran plots. We took into account that the X-ray crystal structure of the lysophosphatidic acid receptor LPA6 is missing part of the ECL2, likely due to low resolution and high flexibility of that region. Nevertheless, we decided to include LPA6 as a template for model generation as it might provide valuable information, e.g., regarding the transmembrane domains and the ligand-binding site. Sequences for the cannabinoid receptors CB1 (P21554) and CB2 (P34972) were retrieved from UniProt.

2.2. Docking Studies

Prior to docking, the homology model of the human GPR18 was prepared using the Protein Preparation Wizard module implemented in Schrödinger [34,35]. In the first step for protein preparation, we preprocessed the structure using the standard protocol at pH 7.4. Docking was performed using Induced Fit Docking (IFD) and Glide as implemented in Schrödinger release 2016 [35,36,37]. In the first step of IFD, Glide ligand docking was performed by removing the side chains of the amino acids in the selected binding pocket. In the second phase of docking, Prime was applied to refine the nearby residues and to optimize the side chains. In the final docking phase, the ligand was re-docked into all induced fit protein structures that were within 30 kcal/mol of the lowest energy structure, by using the Glide XP scoring function. A receptor grid center was specified on the basis of preliminary docking studies, resulting in the highest docking scores for the centroid of Lys174 with a cubic grid side length of 10 Å. Preliminary ensemble docking studies provided highest docking scores and consistent SARs explanation for this selection as well as comparison with published cannabinoid receptor X-ray crystal structures [38].
During the docking simulations, the receptor and the ligands were kept flexible. Following docking, the resulting poses of the best model were selected using the IFD scores and Prime Energy as representative values. The conformations of the docked ligands within an energy window of 2.5 kcal/mol were considered. For Glide docking, the following standard parameters were selected: receptor van der Waals scaling, 0.50, ligand van der Waals scaling, 0.50, and a maximum of 20 poses per ligand. Residues within 5.0 Å of the ligand poses were refined, and the side chains were optimized. The best docking pose was selected based on the IFD score and Prime Energy values.

2.3. Compounds

Synthesis of the compounds which are utilized in this computational study as performed in the Department of Technology and Biotechnology of Drugs Jagiellonian University at Kraków, Poland, and potencies of the compounds were determined at the Department of Pharmaceutical & Medicinal Chemistry, Pharmaceutical Institute, University of Bonn, Germany, as previously reported [21,27]. The synthesis and biological evaluation of the new potent GPR18 antagonist 6 will be published elsewhere.

2.4. Molecular Dynamics Simulation

We selected several successful MD simulations as starting points for our runs [39,40,41,42]. The GPR18 complexes and the unbound GPR18 structure were prepared using the method described above to determine the protein protonation state at pH 7.4. The obtained structures were processed to the CHARMM-GUI molecular simulation program [43,44,45]. The forcefield CHARMM36m was applied for all simulation runs. Ligand parameters were obtained separately from Schrödinger. The orientation of the protein in the phosphatidylcholine lipid bilayer (POPC) was determined using the orientation of proteins in membranes (OPM) database [46]. The cubic water box size was adjusted to the structure size of 20 Å and filled with 0.15 M KCl solution. Water molecules were treated with the transferable intermolecular potential with a 3 points (TIP3P) water model [47]. Equilibration steps for all structures were divided into six steps using NAMD2 [48]. For the first three steps, we selected a runtime of 250 ps in 1 fs intervals. For the last three steps, we selected an equilibration runtime of 2 ns in 2 fs intervals. The system was heated from 0 to 303.15 K during equilibration using the NPT ensemble. During production stages, the system was kept at 303.15 K. Temperature was regulated using the Langevin dynamics thermostat. Production runs were performed for 4 × 50 ns with 4 fs intervals (eventually amounting to 200 ns), and frames were collected every 40 ps using ACEMD by Acellera® with the NVT ensemble [49].

3. Results

So far, no X-ray crystal structure of GPR18 has been published. After performing a BLAST search, three crystal structures with highest sequence identity and overall sequence coverage were chosen as templates: the murine μ-opioid receptor (PDB-ID: 5C1M) in complex with an agonist, the human P2Y1R (PDB-ID: 4XNV) in complex with an allosteric antagonist, and the zebrafish lysophosphatidic acid receptor LPA6 (PDB-ID: 5XSZ) in complex with oleoyl-R-glycerol, showing sequence identities of 24.8%, 25.5% and 27.3%, respectively [28,29,30,50]. Multiple template approaches had been reported to compensate for poor sequence similarity for receptors lacking a template with sequence similarity above 30% [51,52]. Therefore, we decided to include all three templates into the process of homology modeling, although they represent different states of receptor activation. Structures of class A GPCRs belonging to the same δ-branch as GPR18 (P2Y1, LPA6) andone GPCR that is activated by a lipid like GPR18 (LPA6) were selected. BBy this approach, we expected to compensate for gaps and mismatching residues which would be present in a single template approach. The multiple sequence alignment is shown in Figure S1 of Supporting Information.
We subsequently investigated the binding modes of imidazothiazinone antagonists and of the agonist THC in the homology model of the human GPR18. To this end, we used the Induced Fit module implemented in Maestro Schrödinger to propose a binding mode for the selected ligands and to rationalize the potency values obtained in biological studies. Imidazothiazinone derivatives 4 and 5 were selected as representative potent antagonist structures. In addition to the imidazothiazinone core, they both possess a 4-chlorophenoxy substituent connected by a linker which differs in length. For the investigated antagonists 4 and 5, IC50 values of 0.650 and 0.279 µM had been determined [21,27]. In order to investigate whether the proposed antagonist–GPR18 complexes are stable, we performed a 200 ns MD simulation study. Furthermore, we rationalized the SARs of related GPR18 antagonists using a structure-based approach.

3.1. Docking Studies of Antagonists

Recently, several studies on molecular modeling of the human GPR18 have been published [53,54,55,56,57]. However, neither binding mode predictions of THC (or other agonists) nor MD simulations of ≥200 ns of antagonist–GPR18 complexes have been published so far. Kuder et al. reported molecular modeling and docking studies, creating a homology model of the human GPR18 based on the crystal structure of the antagonist-bound human P2Y1 receptor [53]. The imidazothiazinone group of antagonist 5 was predicted to point into a deeper binding pocket towards TM3 where it was hypothesized to form a hydrogen bond with Arg1915.42. The results obtained in the present study indicate a different binding mode which is supported by comprehensive SAR data and based on an ensemble of templates for homology model generation rather than a single, low homology template as in the previous study [53].
The proposed binding mode for antagonist 4 based on the performed docking studies is presented in Figure 2. Antagonist 4 is predicted to bind in the upper third part of the receptor, extending from a hydrophobic cavity formed by ECL2, TM2 and TM3 to an aromatic binding pocket formed by TM6 and TM7, which is a common motif for several GPCRs [58,59,60,61]. Compound 4 likely binds with its imidazothiazinone moiety close to the conserved disulfide bridge of Cys943.25 and Cys172ECL2. Both cysteines and Leu973.28 form a lipophilic binding pocket which is predicted to accommodate the thiazine ring. The keto group of the imidazolone ring likely forms an H-bond with Tyr822.64. Due to the close proximity of Arg782.60, cation–π interactions with the imidazolone system are feasible. The benzylidene ring may extend towards the center of the receptor, where hydrophobic interactions with Thr2727.39 are possible. The hexyloxy linker could bind with several hydrophobic residues (Tyr1604.64, Ile175ECL2, Phe2486.51, Met2757.42) towards an aromatic binding pocket formed by side chains of TM6 and TM7. Additional van der Waals forces for hydrophobic interactions with the benzylidene moiety and the hexyloxy linker may be provided by the alkyl chain of Lys174ECL2. Several aromatic (Phe2486.51, Phe2526.55, Tyr2647.31) and hydrophobic (Cys2516.54 and Leu2556.58) residues of TM6 and TM7 are predicted to form the binding pocket accommodating the 4-chlorophenoxy moiety of compound 5 (see Figure 3).
The smaller antagonist 5 can occupy the same binding cavity as antagonist 4 (see Figure 3). The imidazothiazinone moiety of both compounds can reach the same binding pose. Due to the missing linker, the benzylidine ring is predicted to exhibit an upward shift towards ECL2 where additional cation–π interactions with Lys161ECL2 can be realized. In both cases, the chlorine atoms on the terminal phenyl ring can reach the same binding cavity consisting of aromatic and hydrophobic residues of TM6 and TM7 close to Cys2516.54. Therefore, we expect halogen or methyl substitutions to interact analogously with Cys2516.54. These findings suggest that hydrophobic substituents in position 4 of the terminal phenyl ring of the antagonists are necessary for proper hydrophobic interaction with Cys2516.54 resulting in increased potency.

3.2. MD Simulation Study of Antagonists

Both antagonist–GPR18 complexes were stable during the 200 ns MD simulation runs, which supports our prediction of the binding pocket based on docking studies. The duration of the MD simulation runs was in accordance with similar studies performed for other GPCRs [62,63,64,65]. The behavior of antagonists 4 and 5 in the homology model of GPR18 during the 200 ns MD simulation is presented from a bird’s eye view perspective in Supplementary Information Figures S2 and S3. The 0 ns state refers to the structure of the docked complex after equilibration. The course of the root mean square deviation (RMSD) indicates that the complex of GPR18 with antagonist 5 reached an equilibrated state after approximately 50 ns, and after approximately 100 ns for antagonist 4 (see Figure 4). Compared to the unbound GPR18 structure, the complex of GPR18 with antagonist 5 showed decreased root mean square fluctuation (RMSF) values for TM1, TM2, TM3, TM5 and TM7, and for ECL2 and ECL3, indicating stabilization of these regions upon antagonist binding. Similar results were observed for the complex with the larger antagonist 4, where decreased RMSF values were observed for TM7, ECL2 and ECL3, and ICL2 and ICL3 when compared to the unbound structure. The concept of stabilization of an inactive conformation of the target GPCR upon antagonist binding was postulated for several receptors and supported by mutagenesis experiments, biophysical studies and MD simulations [58,66,67,68]. This had also been observed for the P2Y1 receptor which belongs to the same δ-branch of the class A family of GPCRs. The P2Y1 receptor can be blocked by structurally distinct antagonists that bind to different binding sites, the nucleotide analog MRS2500 and the urea derivative BPTU—both of which stabilize an ionic lock between an aspartic acid residue of ECL2 and an arginine of TM7 [69]. During MD simulations for 2 µs, RSMD values had been significantly lower for the complexes with an antagonist as compared to those with the P2Y1 receptor agonist ADP [69]. A shift in TM3, TM6 and TM7 in the simulation runs with the agonists created a void resulting in receptor activation through a bulk water influx into the binding pocket [69]. Similar observations were reported for several class A family members of GPCRs including µ-opioid receptors and adenosine receptors [70,71,72,73].
To further investigate conformational changes in the receptor, RMSD values for each transmembrane-spanning helix were calculated (see Figure 5). Using the OPM database [46] seven transmembrane region segments were determined: TM1 (Ile231.33–Ser481.58), TM2 (Ile592.41–Phe802.62), TM3 (Glu913.22–Ala1173.53), TM4 (Val1394.43–Tyr1604.64), TM5 (Ala1835.34–Val2095.60), TM6 (Ile2316.34–Phe2546.57) and TM7 (Trp2677.34–Val2897.56). The RMSD values amounted to 2.8, 1.0, 1.2, 1.2, 1.4, 2.0 and 1.9 Å for TM1–TM7, respectively, when comparing the TM regions of the complex of GPR18 and compound 4 at 0 ns and at 200 ns. For the complex with the larger antagonist 4, the RMSD values amounted to 2.2, 2.2, 1.4, 1.9, 2.1, 2.3 and 1.8 Å for TM1–TM7, respectively. The higher RMSD values for antagonist 4 can be explained by the size of the compound when compared to 5: the larger linker requires adaptation of the receptor, resulting in higher RMSD values. In contrast to the behavior of the antagonist-bound complexes, even higher RMSD values were observed for the unbound apo form of GPR18. Here, RMSD values of 4.6, 1.4, 1.8, 2.4, 1.4, 2.9 and 5.4 Å were calculated for TM1–TM7, respectively. Furthermore, the stabilization of TM1, TM6 and TM7 in the presence of an antagonist supports the theory of stabilization of an inactive state of the receptor upon binding of an antagonist [74,75,76].
Potential salt bridges within the receptor were analyzed to further investigate the mode of inhibition. Arg1193.50 of the DRY motif had been proposed to be located in an “arginine cage,” where it forms an ionic lock with Asp1183.49, thus stabilizing the inactive GPR18 [54,77]. Disruption of the ionic lock was postulated to contribute to receptor activation through facilitated movements of TM3 and TM6, resulting in conformational changes towards the intracellular lumen [54]. The authors concluded that stable salt bridges or H-bonds induce a rotamer of Arg1193.50, which is no longer present during receptor activation. The ionic lock between Asp1183.49 and Arg1193.50 was observed in the apo form of GPR18 during our 200 ns MD simulation run, which is consistent with previous studies [54]. Interestingly, we observed differences in the behavior of Arg1193.50 in the apo form as compared to the antagonist-bound complexes: in the apo form, the salt bridge between Arg1193.50 and Asp1183.49 formed after approximately 75 ns and was stable until the end of the MD simulation, while no similar interaction was observed for the antagonist-bound complexes. Asp1183.49 formed a stable salt bridge with Lys133ICL2 in the complexes but not in the apo form. This lysine is neither conserved in the three homology model templates nor in the two CB receptor subtypes. Furthermore, we observed stable ionic interactions of Asp85ECL1 with Lys22N-terminus and of Asp162ECL2 with Lys161ECL2 in the antagonist-bound structures, which were not present in the apo form. Interaction of Glu131ICL2 with Lys1374.41 was observed in all three structures. The salt bridge between Glu2286.31 and Arg2326.35 was stable in the receptor apo form, which was not the case for the antagonist-bound structures. The trajectory for the salt bridge distances is presented in Supplementary Information (Figure S4). We conclude that the binding of an antagonist stabilizes several salt bridges within GPR18, resulting in the stabilization of an inactive conformation of the receptor.
We additionally investigated the binding mode of a new potent antagonist, an analog of 4 and 5, which has a more rigid substituent in position 4 of the phenyl group (a biphenyl derivative). Compound PSB-CB-148 (6) contains a p-cyano-biphenyl group which is larger and at the same time less flexible than the corresponding substituents in antagonists 4 and 5. The imidazothiazinone group is predicted to bind in the same binding cavity as for compounds 4 and 5 (see Figure 6). The trajectory of the linker in the docked structure closely resembles the binding mode of compound 4. Furthermore, the proximity of Arg1915.42 to both oxygen atoms in the linker indicates bidental H-bond interactions. The biphenyl moiety likely binds in a lipophilic binding cavity, where π–π interactions between the phenyl groups and the aromatic residues Phe2486.51, Phe2526.55 and Tyr2647.31 are feasible. Interactions of the terminal phenyl group with Cys2516.54 are not observed for 6. Due to its decreased flexibility, the terminal group does not allow this interaction. The shift in the phenyl group is predicted to place the cyano moiety in close proximity to Asn1855.39. Upon inspection of Asn1855.39, several rotamers were found which could form H-bonds with the nitrile (see Figure S5 in Supplementary Information).
The obtained data of the docking studies were used to re-analyze the SARs of previously published antagonists [21]. The summarized results are presented in Figure 7 (for structures, see Figures S6 and S7 in Supporting Information). The linker size was found to have an impact on the potency of the tested antagonists. The antagonist containing a hexyloxy linker (4) showed an almost 10-fold increase in potency compared to the analog with the shorter ethyloxy linker (7) (IC50 of 0.650 µM versus 5.00 µM). Prolongation of the ethyloxy linker resulted in increased inhibitory potency, with hexyloxy being optimal (IC50 = 0.650 µM), while larger linkers, i.e., heptyloxy (11) and octyloxy (12), led to slightly less potent antagonists (IC50 = 1.71 and 1.15 µM). Our docking results suggest that the hexyloxy linker is required for the 4-chlorophenoxy moiety to reach the aromatic binding pocket and to form hydrophobic interactions with Cys2516.54. The shorter alkyloxy linker is less well stabilized in the hydrophobic cavity formed by Tyr1604.64, Ile175ECL2, Phe2486.51 and Met2757.42. The decrease in potency observed for compounds 11 and 12 despite their higher lipophilicity could be explained by limited space in the binding cavity or unfavorable adaptation of the alkyloxy linker, resulting in a shifted binding position for the 4-chlorophenoxy moiety which prohibits optimal interaction with Cys2516.54. Among the smaller compounds missing an additional linker between the benzylidene and the substituted phenoxy ring, the most potent antagonists contained a hydrophobic substituent in position 4 of the phenyl ring (compounds 5, 1416). Hydrophobic interactions of substituents in position 4 of the phenoxy residue with Cys2516.54 are supported by acceptance of both chlorine and methyl groups in compounds 4 and 13, resulting in comparable IC50 values (0.650 and 0.238 µM). The potency (IC50 values) of the compounds decreased in the following rank order Cl (0.279 µM) > Br (1.73 µM) ≥ CH3 (3.59 µM) > F (> 10 µM), indicating that the size and lipophilicity of the substituent plays a major role. Decreased potency observed for antagonists containing larger substituents in position 4 such as ethyl (17) or isopropyl (18) can be explained by the limited space of the binding pocket in proximity to Cys2516.54. Moreover, the substitution position on the phenyl ring proved to have an effect on the potency of the compounds. Antagonist 19 (o,o-dimethyl-substituted), for example, was inactive (IC50 > 10 µM). Antagonists containing different heterocycles in place of the imidazothiazinone moiety (2034) showed lower potency as compared to antagonist 5. In our homology model, two aromatic residues close to the hydrophobic binding pocket, Tyr812.63 and Trp87ECL1, may form π–π interactions with antagonists containing an additional aromatic group attached to the heterocycle (see Figure S8 in Supporting Information). The ethylthio linker connecting the imidazolone ring with the phenyl ring in compound 32 might be beneficial to enable proper binding for π–π interactions. The results suggest that the imidazothiazinone heterocycle is optimal to allow hydrophobic packing in the binding pocket close to the disulfide bridge of ECL2.
In conclusion, the docking studies, MD simulations and SARs of imidazothiazinones as well as antagonists containing smaller heterocycles further support our suggested binding mode of an aromatic and lipophilic binding pocket of the human GPR18 for antagonists. The most potent antagonists of this series likely interact with Cys2516.54 through lipophilic interactions, and this additional interaction is predicted to be the reason for their high potency.

3.3. Binding Mode of THC

As a next step, we explored the most likely binding pocket for the GPR18 agonist THC (1). The ability of the potent CB receptor agonist THC to activate GPR18 with moderate potency had led to the suggestion to classify GPR18 as a novel CB receptor subtype [19]. Lipophilicity is a feature shared by GPR18 agonists and antagonists [78,79]. THC is regarded as a promiscuous ligand acting not only at cannabinoid but also at several non-cannabinoid receptors [80,81,82,83,84,85]. Studies on the binding mode of cannabinoids at the cannabinoid receptors CB1 and CB2 proposed a binding portal between TM6 and TM7 from the lipid-facing side of the receptor for the entrance of agonists [86,87,88]. Such entry is unique among GPCRs, as ligands typically reach the binding pocket between TM3 and TM7 from the extracellular lumen.
To date, two crystal structures of the CB1 receptor bound to THC-related compounds are available (PDB-ID: 5XR8, 5XRA) [38]. As observed for many other GPCRs, the agonist binding site, which is very lipophilic in the case of the CB1 receptor, is located between a highly conserved Trp6.48 and ECL2 [89,90]. The tricyclic THC ring system is stabilized through lipophilic as well as π–π interactions with an aromatic cluster (Phe1702.57, Phe1742.61, Phe1772.64, Phe1893.25, Phe268ECL2, Phe3797.35). Several previous mutagenesis studies have confirmed the key role of the aromatic residues for the binding of cannabinoids [38,91,92,93]. The alkyl chain of the agonists extends towards a binding cleft formed by several lipophilic residues (Leu1933.29, Val1963.32, Tyr2755.39, Leu2765.40, Leu3596.51 and Met3636.55).
Given the low sequence similarity between the cannabinoid receptors CB1, CB2, and GPR18 (18.7 and 23.7%, respectively), similar binding of the THC ring system in GPR18 cannot be taken for granted. Amino acid residues Val1963.32, Phe268ECL2, Tyr2755.39, Met3636.55 and Phe3797.35 are conserved in both CB receptor subtypes but replaced in GPR18 by leucine, serine, arginine, phenylalanine and glycine, respectively (see Figure S9 in Supplementary Information for multiple sequence alignment). Phe1742.61 and Leu1933.29, but not Leu3596.51, are conserved among all three receptors. The absence of the aromatic network responsible for the binding of the THC ring system in the CB receptors suggests a different binding mode for the agonist THC at GPR18.
Docking studies of THC were performed using the generated homology model of the human GPR18. We observed that THC appears to bind closer to TM4 and TM5 as compared to the cannabinoids in the X-ray crystal structures of the CB1 receptor (see Figure 8). The phenyl group of the tricyclic THC ring system is predicted to bind in a cleft formed by several lipophilic (Val1023.33, Ile175ECL2, Phe2486.51, Phe2526.55) and hydrophilic (Lys161ECL2, Lys174ECL2, Asn1885.39, Arg1915.42, His2496.52) amino acid residues. H-bond interactions are feasible for the oxygen atoms of the chromene moiety and Lys161ECL2, as well as the hydroxy group and Asn1885.39 and Arg1915.42. The cyclohexenyl moiety is likely accommodated in a lipophilic binding pocket formed by Thr1524.56, Pro1554.59, Leu1564.60, Val1845.35 and the alkyl side chain of Arg1915.42. The alkyl group of the agonist likely projects towards TM7, where it can be stabilized through lipophilic interactions with Phe2486.51, Phe2526.55 and Met2757.42. The binding modes of THC in the CB1 receptor as compared to GPR18 are shown in Figure S10 of Supporting Information.
We propose that the tricyclic THC ring system binds in a binding cavity of GPR18 distant to the orthosteric binding site of the CB1 receptor. The absence of aromatic residues in ECL2 of GPR18 may contribute to the proposed shifted binding mode of THC, as π–π stacking with a phenylalanine in position 2.57 is not possible. However, the binding cleft for the alkyl chain is predicted to be overlapping in both receptors. It should be pointed out that THC displays much higher potency at CB1 (and CB2) receptors as compared to GPR18.
Our results suggest that THC shares a common binding pocket with the imidazothiazinone antagonists (see Figure 9). While the imidazothiazinone moiety of the antagonists is predicted to bind in a lipophilic pocket formed by amino acid residues of TM2 and TM7, the benzylidene group is suggested to project towards the putative binding site of the chromene and alkyl group of THC. This is supported by experimental data showing that imidazothiazinone antagonists containing lipophilic residues act as competitive antagonists versus THC [21].

4. Conclusions

Since only approximately 10% of the non-olfactory GPCRs are covered by structural studies, meaningful prediction of ligand-binding modes represents one of the greatest challenges in molecular modeling [94]. In particular, homology modeling assessment of receptors with no resolved closely related crystal structures requires further experimental validation. In the present study, we generated a homology model of the orphan GPR18 and predicted the binding modes of the confirmed agonist THC as well as the most potent class of antagonists containing an imidazothiazinone scaffold. Despite the lack of closely related X-ray crystal structures, we successfully performed docking and MD simulation studies of antagonist complexes which were in agreement with the extensive published SAR data. The investigated potent antagonists are predicted to share the same binding site for the imidazothiazinone core. The linker of the antagonists is likely accommodated in a lipophilic binding cleft shared by the alkyl chain of the agonist THC. The 200 ns MD simulation runs suggested stabilization of a receptor conformation by antagonists which was not observed for the unbound receptor structure. Stabilization of a salt bridge between Asp1183.49 and Lys133ICL2 through imidothiazinone-based antagonists may play a role in the inhibition mechanism. Our docking studies suggest a different binding mode of the agonist THC in GPR18 as compared to that observed in cannabinoid receptors. However, future structural studies will be required to confirm the proposed interactions. The presented data provide a well-founded hypothesis that will support the rational design of new ligands for this poorly investigated receptor which has potential as a future drug target.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2218-273X/10/5/686/s1, Figure S1: Multiple sequence alignment of GPR18 for homology modeling, Figure S2: Time scale of the molecular dynamics simulation of antagonist 4, Figure S3: Time scale of the molecular dynamics simulation of antagonist 5, Figure S4: Trajectories of salt bridges, Figure S5: Possible interaction of antagonist 6 with Asn185, Figure S6: Structures of GPR18 imidazothiazinone antagonists, Figure S7: Structures of GPR18 antagonists with modification of the core structure, Figure S8: Comparison of the binding modes of antagonists 5 and 32, Figure S9: Multiple sequence alignment of human GPR18 and the cannabinoid receptors CB1 and CB2, Figure S10: Comparison of the binding mode of THC to GPR18 with the binding of THC derivatives to the CB1 receptor, Figure S11: Multiple Sequence alignment of top 10 scoring BLAST sequences, Figure S12: Residues considered to be important for the binding of antagonists, Table S1: Results for 10 scoring BLAST sequences for GPR18.

Author Contributions

Conceptualization, A.N. and C.E.M.; methodology, A.N., V.E., A.B.M., C.T.S., V.N. and K.K.-K.; first draft preparation, A.N.; writing of manuscript, A.N. and C.E.M. with contributions of all co-authors; review and editing, all co-authors; supervision, C.E.M. and V.N.; project administration and funding acquisition, C.E.M. All authors have read and agreed to the published version of the manuscript.

Funding

A.N., A.B.M., C.T.S., and C.E.M. were funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 214362475 / GRK1873/2 and the BMBF-funded Bonn International Graduate School Drug Sciences (BIGS DrugS). C.T.S. received a BAYER PhD fellowship through BIGS DrugS. A.B.M. was funded by the Ministry of Finance of Indonesia in the scheme of Indonesia Endowment Fund for Education (Lembaga Pengelola Dana Pendidikan (LPDP).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hilger, D.; Masureel, M.; Kobilka, B.K. Structure and dynamics of GPCR signaling complexes. Nat. Struct. Mol. Biol. 2018, 25, 4–12. [Google Scholar] [CrossRef] [PubMed]
  2. Kobilka, B.K. G protein coupled receptor structure and activation. Biochim. Biophys. Acta 2007, 1768, 794–807. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Tang, X.-L.; Wang, Y.; Li, D.-L.; Luo, J.; Liu, M.-Y. Orphan G protein-coupled receptors (GPCRs): Biological functions and potential drug targets. Acta Pharmacol. Sin. 2012, 33, 363–371. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Gantz, I.; Muraoka, A.; Yang, Y.K.; Samuelson, L.C.; Zimmerman, E.M.; Cook, H.; Yamada, T. Cloning and chromosomal localization of a gene (GPR18) encoding a novel seven transmembrane receptor highly expressed in spleen and testis. Genomics 1997, 42, 462–466. [Google Scholar] [CrossRef] [PubMed]
  5. Sumida, H.; Cyster, J.G. G-Protein Coupled Receptor 18 Contributes to Establishment of the CD8 Effector T Cell Compartment. Front. Immunol. 2018, 9, 660. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Wang, X.; Sumida, H.; Cyster, J.G. GPR18 is required for a normal CD8αα intestinal intraepithelial lymphocyte compartment. J. Exp. Med. 2014, 211, 2351–2359. [Google Scholar] [CrossRef] [Green Version]
  7. Becker, A.M.; Callahan, D.J.; Richner, J.M.; Choi, J.; DiPersio, J.F.; Diamond, M.S.; Bhattacharya, D. GPR18 controls reconstitution of mouse small intestine intraepithelial lymphocytes following bone marrow transplantation. PLoS ONE 2015, 10, e0133854. [Google Scholar] [CrossRef]
  8. Pridgeon, J.W.; Klesius, P.H. G-protein coupled receptor 18 (GPR18) in channel catfish: Expression analysis and efficacy as immunostimulant against Aeromonas hydrophila infection. Fish Shellfish Immunol. 2013, 35, 1070–1078. [Google Scholar] [CrossRef]
  9. Morales, P.; Reggio, P.H. An update on non-CB1, non-CB2 cannabinoid related G-protein-coupled receptors. Cannabis Cannabinoid Res. 2017, 2, 265–273. [Google Scholar] [CrossRef] [Green Version]
  10. Reyes-Resina, I.; Navarro, G.; Aguinaga, D.; Canela, E.I.; Schoeder, C.T.; Załuski, M.; Kieć-Kononowicz, K.; Saura, C.A.; Müller, C.E.; Franco, R. Molecular and functional interaction between GPR18 and cannabinoid CB2 G-protein-coupled receptors. Relevance in neurodegenerative diseases. Biochem. Pharmacol. 2018, 157, 169–179. [Google Scholar] [CrossRef] [Green Version]
  11. McHugh, D. GPR18 in microglia: Implications for the CNS and endocannabinoid system signalling. Br. J. Pharmacol. 2012, 167, 1575–1582. [Google Scholar] [CrossRef] [Green Version]
  12. Walter, L.; Franklin, A.; Witting, A.; Wade, C.; Xie, Y.; Kunos, G.; Mackie, K.; Stella, N. Nonpsychotropic cannabinoid receptors regulate microglial cell migration. J. Neurosci. 2003, 23, 1398–1405. [Google Scholar] [CrossRef] [Green Version]
  13. Haugh, O.; Penman, J.; Irving, A.J.; Campbell, V.A. The emerging role of the cannabinoid receptor family in peripheral and neuro-immune interactions. Curr. Drug Targets 2016, 17, 1834–1840. [Google Scholar] [CrossRef]
  14. Miller, S.; Leishman, E.; Oehler, O.; Daily, L.; Murataeva, N.; Wager-Miller, J.; Bradshaw, H.; Straiker, A. Evidence for a GPR18 role in diurnal regulation of intraocular pressure. Investig. Ophthalmol. Vis. Sci. 2016, 57, 6419–6426. [Google Scholar] [CrossRef] [Green Version]
  15. Caldwell, M.D.; Hu, S.S.-J.; Viswanathan, S.; Bradshaw, H.; Kelly, M.E.M.; Straiker, A. A GPR18-based signalling system regulates IOP in murine eye. Br. J. Pharmacol 2013, 169, 834–843. [Google Scholar] [CrossRef] [Green Version]
  16. Haskó, J.; Fazakas, C.; Molnár, J.; Nyúl-Tóth, Á.; Herman, H.; Hermenean, A.; Wilhelm, I.; Persidsky, Y.; Krizbai, I.A. CB2 receptor activation inhibits melanoma cell transmigration through the blood-brain barrier. Int. J. Mol. Sci. 2014, 15, 8063–8074. [Google Scholar] [CrossRef] [Green Version]
  17. Qin, Y.; Verdegaal, E.M.E.; Siderius, M.; Bebelman, J.P.; Smit, M.J.; Leurs, R.; Willemze, R.; Tensen, C.P.; Osanto, S. Quantitative expression profiling of G-protein-coupled receptors (GPCRs) in metastatic melanoma: The constitutively active orphan GPCR GPR18 as novel drug target. Pigment Cell Melanoma Res. 2011, 24, 207–218. [Google Scholar] [CrossRef]
  18. Noreen, N.; Muhammad, F.; Akhtar, B.; Azam, F.; Anwar, M.I. Is cannabidiol a promising substance for new drug development? A review of its potential therapeutic applications. Crit. Rev. Eukaryot. Gene Expr. 2018, 28, 73–86. [Google Scholar] [CrossRef]
  19. Console-Bram, L.; Brailoiu, E.; Brailoiu, G.C.; Sharir, H.; Abood, M.E. Activation of GPR18 by cannabinoid compounds: A tale of biased agonism. Br. J. Pharmacol. 2014, 171, 3908–3917. [Google Scholar] [CrossRef] [Green Version]
  20. Pertwee, R.G.; Howlett, A.C.; Abood, M.E.; Alexander, S.P.H.; Di Marzo, V.; Elphick, M.R.; Greasley, P.J.; Hansen, H.S.; Kunos, G.; Mackie, K.; et al. International Union of Basic and Clinical Pharmacology. LXXIX. Cannabinoid receptors and their ligands: Beyond CB1 and CB2. Pharmacol. Rev. 2010, 62, 588–631. [Google Scholar] [CrossRef] [Green Version]
  21. Schoeder, C.T.; Kaleta, M.; Mahardhika, A.B.; Olejarz-Maciej, A.; Łażewska, D.; Kieć-Kononowicz, K.; Müller, C.E. Structure-activity relationships of imidazothiazinones and analogs as antagonists of the cannabinoid-activated orphan G protein-coupled receptor GPR18. Eur. J. Med. Chem. 2018, 155, 381–397. [Google Scholar] [CrossRef] [PubMed]
  22. Schoeder, C.T.; Hess, C.; Madea, B.; Meiler, J.; Müller, C.E. Pharmacological evaluation of new constituents of “Spice”: Synthetic cannabinoids based on indole, indazole, benzimidazole and carbazole scaffolds. Forensic Toxicol. 2018, 36, 385–403. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Kohno, M.; Hasegawa, H.; Inoue, A.; Muraoka, M.; Miyazaki, T.; Oka, K.; Yasukawa, M. Identification of N-arachidonylglycine as the endogenous ligand for orphan G-protein-coupled receptor GPR18. Biochem. Biophys. Res. Commun. 2006, 347, 827–832. [Google Scholar] [CrossRef]
  24. Chiang, N.; Dalli, J.; Colas, R.A.; Serhan, C.N. Identification of resolvin D2 receptor mediating resolution of infections and organ protection. J. Exp. Med. 2015, 212, 1203–1217. [Google Scholar] [CrossRef] [Green Version]
  25. Yin, H.; Chu, A.; Li, W.; Wang, B.; Shelton, F.; Otero, F.; Nguyen, D.G.; Caldwell, J.S.; Chen, Y.A. Lipid G protein-coupled receptor ligand identification using beta-arrestin PathHunter assay. J. Biol. Chem. 2009, 284, 12328–12338. [Google Scholar] [CrossRef] [Green Version]
  26. van Lu, B.; Puhl, H.L.; Ikeda, S.R. N-Arachidonyl glycine does not activate G protein-coupled receptor 18 signaling via canonical pathways. Mol. Pharmacol. 2013, 83, 267–282. [Google Scholar]
  27. Rempel, V.; Atzler, K.; Behrenswerth, A.; Karcz, T.; Schoeder, C.; Hinz, S.; Kaleta, M.; Thimm, D.; Kieć-Kononowicz, K.; Müller, C.E. Bicyclic imidazole-4-one derivatives: A new class of antagonists for the orphan G protein-coupled receptors GPR18 and GPR55. Med. Chem. Commun. 2014, 5, 632–649. [Google Scholar] [CrossRef]
  28. Huang, W.; Manglik, A.; Venkatakrishnan, A.J.; Laeremans, T.; Feinberg, E.N.; Sanborn, A.L.; Kato, H.E.; Livingston, K.E.; Thorsen, T.S.; Kling, R.C.; et al. Structural insights into µ-opioid receptor activation. Nature 2015, 524, 315–321. [Google Scholar] [CrossRef] [Green Version]
  29. Taniguchi, R.; Inoue, A.; Sayama, M.; Uwamizu, A.; Yamashita, K.; Hirata, K.; Yoshida, M.; Tanaka, Y.; Kato, H.E.; Nakada-Nakura, Y.; et al. Structural insights into ligand recognition by the lysophosphatidic acid receptor LPA6. Nature 2017, 548, 356–360. [Google Scholar] [CrossRef]
  30. Zhang, D.; Gao, Z.-G.; Zhang, K.; Kiselev, E.; Crane, S.; Wang, J.; Paoletta, S.; Yi, C.; Ma, L.; Zhang, W.; et al. Two disparate ligand-binding sites in the human P2Y1 receptor. Nature 2015, 520, 317–321. [Google Scholar] [CrossRef]
  31. UniProt Consortium. UniProt: A hub for protein information. Nucleic Acids Res. 2015, 43, D204–D212. [Google Scholar] [CrossRef] [PubMed]
  32. Sievers, F.; Wilm, A.; Dineen, D.; Gibson, T.J.; Karplus, K.; Li, W.; Lopez, R.; McWilliam, H.; Remmert, M.; Söding, J.; et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol. Syst. Biol. 2011, 7, 539. [Google Scholar] [CrossRef] [PubMed]
  33. Sali, A.; Blundell, T.L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993, 234, 779–815. [Google Scholar] [CrossRef]
  34. Webb, B.; Sali, A. Protein structure modeling with MODELLER. Methods Mol. Biol. 2014, 1137, 1–15. [Google Scholar]
  35. Friesner, R.A.; Banks, J.L.; Murphy, R.B.; Halgren, T.A.; Klicic, J.J.; Mainz, D.T.; Repasky, M.P.; Knoll, E.H.; Shelley, M.; Perry, J.K.; et al. Glide: A new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J. Med. Chem. 2004, 47, 1739–1749. [Google Scholar] [CrossRef]
  36. Halgren, T.A.; Murphy, R.B.; Friesner, R.A.; Beard, H.S.; Frye, L.L.; Pollard, W.T.; Banks, J.L. Glide: A new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening. J. Med. Chem. 2004, 47, 1750–1759. [Google Scholar] [CrossRef]
  37. Sherman, W.; Day, T.; Jacobson, M.P.; Friesner, R.A.; Farid, R. Novel procedure for modeling ligand/receptor induced fit effects. J. Med. Chem. 2006, 49, 534–553. [Google Scholar] [CrossRef]
  38. Hua, T.; Vemuri, K.; Nikas, S.P.; Laprairie, R.B.; Wu, Y.; Qu, L.; Pu, M.; Korde, A.; Jiang, S.; Ho, J.-H.; et al. Crystal structures of agonist-bound human cannabinoid receptor CB1. Nature 2017, 547, 468–471. [Google Scholar] [CrossRef]
  39. Abdelrahman, A.; Yerande, S.G.; Namasivayam, V.; Klapschinski, T.A.; Alnouri, M.W.; El-Tayeb, A.; Müller, C.E. Substituted 4-phenylthiazoles: Development of potent and selective A1, A3 and dual A1/A3 adenosine receptor antagonists. Eur. J. Med. Chem. 2020, 186, 111879. [Google Scholar] [CrossRef]
  40. Ciancetta, A.; O’Connor, R.D.; Paoletta, S.; Jacobson, K.A. Demystifying P2Y1 Receptor Ligand Recognition through Docking and Molecular Dynamics Analyses. J. Chem. Inf. Model. 2017, 57, 3104–3123. [Google Scholar] [CrossRef]
  41. Liang, Y.-L.; Belousoff, M.J.; Fletcher, M.M.; Zhang, X.; Khoshouei, M.; Deganutti, G.; Koole, C.; Furness, S.G.B.; Miller, L.J.; Hay, D.L.; et al. Structure and Dynamics of Adrenomedullin Receptors AM1 and AM2 Reveal Key Mechanisms in the Control of Receptor Phenotype by Receptor Activity-Modifying Proteins. ACS Pharmacol. Transl. Sci. 2020, 3, 263–284. [Google Scholar] [CrossRef]
  42. Tosh, D.K.; Rao, H.; Bitant, A.; Salmaso, V.; Mannes, P.; Lieberman, D.I.; Vaughan, K.L.; Mattison, J.A.; Rothwell, A.C.; Auchampach, J.A.; et al. Design and in Vivo Characterization of A1 Adenosine Receptor Agonists in the Native Ribose and Conformationally Constrained (N)-Methanocarba Series. J. Med. Chem. 2019, 62, 1502–1522. [Google Scholar] [CrossRef] [PubMed]
  43. Brooks, B.R.; Brooks, C.L.; Mackerell, A.D.; Nilsson, L.; Petrella, R.J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545–1614. [Google Scholar] [CrossRef] [PubMed]
  44. Jo, S.; Kim, T.; Iyer, V.G.; Im, W. CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859–1865. [Google Scholar] [CrossRef] [PubMed]
  45. Lee, J.; Cheng, X.; Swails, J.M.; Yeom, M.S.; Eastman, P.K.; Lemkul, J.A.; Wei, S.; Buckner, J.; Jeong, J.C.; Qi, Y.; et al. CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field. J. Chem. Theory Comput. 2016, 12, 405–413. [Google Scholar] [CrossRef] [PubMed]
  46. Lomize, M.A.; Pogozheva, I.D.; Joo, H.; Mosberg, H.I.; Lomize, A.L. OPM database and PPM web server: Resources for positioning of proteins in membranes. Nucleic Acids Res. 2012, 40, D370–D376. [Google Scholar] [CrossRef] [PubMed]
  47. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. [Google Scholar] [CrossRef]
  48. Phillips, J.C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R.D.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. [Google Scholar] [CrossRef] [Green Version]
  49. Harvey, M.J.; Giupponi, G.; Fabritiis, G.D. ACEMD: Accelerating Biomolecular Dynamics in the Microsecond Time Scale. J. Chem. Theory Comput. 2009, 5, 1632–1639. [Google Scholar] [CrossRef] [Green Version]
  50. Altschul, S.F.; Gish, W.; Miller, W.; Myers, E.W.; Lipman, D.J. Basic local alignment search tool. J. Mol. Biol. 1990, 215, 403–410. [Google Scholar] [CrossRef]
  51. Chaudhari, R.; Heim, A.J.; Li, Z. Improving homology modeling of G-protein coupled receptors through multiple-template derived conserved inter-residue interactions. J. Comput. Aided Mol. Des. 2015, 29, 413–420. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Larsson, P.; Wallner, B.; Lindahl, E.; Elofsson, A. Using multiple templates to improve quality of homology models in automated homology modeling. Protein Sci. 2008, 17, 990–1002. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Kuder, K.J.; Karcz, T.; Kaleta, M.; Kieć-Kononowicz, K. Molecular Modeling of an Orphan GPR18 Receptor. Lett. Drug Des. Discov. 2019, 16, 1167–1174. [Google Scholar] [CrossRef]
  54. Sotudeh, N.; Morales, P.; Hurst, D.P.; Lynch, D.L.; Reggio, P.H. Towards a molecular understanding of the cannabinoid related orphan receptor GPR18: A focus on its constitutive activity. Int. J. Mol. Sci. 2019, 20, 2300. [Google Scholar] [CrossRef] [Green Version]
  55. Schmeisser, M.G.; Pearsall, E.A.; Reggio, P.H. Construction of a GPR18 receptor model using conformational memories. Biophys. J. 2013, 104, 409a. [Google Scholar] [CrossRef] [Green Version]
  56. Reynolds, J.L. Inhibition of GPR18 through docking of known antagonists using a homology model. Biophys. J. 2015, 108, 513a. [Google Scholar] [CrossRef] [Green Version]
  57. Kothandan, G.; Cho, S.J. Homology modeling of GPR18 Receptor, an orphan G-protein-coupled receptor. J. Chosun Nat. Sci. 2013, 6, 16–20. [Google Scholar] [CrossRef] [Green Version]
  58. Rosenbaum, D.M.; Rasmussen, S.G.F.; Kobilka, B.K. The structure and function of G-protein-coupled receptors. Nature 2009, 459, 356–363. [Google Scholar] [CrossRef] [Green Version]
  59. Weis, W.I.; Kobilka, B.K. The molecular basis of G protein-coupled receptor activation. Annu. Rev. Biochem. 2018, 87, 897–919. [Google Scholar] [CrossRef]
  60. Gacasan, S.B.; Baker, D.L.; Parrill, A.L. G protein-coupled receptors: The evolution of structural insight. AIMS Biophys. 2017, 4, 491–527. [Google Scholar] [CrossRef]
  61. Zhang, D.; Zhao, Q.; Wu, B. Structural studies of G protein-coupled receptors. Mol. Cells 2015, 38, 836–842. [Google Scholar] [PubMed] [Green Version]
  62. Grossfield, A. Recent progress in the study of G protein-coupled receptors with molecular dynamics computer simulations. Biochim. Biophys. Acta 2011, 1808, 1868–1878. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Chen, H.; Fu, W.; Wang, Z.; Wang, X.; Lei, T.; Zhu, F.; Li, D.; Chang, S.; Xu, L.; Hou, T. Reliability of docking-based virtual screening for GPCR ligands with homology modeled structures: A case study of the angiotensin II type I receptor. ACS Chem. Neurosci. 2019, 10, 677–689. [Google Scholar] [CrossRef]
  64. Ciancetta, A.; Rubio, P.; Lieberman, D.I.; Jacobson, K.A. A3 adenosine receptor activation mechanisms: Molecular dynamics analysis of inactive, active, and fully active states. J. Comput. Aided Mol. Des. 2019, 33, 983–996. [Google Scholar] [CrossRef] [PubMed]
  65. Kashani-Amin, E.; Sakhteman, A.; Larijani, B.; Ebrahim-Habibi, A. Presence of carbohydrate binding modules in extracellular region of class C G-protein coupled receptors (C GPCR): An in silico investigation on sweet taste receptor. J. Biosci. 2019, 44, 138. [Google Scholar] [CrossRef]
  66. Wacker, D.; Stevens, R.C.; Roth, B.L. How ligands illuminate GPCR molecular pharmacology. Cell 2017, 170, 414–427. [Google Scholar] [CrossRef] [Green Version]
  67. Heydenreich, F.M.; Vuckovic, Z.; Matkovic, M.; Veprintsev, D.B. Stabilization of G protein-coupled receptors by point mutations. Front. Pharmacol. 2015, 6, 82. [Google Scholar] [CrossRef] [Green Version]
  68. Deganutti, G.; Moro, S.; Reynolds, C.A. Peeking at G-protein-coupled receptors through the molecular dynamics keyhole. Future Med. Chem. 2019, 11, 599–615. [Google Scholar] [CrossRef] [Green Version]
  69. Yuan, S.; Chan, H.C.S.; Vogel, H.; Filipek, S.; Stevens, R.C.; Palczewski, K. The molecular mechanism of P2Y1 receptor activation. Angew. Chem. 2016, 55, 10331–10335. [Google Scholar] [CrossRef]
  70. Ribeiro, J.M.L.; Filizola, M. Insights from molecular dynamics simulations of a number of G-Protein coupled receptor targets for the treatment of pain and opioid use disorders. Front. Mol. Neurosci. 2019, 12, 207. [Google Scholar] [CrossRef] [Green Version]
  71. An, X.; Bai, Q.; Bing, Z.; Zhou, S.; Shi, D.; Liu, H.; Yao, X. How does agonist and antagonist binding lead to different conformational ensemble equilibria of the κ-opioid receptor: Insight from long-time gaussian accelerated molecular dynamics simulation. ACS Chem. Neurosci. 2019, 10, 1575–1584. [Google Scholar] [CrossRef] [PubMed]
  72. Tautermann, C.S.; Seeliger, D.; Kriegl, J.M. What can we learn from molecular dynamics simulations for GPCR drug design? Comput. Struct. Biotechnol. J. 2015, 13, 111–121. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Renault, P.; Louet, M.; Marie, J.; Labesse, G.; Floquet, N. Molecular dynamics simulations of the allosteric modulation of the adenosine A2a receptor by a mini-G Protein. Sci. Rep. 2019, 9, 1–12. [Google Scholar] [CrossRef] [PubMed]
  74. Bartuzi, D.; Kaczor, A.A.; Matosiuk, D. Molecular mechanisms of allosteric probe dependence in μ opioid receptor. J. Biomol. Struct. Dyn. 2019, 37, 36–47. [Google Scholar] [CrossRef]
  75. Dalton, J.A.R.; Lans, I.; Giraldo, J. Quantifying conformational changes in GPCRs: Glimpse of a common functional mechanism. BMC Bioinform. 2015, 16, 124. [Google Scholar] [CrossRef] [Green Version]
  76. Lee, Y.; Basith, S.; Choi, S. Recent advances in structure-based drug design targeting class A G protein-coupled receptors utilizing crystal structures and computational simulations. J. Med. Chem. 2018, 61, 1–46. [Google Scholar] [CrossRef]
  77. Ballesteros, J.; Kitanovic, S.; Guarnieri, F.; Davies, P.; Fromme, B.J.; Konvicka, K.; Chi, L.; Millar, R.P.; Davidson, J.S.; Weinstein, H.; et al. Functional microdomains in G-protein-coupled receptors. The conserved arginine-cage motif in the gonadotropin-releasing hormone receptor. J. Biol. Chem. 1998, 273, 10445–10453. [Google Scholar] [CrossRef] [Green Version]
  78. Wiley, J.L.; Marusich, J.A.; Thomas, B.F. Combination chemistry: Structure-activity relationships of novel psychoactive cannabinoids. Curr. Top. Behav. Neurosci. 2017, 32, 231–248. [Google Scholar]
  79. Morales, P.; Hurst, D.P.; Reggio, P.H. Molecular targets of the phytocannabinoids-a complex picture. Prog. Chem. Org. Nat. Prod. 2017, 103, 103–131. [Google Scholar]
  80. Ryberg, E.; Larsson, N.; Sjögren, S.; Hjorth, S.; Hermansson, N.-O.; Leonova, J.; Elebring, T.; Nilsson, K.; Drmota, T.; Greasley, P.J. The orphan receptor GPR55 is a novel cannabinoid receptor. Br. J. Pharmacol. 2007, 152, 1092–1101. [Google Scholar] [CrossRef]
  81. de Petrocellis, L.; Di Marzo, V. Non-CB1, non-CB2 receptors for endocannabinoids, plant cannabinoids, and synthetic cannabimimetics: Focus on G-protein-coupled receptors and transient receptor potential channels. J. Neuroimmune Pharmacol. 2010, 5, 103–121. [Google Scholar] [CrossRef] [PubMed]
  82. Qin, N.; Neeper, M.P.; Liu, Y.; Hutchinson, T.L.; Lubin, M.L.; Flores, C.M. TRPV2 is activated by cannabidiol and mediates CGRP release in cultured rat dorsal root ganglion neurons. J. Neurosci. 2008, 28, 6231–6238. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. de Petrocellis, L.; Vellani, V.; Schiano-Moriello, A.; Marini, P.; Magherini, P.C.; Orlando, P.; Di Marzo, V. Plant-derived cannabinoids modulate the activity of transient receptor potential channels of ankyrin type-1 and melastatin type-8. J. Pharmacol. Exp. Ther. 2008, 325, 1007–1015. [Google Scholar] [CrossRef] [PubMed]
  84. Xiong, W.; Cheng, K.; Cui, T.; Godlewski, G.; Rice, K.; Xu, Y.; Zhang, L. Cannabinoid potentiation of glycine receptors contributes to cannabis-induced analgesia. Nat. Chem. Biol. 2011, 7, 296–303. [Google Scholar] [CrossRef] [Green Version]
  85. Anavi-Goffer, S.; Baillie, G.; Irving, A.J.; Gertsch, J.; Greig, I.R.; Pertwee, R.G.; Ross, R.A. Modulation of L-α-lysophosphatidylinositol/GPR55 mitogen-activated protein kinase (MAPK) signaling by cannabinoids. J. Biol. Chem. 2012, 287, 91–104. [Google Scholar] [CrossRef] [Green Version]
  86. Hurst, D.P.; Grossfield, A.; Lynch, D.L.; Feller, S.; Romo, T.D.; Gawrisch, K.; Pitman, M.C.; Reggio, P.H. A lipid pathway for ligand binding is necessary for a cannabinoid G protein-coupled receptor. J. Biol. Chem. 2010, 285, 17954–17964. [Google Scholar] [CrossRef] [Green Version]
  87. Picone, R.P.; Khanolkar, A.D.; Xu, W.; Ayotte, L.A.; Thakur, G.A.; Hurst, D.P.; Abood, M.E.; Reggio, P.H.; Fournier, D.J.; Makriyannis, A. (-)-7′-Isothiocyanato-11-hydroxy-1′,1′-dimethylheptylhexahydrocannabinol (AM841), a high-affinity electrophilic ligand, interacts covalently with a cysteine in helix six and activates the CB1 cannabinoid receptor. Mol. Pharmacol. 2005, 68, 1623–1635. [Google Scholar] [CrossRef] [Green Version]
  88. Pei, Y.; Mercier, R.W.; Anday, J.K.; Thakur, G.A.; Zvonok, A.M.; Hurst, D.; Reggio, P.H.; Janero, D.R.; Makriyannis, A. Ligand-binding architecture of human CB2 cannabinoid receptor: Evidence for receptor subtype-specific binding motif and modeling GPCR activation. Chem. Biol. 2008, 15, 1207–1219. [Google Scholar] [CrossRef] [Green Version]
  89. Palczewski, K.; Kumasaka, T.; Hori, T.; Behnke, C.A.; Motoshima, H.; Fox, B.A.; Le Trong, I.; Teller, D.C.; Okada, T.; Stenkamp, R.E.; et al. Crystal structure of rhodopsin: A G protein-coupled receptor. Science 2000, 289, 739–745. [Google Scholar] [CrossRef] [Green Version]
  90. Ring, A.M.; Manglik, A.; Kruse, A.C.; Enos, M.D.; Weis, W.I.; Garcia, K.C.; Kobilka, B.K. Adrenaline-activated structure of β2-adrenoceptor stabilized by an engineered nanobody. Nature 2013, 502, 575–579. [Google Scholar] [CrossRef] [Green Version]
  91. Shim, J.-Y.; Bertalovitz, A.C.; Kendall, D.A. Identification of essential cannabinoid-binding domains: Structural insights into early dynamic events in receptor activation. J. Biol. Chem. 2011, 286, 33422–33435. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  92. McAllister, S.D.; Tao, Q.; Barnett-Norris, J.; Buehner, K.; Hurst, D.P.; Guarnieri, F.; Reggio, P.H.; Nowell Harmon, K.W.; Cabral, G.A.; Abood, M.E. A critical role for a tyrosine residue in the cannabinoid receptors for ligand recognition. Biochem. Pharmacol. 2002, 63, 2121–2136. [Google Scholar] [CrossRef]
  93. Murphy, J.W.; Kendall, D.A. Integrity of extracellular loop 1 of the human cannabinoid receptor 1 is critical for high-affinity binding of the ligand CP 55,940 but not SR 141716A. Biochem. Pharmacol. 2003, 65, 1623–1631. [Google Scholar] [CrossRef]
  94. Munk, C.; Mutt, E.; Isberg, V.; Nikolajsen, L.F.; Bibbe, J.M.; Flock, T.; Hanson, M.A.; Stevens, R.C.; Deupi, X.; Gloriam, D.E. An online resource for GPCR structure determination and analysis. Nat. Methods 2019, 16, 151–162. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Structures of proposed GPR18 agonists (13) and antagonists (46).
Figure 1. Structures of proposed GPR18 agonists (13) and antagonists (46).
Biomolecules 10 00686 g001
Figure 2. Proposed binding mode of antagonist 4. (A) Docked pose of 4 in complex with the homology model of the human GPR18 shown with the residues forming the binding pocket. The receptor is displayed in cartoon representation, the amino acid residues (white) and compound 4 (orange) are shown as stick models. Oxygen atoms are colored in red, nitrogen atoms in blue, chlorine in green and sulfur atoms in yellow. (B) Schematic 2D representation of the binding pocket. Lipophilic amino acids are colored in yellow, hydrophilic ones in blue, aromatic ones in red, amino acid residues with mixed properties in green. (C) Schematic presentation of the homology model of GPR18 in complex with antagonist 4. The imidazothiazinone moiety is predicted to bind in the hydrophobic binding pocket consisting of residues of TM3 and ECL2. The 4-chlorophenyl moiety binds in the aromatic binding pocket consisting of residues of TM6 and TM7. Cys2516.54 in the aromatic binding pocket most likely interacts with hydrophobic substituents in position 4 of the phenoxy (4) or benzyloxy (5) moiety of the antagonists.
Figure 2. Proposed binding mode of antagonist 4. (A) Docked pose of 4 in complex with the homology model of the human GPR18 shown with the residues forming the binding pocket. The receptor is displayed in cartoon representation, the amino acid residues (white) and compound 4 (orange) are shown as stick models. Oxygen atoms are colored in red, nitrogen atoms in blue, chlorine in green and sulfur atoms in yellow. (B) Schematic 2D representation of the binding pocket. Lipophilic amino acids are colored in yellow, hydrophilic ones in blue, aromatic ones in red, amino acid residues with mixed properties in green. (C) Schematic presentation of the homology model of GPR18 in complex with antagonist 4. The imidazothiazinone moiety is predicted to bind in the hydrophobic binding pocket consisting of residues of TM3 and ECL2. The 4-chlorophenyl moiety binds in the aromatic binding pocket consisting of residues of TM6 and TM7. Cys2516.54 in the aromatic binding pocket most likely interacts with hydrophobic substituents in position 4 of the phenoxy (4) or benzyloxy (5) moiety of the antagonists.
Biomolecules 10 00686 g002
Figure 3. Proposed binding mode of antagonist 5. (A) Docked pose of 5 in complex with the homology model of GPR18 shown with the residues forming the binding pocket. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2.
Figure 3. Proposed binding mode of antagonist 5. (A) Docked pose of 5 in complex with the homology model of GPR18 shown with the residues forming the binding pocket. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2.
Biomolecules 10 00686 g003
Figure 4. (a,b) Root mean square deviation (RMSD) curves for the 200 ns MD simulation runs of the GPR18 complex with antagonist 4 (a) and antagonist 5 (b). (c,d) Root mean square fluctuation (RMSF) curves of the molecular dynamics (MD) simulation for complexes with antagonist 4 (c) and 5 (d). Curves of the complexes are colored in orange, and the curve of the apo form of the receptor in black.
Figure 4. (a,b) Root mean square deviation (RMSD) curves for the 200 ns MD simulation runs of the GPR18 complex with antagonist 4 (a) and antagonist 5 (b). (c,d) Root mean square fluctuation (RMSF) curves of the molecular dynamics (MD) simulation for complexes with antagonist 4 (c) and 5 (d). Curves of the complexes are colored in orange, and the curve of the apo form of the receptor in black.
Biomolecules 10 00686 g004
Figure 5. RMSD values for the transmembrane (TM) regions during the 200 ns MD simulation runs. Values were calculated based on the initial complex state after equilibration (0 ns).
Figure 5. RMSD values for the transmembrane (TM) regions during the 200 ns MD simulation runs. Values were calculated based on the initial complex state after equilibration (0 ns).
Biomolecules 10 00686 g005
Figure 6. Proposed binding mode of antagonist 6. (A) Docked pose of 6 in complex with the homology model of human GPR18 shown with the residues forming the binding pocket. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2. (C) Overlay of the proposed binding modes of GPR18 antagonists. Antagonist 4 is colored in orange, antagonist 5 in red, antagonist 6 in blue.
Figure 6. Proposed binding mode of antagonist 6. (A) Docked pose of 6 in complex with the homology model of human GPR18 shown with the residues forming the binding pocket. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2. (C) Overlay of the proposed binding modes of GPR18 antagonists. Antagonist 4 is colored in orange, antagonist 5 in red, antagonist 6 in blue.
Biomolecules 10 00686 g006
Figure 7. Schematic representation of the structure–activity relationships (SARs) of GPR18 antagonists. The different heterocycles shown contain a 4-chlorophenoxy group, while the compounds with varying aryl substituents and linker lengths contain the imidazothiazinone heterocycle. Compounds were categorized into three groups: highest potency (IC50 < 1 µM), moderate potency (1 µM < IC50 < 10 µM) and low potency (IC50 > 10 µM) based on their antagonistic activity.
Figure 7. Schematic representation of the structure–activity relationships (SARs) of GPR18 antagonists. The different heterocycles shown contain a 4-chlorophenoxy group, while the compounds with varying aryl substituents and linker lengths contain the imidazothiazinone heterocycle. Compounds were categorized into three groups: highest potency (IC50 < 1 µM), moderate potency (1 µM < IC50 < 10 µM) and low potency (IC50 > 10 µM) based on their antagonistic activity.
Biomolecules 10 00686 g007
Figure 8. Proposed binding mode of ∆9-tetrahydrocannabinol (THC) in the homology model of human GPR18. (A) The receptor is displayed in cartoon representation, the amino acid residues (white) and THC (1, green) are shown as stick models. Oxygen atoms are colored in red, nitrogen atoms in blue, sulfur atoms in yellow. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2.
Figure 8. Proposed binding mode of ∆9-tetrahydrocannabinol (THC) in the homology model of human GPR18. (A) The receptor is displayed in cartoon representation, the amino acid residues (white) and THC (1, green) are shown as stick models. Oxygen atoms are colored in red, nitrogen atoms in blue, sulfur atoms in yellow. (B) Schematic 2D representation of the binding pocket. For color code, see Figure 2.
Biomolecules 10 00686 g008
Figure 9. (A) Comparison of the proposed binding mode of THC (green) and antagonist 4 (orange) at human GPR18. (B) Comparison of the proposed binding modes of THC and antagonist 5 (red).
Figure 9. (A) Comparison of the proposed binding mode of THC (green) and antagonist 4 (orange) at human GPR18. (B) Comparison of the proposed binding modes of THC and antagonist 5 (red).
Biomolecules 10 00686 g009

Share and Cite

MDPI and ACS Style

Neumann, A.; Engel, V.; Mahardhika, A.B.; Schoeder, C.T.; Namasivayam, V.; Kieć-Kononowicz, K.; Müller, C.E. Computational Investigations on the Binding Mode of Ligands for the Cannabinoid-Activated G Protein-Coupled Receptor GPR18. Biomolecules 2020, 10, 686. https://0-doi-org.brum.beds.ac.uk/10.3390/biom10050686

AMA Style

Neumann A, Engel V, Mahardhika AB, Schoeder CT, Namasivayam V, Kieć-Kononowicz K, Müller CE. Computational Investigations on the Binding Mode of Ligands for the Cannabinoid-Activated G Protein-Coupled Receptor GPR18. Biomolecules. 2020; 10(5):686. https://0-doi-org.brum.beds.ac.uk/10.3390/biom10050686

Chicago/Turabian Style

Neumann, Alexander, Viktor Engel, Andhika B. Mahardhika, Clara T. Schoeder, Vigneshwaran Namasivayam, Katarzyna Kieć-Kononowicz, and Christa E. Müller. 2020. "Computational Investigations on the Binding Mode of Ligands for the Cannabinoid-Activated G Protein-Coupled Receptor GPR18" Biomolecules 10, no. 5: 686. https://0-doi-org.brum.beds.ac.uk/10.3390/biom10050686

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop