A simple open source bioinformatic methodology for initial exploration of GPCR ligands’ agonistic/antagonistic properties

Abstract Drug development is an arduous procedure, necessitating testing the interaction of a large number of potential candidates with potential interacting (macro)molecules. Therefore, any method which could provide an initial screening of potential candidate drugs might be of interest for the acceleration of the procedure, by highlighting interesting compounds, prior to in vitro and in vivo validation. In this line, we present a method which may identify potential hits, with agonistic and/or antagonistic properties on GPCR receptors, integrating the knowledge on signaling events triggered by receptor activation (GPCRs binding to Gα,β,γ proteins, and activating Gα, exchanging GDP for GTP, leading to a decreased affinity of the Gα for the GPCR). We show that, by integrating GPCR‐ligand and Gα‐GDP or ‐GTP binding in docking simulation, which correctly predicts crystallographic data, we can discriminate agonists, partial agonists, and antagonists, through a linear function, based on the ΔG (Gibbs‐free energy) of liganded‐GPCR/Gα‐GDP. We built our model using two Gαs (β2‐adrenergic and prostaglandin‐D2), four Gαi (μ‐opioid, dopamine‐D3, adenosine‐A1, rhodopsin), and one Gαo (serotonin) receptors and validated it with a series of ligands on a recently deorphanized Gαi receptor (OXER1). This approach could be a valuable tool for initial in silico validation and design of GPRC‐interacting ligands.

testing of an increasing number of chemical libraries for positive hits and the subsequent biological validation of promising candidates. Therefore, any progress leading to an initial discrimination of novel potential bioactive compounds could lead to a more accurate identification of possible novel drug candidates.
Of the 59 drugs approved in 2018, 23 target membrane components/receptors, whereas 11 target G-protein-coupled receptors (GPCR). 1 Membrane receptors act as conveyors of extracellular signals into the cell. Membrane receptors can be distinguished as onepass single-chain proteins, acting as mono-or oligomers and multiple passes proteins. Among the latter, the seven transmembrane helix (7TM) GPCR family contains ~800 members overall (of which ~400 are olfactory receptors). GPCRs are involved in different signal transduction pathways, triggered by a plurality of extracellular signals (including photons, light-sensitive compounds, photons, odorants, pheromones, hormones, neurotransmitters, and a number of ligands, varying in size from small molecules to peptides to large (glycol) proteins). GPCR-initiated signal transduction results in many physiological processes, interfering with the (patho)physiology of many systems, such as the endocrine (including the reproductive), neurological or cardiovascular systems. Such a wide impact, makes GPCR a preferential drug target candidate group. [2][3][4] Indeed, GPCR-interacting drugs account for ~34% of the global market share. 3 However, only a small fraction of GPCRs (206 entries according to https://gpcrdb.org/ struc ture/stati stics) have been crystallized to date, making difficult the prediction of novel pharmacological substances.

Molecular docking plays a major role in identifying molecules
that might fulfill the requirements of drug development. The applied methodologies simulate the interaction of ligands (small molecules or peptides) with corresponding receptors, in monomeric or oligomeric states. The derived solutions are represented as scoring function (usually reported as the difference in Gibbs-free energy for molecular association, denoted as ΔG, in kcal/mol, relying on the enthalpy, the entropy and the temperature of the complex), which allows the evaluation of the ligand interaction with the receptor. 5 In recent years, an increased number of commercial and open source software has been released (see, 6 for a recent discussion of available resources, and, 7 for open source solutions). Furthermore, the existence and release of open libraries with chemical structures also accelerated the implementation of this process. 8 In the field of GPCR pharmacology, the integration of GPCRligand interactions (see, 9 for a recent example analysis) resulted in a high success rate of GPCR-targeting ligands, translated in successful drug design and achieving 78% success rate in Phase I, 39% in Phase II and 29% in Phase III clinical trials 2 (see also, 10 for a successful recent paradigm). However, in spite of the identification of GPCR downstream signaling events (G α,β,γ complex, or arrestin signaling) triggered by receptor activation (see, [11][12][13] for recent reviews), no attempts have been made to integrate such a knowledge into the search of novel pharmacophores or drugs. Here, we have developed and present a pipeline for GPCR-ligand interactions and candidate identification, based on free online resources and programs, that also integrates the subsequent steps of molecular docking to GDP-and/or GTP-linked G α protein binding. We show that it can correctly predict small ligand putative agonistic or antagonistic nature, presenting a valuable tool that could significantly accelerate the search of novel molecules in GPCR pharmacology.

| In silico methods
Our approach consists of three sequential phases: (A) Ligand and receptor preparation, (B) ligand-receptor docking and (C) G α -protein interaction ( Figure 1 and Supporting Information 2, which provides an illustrated User's Manual).

| Ligand and Receptor Preparation
• For receptor preparation: the sequence of human receptors, in fasta format, were retrieved from the NCBI protein database (https://www.ncbi.nlm.nih.gov/prote in/) and introduced to the Swiss Model Biospace (http://swiss model.expasy.org/inter active). 14 • If a crystalized receptor file (bound or not with an agonist or an antagonist) was available, the system returned the code (and corresponding pdb 3D coordinates file of the receptor or its complex with an agonist or an antagonist, from data stored in the Protein Data Bank (https://www.rcsb.org/). 15 Whenever the structure contained a ligand, the receptor structure was manually extracted (using a text editor) from the returned pdb file. This did not interfere with the subsequent (flexible) binding, as a full backbone and side chain flexible binding was performed (see below). All protein pdb codes, used in this study, are presented in Table S1.
• If a crystalized structure was not available, the 3D structure of the receptor was simulated by molecular modeling calculations.
In this case, the fasta receptor file was introduced in the Swiss Model Biospace (http://swiss model.expasy.org/inter active), 14,16 which returned a series of files/crystalized solutions with a variable coverage homology based on the sequence identity to the file in question. We have retained solutions with a coverage homology ≥ 50%-70%

| G-protein interactions
A subsequent step of GPCR-ligand activation is the binding with G-proteins, specific sites. More specifically, G α is bound to intracellular activated receptor loops 2 and 3, and G β,γ is bound to intracellular chain 8, 25,26 initiating specific intracellular signaling events. In this work, we have examined the interaction of GPCRs with G α proteins. It is to note that, after binding of G α protein (bound to GDP and denoted here as G α -GDP), the nucleotide is exchanged, after receptor activation, to GTP (G α -GTP), and the G α -GTP protein is liberated (due to a decrease of affinity for the GPCR) and subsequently triggers specific signaling pathways.
Here, we simulated the interaction of known and novel GPCRs with G α proteins.
At a first step, we have retrieved the sequences of G α -proteins, in fasta format, from the NCBI protein database, and after a SwissDock generation of a 3D structure, with known crystal templates, were introduced to the GalaxyWEB to generate and refine the structures, as discussed above for the receptor files. This step was necessary, as the reported crystal structures of the different Hex returns, through a graphical user interface, a set of > 100 solutions, with the corresponding ΔG values. We have manually inspected and retained only solutions (usually scored first) in which G α molecules bind to GPCRs intracellular loops 2 and 3.

| Validation of the obtained solutions
The obtained GPCR-G α models, obtained from the above procedure, were compared with the reported structures of liganded receptor-G α proteins. We have retrieved data for the liganded β-adrenergic receptor (PDB code 3SN6), 28 the μ-opioid receptor (PDB code 6DDE), 29 the rhodopsin receptor (PDB code 6CMO), 30 the serotonin receptor (PDB code 6G79), 31 the adenosine A1 receptor (PDB code 6D9H), 32 co-crystalized with corresponding Gα proteins. Data were inspected in UCSF Chimera program, by superpositioning of the two structures and the corresponding total and local RMSD value (in Å) were retrieved, with Needleman-Wunsch alignment 33 and with the use of BLOSUM-62 matrix. 33

| In vitro validation assay
As our goal was to use the proposed algorithm as a prediction tool for the agonistic or antagonistic character of novel ligands on specific GPCRs, we have further validated our in silico results, by exploring the interaction of a series of pregnenolone analogs 34 and polyphenol molecules 35 , as agonists or antagonists of the novel deorphanized GPCR OXER1. 36,37 OXER1 is an oxo-eicosanoid receptor, on which 5-oxo-ETE is reported to be the physiological agonist and which can also bind other oxo-eicosanoids, products of arachidonic acid cellular transformation. Recently, we have identified this receptor as a membrane androgen binding site, 38 with testosterone acting antagonistically on cAMP production and kinases signaling. It is to note that OXER1 binds to a G αi protein and decreases intracellular cAMP production, whereas testosterone, in equimolar concentration, reverts this inhibition by ~50%. 38 Therefore, in order to validate our in silico data, we have assayed cAMP production in DU145 human prostate cancer cells, bearing OXER1, according to our previous report. Results were expressed as % reversion of the 5-oxo-ETE effect, in the presence of forskolin. 38

| Statistical Analysis
Discriminant analysis was performed with the SPSS V21 program (IBM, SPSS Statistics), whereas group comparisons were made with the GraphPad Prism V6.0.5 (GraphPad Software Inc). A statistical threshold of P < .05 was retained for significance. Results (as changes of the Gibbs-free energy changes, ΔG, in kcal/mol) are shown in Tables 1-3, column GPCR-Ligand.

Interaction of the ligand-receptor complex with G α -proteins
A subsequent step following ligand-GPCR interaction is the binding of the liganded receptor to the heteroprotein complexG α,β,γ , 11-13 triggering specific signaling events. G α proteins interact with intracellular loops 2 and 3, whereas G β,γ with intracellular loop 3 and the intracel-  TA B L E 2 Fully flexible ligand binding results on dopamine D 3 (pdb 3PBL), μ-opioid (pdb 5C1M), Adenosine A 1 (pdb 6D9H), and Rhodopsin (pdb 6CMO) receptors, together with the liganded GPCR (GPCR(L)) binding to G αi in its GDP and GTP-bound forms. All data are reported as differences in the Gibbs-free energy (ΔG), expressed in kcal/mol. The effect column presents the reported action of the compound (bibliography) and the predicted effect by the proposed model (Model the G α proteins (Tables 1-3), the interaction of G α proteins with the receptor was significantly decreased. In this case, the obtained values do not differ among agonists, partial agonists, and antagonists ( Figure 2).

Verification of the GPCR-ligand binding and G α interaction
In In the latter, a very good match was found at the interacting part of the G α protein, whereas the observed differences in G α proteins might be attributed to the recalculation of these proteins' 3D structure, due to major disruptions of the protein structures in the G α crystal.
Finally, we have tried to simulate G α -GDP binding to an unrelated multipass membrane protein (aquaporin monomer, extracted

TA B L E 3
Fully flexible ligand binding results on the serotonin receptor (pdb 6G79), with the liganded GPCR (GPCR(L)) binding to G αo in its GDP and GTP-bound forms. All data are reported as differences in the Gibbs-free energy (ΔG), expressed in kcal/mol. The effect column presents the reported action of the compound (bibliography) and the predicted effect by the proposed model (Model). See text for details from pdb 6KXW). All three G α -GDP complexes did not interact with this protein monomer, corroborating about the specificity of the simulated interaction.

EFFECT Bibliography/Model
In view of the above, we have concluded that our approach may indeed correctly simulate the interaction of known GPCRs with their corresponding G α -proteins.

G α binding to the liganded receptor can discriminate between GPCR agonists and antagonists
After determining the different ΔGs for ligand-GPCR binding and for liganded GPCR-liganded G α interaction (presented in   (presented in Tables 1-3), upon agonist, partial agonist (PA) and antagonist binding. In (A), the negative ΔG GPCR(L)-G α GDP value is shown, whereas in (B) the corresponding negative ΔG GPCR(L)-G α GTP value is depicted. Post hoc group comparisons were made after ANOVA, with the Turkey's multiple comparison test, in GraphPad Prism V6 protein used is denoted in parentheses in the first column of Table 4).

Applying our cut-off values −666 and −887 kcal/mol for GDP-bound
Gα protein, we show that we have correctly identified the three agonistic drugs.

| Detection of agonists and antagonists for a novel receptor (OXER1)
A valid prediction method should provide useful hints about the agonistic or antagonistic properties of both crystalized or not GPCRs.
Hence, we used our approach to predict the agonistic-antagonistic properties of a number of substances, of very different molecular structure (lipids, steroids, polyphenols), on the oxo-eicosanoid receptor OXER1.
OXER1 was deorphanized in 2002-3 and was found to be the endogenous receptor for the arachidonic acid metabolic product 5-oxo-ETE, produced through the action of 5-lipoxygenase (5-LOX) and peroxidase. 36,37 However, recently, we have reported that OXER1, coupled to G αi protein, also mediates membrane-initiated androgen actions (see, 38 (Table 5). We have also calculated the affinity of a series of derivatives of arachidonic acid biotransformation, which have been previously reported to act as partial OXER1 agonists (see https://genec ards.weizm ann.ac.il/v3/cgi-bin/cardd isp. pl?gene=OXER1 and references therein). Obtained ΔG GPCR(L)-G α GDP values are intermediate between 5-oxo-ETE and testosterone, verifying their partial agonistic nature (Table 5).
In addition to the above compounds, we have tested a series of pregnenolone analogs, with reported antiproliferative activity in different cancer cell lines. 34 As shown in Table 5 In order to verify our prediction, we have experimentally tested whether these compounds can antagonize 5-oxo-ETE action on cAMP production, like testosterone 38 OXER1-G αi interaction results in an inhibition of cAMP. [36][37][38] This is experimentally tested by stimulating cAMP production in cells by forskolin and detecting the cAMP inhibition after incubation of cells with the corresponding ligands. We have previously shown that testosterone incubation of prostate cancer cells reverts the 5-oxo-ETE-induced inhibition, in a dose-dependent manner. 38 Here, we have applied the same protocol using pregnenolone analogs and polyphenols, after forskolin stimulation of DU145 human prostate cancer cells and application of 5-oxo-ETE. Table 5 presents the normalized cAMP inhibition (5-oxo-ETE inhibition = 100%). Testosterone reverts this inhibition by 51%, at a concentration 1 μM, as reported previously. 38 Of the tested compounds, all reverted 5-oxo-ETE cAMP inhibition by 48%-67%, at the same 1 μM concentration, classifying them as antagonists of OXER1, with the notable exception of B5 procyanidin, which reverted 5-oxo-ETE cAMP inhibition by only 25%, classifying it as a partial agonist, as also suggested by the in silico binding data.

| D ISCUSS I ON
Drug development is a laborious procedure, necessitating the testing the interaction of a large number of potential candidates, with potential target (macro) molecules. Therefore, any method which TA B L E 4 Simulation data of four novel compounds approved by the FDA in 2018 (1). For each compound its affinity (Galaxy Docking) and its interaction with the corresponding G α protein are shown. The interaction of each drug with the GPCR analyzed here and with its cognate receptor (for which an FDA approval was provided) is shown. X denotes nonassociation could provide an initial screening of chemicals as positive hits, might be of interest for the selection of interesting compounds, which could decrease the time-frame in drug discovery, prior to in vitro and in vivo validation. Here, we report a method (see Figure 1 for a schematic representation) which may be used for the initial, in silico screening of potentially active compounds, taking into account the binding of the ligand on the corresponding receptor and its subsequent simulated affinity for G α -GDP. We report that the latter may correctly discriminate ~90% of substances between agonists, partial agonists and antagonists.
GPCRs-related drugs account for 34% of all drug targets. [2][3]43 In addition, several pharmacological substances, designed to interact with a single target, were found to mediate effects via several The novelty of our approach relies on exploiting, in addition to ligand-GPCR fully flexible docking, an initial step of the subsequent signaling event, their interaction with G α -proteins, 25 to provide a quick initial estimate of ligand agonistic or antagonistic properties. In our analysis, agonistic ligands induce a significantly higher affinity for the liganded receptor G α -GDP interaction.
This affinity decreases substantially when the same G-protein is bound to GTP, expressing the biologically relevant dissociation of the GTP-bound G-protein from the receptor and the initiation of intracellular signaling events. 25 Our approach is based on bibliographic data from known ligand interactors of crystalized or TA B L E 5 Fully flexible ligand binding results on the OXER1 receptor, together with the liganded GPCR (GPCR(L)) binding to G αi in its GDP-and GTP-bound forms. All data are reported as differences in the Gibbs-free energy (ΔG), expressed in kcal/mol. Data from 5-HETE, 12-HpETE, 15-HpETE, 12-HETE, and 15-HETE were from previous studies, and extracted from the Gene Cards web site, whereas data for all other compounds were experimentally verified, through an inhibition of 5-oxo-ETE effect on cAMP production. Here, the maximum inhibition of forskolin stimulated inhibition of cAMP production by 1 μM 5-oxo-ETE (the natural ligand of OXER1 receptor) was set as 100% inhibition, and data obtained by all other compounds were compared to this maximum value at a similar 1 μM concentration, added simultaneously with 5-oxo-ETE. Please refer to the Material and Methods section, to Figure 3E and text of reference (38), and to Figure S3 for further details. The effect column presents the reported action of the compound (bibliography), the experimental validation (experimental), and the predicted effect by the proposed model (Model Abbreviations: Ago, Agonist; Antago, Antagonist; GPCR(L), Ligand-bound GPCR; NA, Non-available; PA, Partial Agonist. a Mean ± SE, n = 3.

| CON CLUS ION
Our data clearly show that, by integrating sequential steps of receptor downstream signaling in ligand-GPCR simulations, as expressed by GDP-G α binding, we can correctly predict the nature (agonist, antagonist, partial agonist) of a given small molecule. This approach, combined to properly implemented and successfully validated QSAR methods, 45 may represent a useful addition to current research processes for the initial prediction and design of novel GPRC-interacting molecules. It might be of interest to explore further whether similar initial estimates might be also applied on other, non-GPCR, receptors, which could provide a generalization of our approach.

ACK N OWLED G M ENTS
This work was partially supported by Greece and the European Union

D I SCLOS U R E S
None declared.

AUTH O R S' CO NTR I B UTI O N S
MK and EC conceived and designed the study, and wrote the paper.
AP, CP, KK, and PM performed the analyses and experiments. AP, TC, GN, and PAT participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.

E TH I C A L S TATEM ENT
This article does not contain any studies involving animals or human participants performed by any of the authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request. Some data may not be made available because of privacy or ethical restrictions.