- Research article
- Open Access
Discriminating agonist and antagonist ligands of the nuclear receptors using 3D-pharmacophores
Journal of Cheminformaticsvolume 8, Article number: 43 (2016)
Nuclear receptors (NRs) constitute an important class of therapeutic targets. We evaluated the performance of 3D structure-based and ligand-based pharmacophore models in predicting the pharmacological profile of NRs ligands using the NRLiSt BDB database. We could generate selective pharmacophores for agonist and antagonist ligands and we found that the best performances were obtained by combining the structure-based and the ligand-based approaches. The combination of pharmacophores that were generated allowed to cover most of the chemical space of the NRLiSt BDB datasets. By screening the whole NRLiSt BDB on our 3D pharmacophores, we demonstrated their selectivity towards their dedicated NRs ligands. The 3D pharmacophores herein presented can thus be used as a predictor of the pharmacological activity of NRs ligands.
Nuclear receptors (NRs) are involved in a wide range of physiological key functions. They are potential targets for numerous diseases and constitute an important class of therapeutic targets [1, 2]. NRs are transcription factors naturally switched on and off by small-molecule hormones, and artificially by synthetic ligands. Taking advantage of the biological potency of the NRs, a large amount of compounds has been proposed to modulate their activity and some of them are still marketed [3, 4]. The NRs ligands can be classified according to their pharmacological profiles, the two main classes being agonist and antagonist ligands. These two classes of compounds act through the binding to a NR and the activation (agonist ligands) or the inhibition (antagonist ligands) of its activity. The drug discovery process is thus not limited to the search of the best ligand of a given target, but consists in the search of a ligand with a pharmacological profile that his compatible to the required activity. In this context, the ability to predict the agonist or antagonist behaviour of a NR ligand is of major importance. In recent years, virtual screening methods have proven their ability to predict the activity of small compounds [5–7] and can be used to predict the pharmacological profile of NRs ligands. Numerous ligand-based (LB) and structure-based (SB) virtual screening studies dedicated to NRs were conducted but only few focused on the agonism/antagonism issues [8–17]. Despite these several prediction attempts and the elucidation of the molecular bases of agonism and antagonism [18–21], discriminating agonist from antagonist ligands based on their sole structure remains a challenge. In this study, we describe a 3D pharmacophore modeling study performed on 27 NRs, with the aim to provide separate and selective agonist and antagonist pharmacophores for each NR. To our knowledge, this is the first large-scale study conducted to predict the agonist and antagonist behaviour of NRs ligands using a 3D pharmacophore modeling method. 3D pharmacophores are nowadays widely used as filters in virtual screening protocols and several studies successfully identified new NRs ligands using pharmacophore models [22–30]. Pharmacophore models display two main advantages: reduced computational times associated to the simplified pharmacophoric representations and a large diversity of potential hits with scaffolds and functional groups distinct to the original ligands [31, 32]. To design our study, we used the 27 NRLiSt BDB datasets . For each dataset, we created both SB and LB 3D pharmacophores and compared the ability of these two approaches to generate agonist selective pharmacophores and antagonist selective pharmacophores covering the whole NRLiSt BDB ligands chemical space. We also studied the performance obtained using a combination of SB and LB pharmacophores and analyzed the composition and the selectivity against all NRs datasets of these combinations. In the present study, we describe our attempt to develop selective pharmacophores for agonist or antagonist ligands that could be used to predict the pharmacological activity of NRs.
Nuclear receptors ligands and structures benchmarking DataBase (NRLiSt BDB)
The NRLiSt BDB  is a freely available benchmarking database for both SB and LB methods evaluation and dedicated to the NRs. The NRLiSt BDB presents separated agonist and antagonist datasets for the 27 targets (out of the 48 known NRs) for which more than one agonist ligand, one antagonist ligand, and at least one experimental structure was available. All of the ligands found to be agonist or antagonist in the scientific literature are provided in two separated datasets and all of the available human holo PDB structures (except for RXR_gamma, for which only one apo structure was available). A total of 7853 actives, 458,981 decoys, and 339 structures are divided into 54 datasets. The NRLiSt BDB was downloaded from the Web site http://nrlist.drugdesign.fr.
3D pharmacophores were generated using the software LigandScout  (version 4.0) in SB and LB approaches.
3D SB pharmacophores were automatically generated using the PDB structures included in the NRLiSt BDB. This approach is only possible with holo structures, thus no RXR gamma 3D SB pharmacophore could be computed. In this approach, the LigandScout algorithm tags the key features of the ligand that are interacting with the residues of the receptor: aromatic ring, hydrophobic area, hydrogen bond donor or acceptor, negative or positive ionisable atom and metal binding location. To complete the pharmacophore, an ensemble of exclusion volume spheres is generated to represent the shape of the active site.
All ligands of each dataset were clustered with LigandScout using default settings except for the cluster distance that was adjusted for each NR to obtain balanced clusters. For each cluster, a 3D LB pharmacophore was generated using the “merged feature pharmacophore approach” with the number of omitted features for a given merged pharmacophore set to 4 and optional partially matching features with a threshold set to 10 %. In this approach, all the features observed in each ligand of the training datasets are identified, scored and removed according to the threshold number of omitted features. We chose to enable the creation of exclusion volume spheres around the alignment of ligands. In some cases, we added manually exclusion volume spheres to remove decoys compounds since inactive compounds can map all the pharmacophore require features, their inactivity being explained by steric clashes with the binding site [35, 36]. For each pharmacophore, the ligands of the cluster used to generate the pharmacophore constituted the training set and the test set was formed by all agonist ligands and all antagonist ligands of the corresponding NR. During the pharmacophore generation, the ligands of the training set were automatically aligned with the LigandScout pharmacophore-based alignment algorithm .
Model optimization protocol
The generated 3D pharmacophores were used to screen the NRs datasets. All of the ligands provided in SMILES format in the NRLiSt BDB were converted in .ldb format using the idbgen tool provided with LigandScout with the omega-fast option. Two databases were used for each screening, a screening database of active compounds and a screening database of decoys. Agonist ligands were used as decoys for antagonist pharmacophores and reciprocally antagonist ligands were used as decoys for agonist pharmacophores. We developed an original model optimization protocol for this study (Fig. 1), to sequentially refine the pharmacophore models according to several literature recommendations [32, 38, 39]. For each pharmacophore, a first screening was made with LigandScout default settings and particularly the Max. number of omitted features set to 0. If the hits retrieved with this first screening contained both agonist and antagonist ligands, the pharmacophore was not validated and was not retained. If only agonist or antagonist ligands were retrieved in this first screening, the pharmacophore was validated and a second screening was performed with this pharmacophore, but with the Max. number of omitted features parameter set to 1. This second screening was carried out to identify possible non-essential pharmacophore features, i.e. features that can be disabled to obtain less stringent pharmacophores able to retrieve more active ligands (agonist ligands when using agonist datasets or antagonist ligands when using antagonists datasets), but no decoys (antagonists when using agonist datasets and agonists when using antagonist datasets). When a non-essential pharmacophore feature was identified, a third screening was performed with the non-essential feature marked as disabled and the Max. number of omitted features parameter set to 0. If the hits retrieved with this third screening were both agonist and antagonist ligands, this second pharmacophore was not validated and another round of identification of non-essential features was performed. If only active ligands were retrieved, the pharmacophore was validated and other non-essential features were studied. This protocol was applied to each pharmacophore until 3 pharmacophore features were retained or until no non-essential feature could be identified.
Combination of SB pharmacophores, combination of LB pharmacophores and combination of SBLB pharmacophores
Using the SB approach, for each NR, all the selective pharmacophores generated, i.e. all pharmacophores that retrieved only agonist ligands or antagonist ligands, were gathered into two groups: “SB agonist selective pharmacophores” and “SB antagonist selective pharmacophores”; redundant pharmacophores were removed. Similarly, all the selective pharmacophores obtained with the LB approach for each NR were gathered into two groups: “LB agonist selective pharmacophores” and “LB antagonist selective pharmacophores”; redundant pharmacophores were removed. Finally, the SB and LB selective pharmacophores previously generated were gathered in two pharmacophore ensembles: “SBLB agonist selective pharmacophores” and “SBLB antagonist selective pharmacophores”; redundant pharmacophores were removed.
Redundant pharmacophores are pharmacophores that could be removed without decreasing the recall of the set of combined pharmacophores i.e. pharmacophores that only retrieved ligands that were also retrieved with other pharmacophores of the set. To remove these redundant pharmacophores, all generated pharmacophores were ranked according to the number of hits they retrieved. Then, each pharmacophore was removed sequentially, starting from the pharmacophore associated with the smallest number of hits. For each removal, the impact on the recall was evaluated. If the recall was not affected, the pharmacophore was dismissed and in the opposite, if the recall decreased, the pharmacophore was conserved.
All the graphs were produced with the statistical and graphical tool R (http://www.r-project.org/). The ggplot2 package was used to produce the barplot of Figs. 2, 3, 4, 5 and 6. The corrplot and RColorBrewer packages were used to produce the graph of pharmacophores selectivity using the recall (R) value (Fig. 9). For each dataset, the recall (R), the specificity (Sp) and the Matthew’s correlation coefficient (MCC) were computed as follows:
with TP the number of true positives (number of active compounds of the dataset retrieved as screening hits), FN the number of false negatives (number of active compounds of the dataset not retrieved as screening hits), TN the number of true negatives (number of inactive compounds of the dataset not retrieved as screening hits), FP (number of inactive compounds of the dataset retrieved as screening hits). As we chose to generate only selective agonist or antagonist pharmacophores, the number of FP was always equal to 0 and thus the Sp value was always equal to 1. Similarly, in the SB approach, when the number of TP was equal to 0, it was not possible to compute the MCC value (because the denominator value is equal to 0), and the MCC value was qualified as not determined (ND).
Structure-based pharmacophore modeling
338 3D SB pharmacophores were generated from the 339 PDB structures included in the NRLiSt BDB. The protein structures included in the NRLiSt BDB are classified according to the pharmacological profile of the ligand bound in the active site: 266 agonist-bound structures, 17 antagonist-bound structures, 55 other-bound structures (partial agonists, modulators, inverse agonists etc.). Since only 1 apo structure was available for RXR_gamma, no 3D SB pharmacophore could be generated and this NR was excluded for this part of the study. For respectively 25 and 10 out the 26 remaining NRs, at least one agonist-bound or one antagonist-bound structure was available. Using the screening protocol described in the “Methods” section, we succeeded in generated at least one pharmacophore that was selective for agonist ligands for 25 NRs out of the 26 used, and at least one pharmacophore selective for antagonist ligands for 9 NRs out of 26. As presented in the “Methods” section, all these pharmacophores were gathered into two groups: “SB agonist selective pharmacophores” and “SB antagonist selective pharmacophores”, and redundant pharmacophores were removed. The average recall for the “SB agonist selective pharmacophores” was of 55 %, ranging from 0 % for SF1 to 98 % for PPAR_gamma whereas the average recall for “SB antagonist selective pharmacophores” was of 8 %, ranging from 0 % for AR, CAR, ERR_alpha, FXR_alpha, LXR_alpha, LXR_beta PPAR_alpha, PPAR_beta, PPAR_gamma, PXR, RAR_beta, RAR_gamma, RXR_beta, SF1, TR_alpha, TR_beta and VDR to 61 % for RAR_alpha (Fig. 2; Table 1). The average MCC value for the “SB agonist selective pharmacophores” was of 0.484, ranging from 0.088 for CAR (the SF1 MCC value was ND) to 0.881 for RXR_alpha whereas the average MCC value for “SB antagonist selective pharmacophores” was of 0.326, ranging from 0.097 for RXR_alpha (the MCC value was ND for AR, CAR, ERR_alpha, FXR_alpha, LXR_alpha, LXR_beta PPAR_alpha, PPAR_beta, PPAR_gamma, PXR, RAR_beta, RAR_gamma, RXR_beta, SF1, TR_alpha, TR_beta and VDR) to 0.712 for RAR_alpha (Table 1).
Ligand-based pharmacophore modeling
To perform the LB pharmacophore modeling approach, the ligands of each NRLiSt BDB dataset were clustered using the Pharmacophore RDF-Code similarity. The cluster distance was set to 0.4 for the majority of the datasets but was lowered to 0.3 for 15 datasets (AR_agonist, ERR_alpha_agonist, GR_agonist, LXR_alpha_agonist, LXR_beta_agonist, PR_agonist, RXR_alpha_agonist, RXR_beta_agonist, RXR_beta_antagonist, RXR_gamma_agonist, TR_alpha_agonist, TR_alpha_antagonist, TR_beta_agonist, TR_beta_antagonist, VDR_antagonist) and to 0.2 for 1 dataset (RXR_alpha_antago). From 1 cluster (for the ERR_alpha_agonist, ROR_gamma_antagonist, and RXR_gamma_antagonist datasets) to 65 clusters (for the ER_alpha_agonist dataset) were generated, with an average of 18 clusters per dataset and a mean value of 7.8 ligands per cluster.
3D ligand-based pharmacophores
Using the screening protocol described in the “Methods” section, we succeeded in generated pharmacophores that were selective for agonist ligands and pharmacophores selective for antagonist ligands for each of the 27 NRs of the NRLiSt BDB. All these pharmacophores were gathered into two groups according to their selectivity for agonist or antagonist ligands, “LB agonist selective pharmacophores” and “LB antagonist selective pharmacophores”. Redundant pharmacophores were eliminated. The “LB agonist selective pharmacophores” were associated with an average recall of 97 % and a mean value of 0.918. The lower recall and MCC value were respectively of 88 % for TR_alpha and 0.253 for PPAR_gamma; the higher recall and MCC values respectively reached 100 % and 1 for CAR, ERR_alpha, PPAR_beta, RAR_beta, ROR_alpha, ROR_gamma, RXR_gamma and SF1. The “LB antagonist selective pharmacophores” presented an average recall of 99 % and a mean MCC value of 0.99, and the individual recall and MCC values were equal to 100 % and 1 for all antagonist datasets but 5 (ER_alpha, ER_beta, GR, PR, RXR_alpha) (Fig. 3; Table 1).
Combination of structure-based and ligand-based pharmacophores
3D SBLB pharmacophores performance
The “SB agonist selective pharmacophores” and “LB agonist selective pharmacophores” on the one hand and the “SB antagonist selective pharmacophores” and “LB antagonist selective pharmacophores” on the other hand were respectively concatenated into two groups: “SBLB agonist selective pharmacophores” and “SBLB antagonist selective pharmacophores”. Using these combinations of SB and LB pharmacophores, average recalls of 99.7 and 99.9 % and mean MCC values of 0.993 and 0.999 were obtained for agonist and antagonist datasets respectively. The “SBLB agonist selective pharmacophores” were able to retrieve all agonist ligands and no antagonist ligands (i.e. recall of 100 % and MCC values of 1) for all NRs but 5 (LXR_alpha, LXR_beta, PR, RAR_alpha and RAR_gamma). Similarly, the “SBLB antagonist selective pharmacophores” were able to retrieve all antagonist ligands and no agonist ligands (i.e. recall of 100 % and MCC values of 1) for all NRs but 3 (ER_alpha, GR, PR) (Fig. 4; Table 1).
The “SBLB agonist selective pharmacophores” group contained 413 pharmacophores (from 1 pharmacophore for ROR_gamma to 52 pharmacophores for PR) whereas the “SBLB antagonist selective pharmacophores” group contained 305 pharmacophores (from 1 pharmacophore for CAR, ERR_alpha, ROR_gamma and RXR_gamma to 64 pharmacophores for PR) (Fig. 5). The number of pharmacophores that were necessary to cover a given dataset is significantly correlated with the number of ligands in the dataset (Kendall’s tau coefficient, p value = 9.55e−15, Additional file 1: Figure S1). These pharmacophores were composed of 3–16 features, with a median value of 5 features per pharmacophore (Fig. 6; Additional file 1: Figure S2A−N, Additional file 1: Tables S1−S54). Pharmacophore features were mainly hydrophobic groups and hydrogen bond acceptors (39.3 and 32.5 % of the total of all pharmacophores features of the 718 SBLB pharmacophores), but aromatic rings and hydrogen bond donors represented also an important part of the pharmacophore features (14.0 and 9.2 % respectively) far ahead negative and positive ionisable area (Fig. 7). These proportions were similar when agonist and antagonist data sets were considered separately (Fig. 7). However, when comparing the SBLB agonist and antagonist pharmacophores for each NR (Fig. 8), some significant differences (p-value <0.05) appeared in the pharmacophore features distribution (Additional file 1: Figure S3A−C). Thus, for respectively 9, 5, 4, 2 and 1 NRs, the SBLB agonist selective pharmacophores included significantly less HBA, hydrophobic, AR, PI and NI features than the corresponding SBLB antagonist selective pharmacophores. Similarly, for 3NRs, the SBLB antagonist selective pharmacophores included significantly less HBD features than the SBLB agonist selective pharmacophores. Each pharmacophore allowed to retrieve from 1 to 1299 ligands, with an average value of 32 ligands retrieved per pharmacophore (Additional file 1: Figure S2A−N, Additional file 1: Tables S1−S54).
To evaluate the pharmacophores selectivity for their dedicated NR ligands, each “SBLB agonist selective pharmacophores” and “SBLB antagonist selective pharmacophores” combinations were screened against all the other NRLiSt BDB datasets of ligands. The corresponding recalls are displayed in Fig. 9. The average recall of this large scale cross-screening was of 19.8 %. The “SBLB agonist selective pharmacophores” were associated with higher recalls with an average value of 28.8 versus 10.8 % for the “SBLB antagonist selective pharmacophores”. The most selective combination of pharmacophores was the PPAR_beta “SBLB antagonist selective pharmacophores” with an average recall of 0.001 %, and the less selective pharmacophores were the PPAR_gamma “SBLB agonist selective pharmacophores” with an average recall of 76 %. For 29 combinations of pharmacophores, the average recall was below 10 %. For only 8 combinations of pharmacophores, the average recall was above 50 %. This selectivity was significantly correlated with the number of ligands in the dataset that was used to generate the pharmacophores (Kendall’s tau coefficient, p-value = 3.476e−8, Additional file 1: Figure S4) and with the number of pharmacophores included in the combination for the considered dataset (Kendall’s tau coefficient, p-value = 5.915e−5, Additional file 1: Figure S5). The selectivity could also be correlated with the active ligands over decoys ratio (Kendall’s tau coefficient, p-value = 4.461e-11, Additional file 1: Figure S6).
Structure-based pharmacophore modeling
In the SB pharmacophore modeling approach, pharmacophores are intuitively derived from the analysis of experimentally determined (X-ray or NMR) target-ligand complexes . The identified pharmacophore features represent chemical features directly involved in the ligand-binding site interactions .
The PDB structures included in the NRLiSt BDB were used to generate SB pharmacophores. RXR_gamma was excluded of this part of the study due to the absence of a holo PDB structure. For the remaining 26 NRs, we were able to create pharmacophores that were selective for agonist ligands and pharmacophores that were selective for antagonist ligands. However this selectivity was not always achieved since the generated SB pharmacophores only covered a small part of the NRLiSt BDB ligands chemical space. Indeed, the recalls for “SB agonist selective pharmacophores” and “SB antagonist selective pharmacophores” were respectively of 55 and 8 % and the mean MCC values were 0.484 and 0.326, both varied greatly with the datasets. Very poor performance (recall of 0 % for 1 agonist dataset and 17 antagonist datasets), to very good performance (recall superior or equal to 90 % for 5 agonist datasets and MCC value superior to 0.8 for 3 agonist dataset) were obtained by screening the SB pharmacophores against the corresponding NRLiSt BDB datasets.
These performance variations were highly linked to the availability of PDB structures and to the structural diversity of the ligands that were co-cristallized in these structures. In particular, no antagonist-bound structure was available for 15 out of the 17 NRs for which no SB antagonist selective pharmacophore could be generated. For the 2 remaining NRs, PPAR_alpha and PPAR_gamma, one antagonist-bound structure was available but the SB pharmacophores created with these structures were non selective for antagonist ligands. In the same way, two agonist-bound structures were available for SF1 but no selective pharmacophore could be obtained. Surprisingly, for ROR_alpha, no antagonist-bound structure was available, but we were able to obtain an antagonist selective SB pharmacophore. Similarly, the higher recall was associated with the PPAR_gamma agonist dataset, the larger NRLiSt BDB dataset in terms of number of PDB complexes with 84 structures available.
Another interesting point was that SB pharmacophores obtained from similar NR–ligand complexes of different PDB structures could be different. For example, 11 RXR_alpha PDB structures co-crystallized with 9-cis-retinoic acid were included in the NRLiSt BDB and their corresponding 11 SB pharmacophores differed in the composition and the distribution of pharmacophore features and in the resulting hits (Fig. 10).
Ligand-based pharmacophore modeling
The LB pharmacophore modeling identifies the maximum common set of chemical features of an ensemble of ligands supposed to bind to the same active site . The pharmacophore features are presumed to be essential for the ligand-binding site interactions but no structural experimental data of these interactions are used .
The NRLiSt BDB datasets of ligands were too large and too structurally dissimilar to be represented by one unique pharmacophore. Thus, we performed a first step of ligand clustering using the LigandScout Ligand-set Clustering tool with default settings except for the cluster distance threshold that was manually adapted for each NRLiSt BDB dataset. The distance, by default set to 0.4, was lowered for some datasets to obtain a homogeneous ligands distribution. The number of resulting clusters varied for each dataset from 1 to 65. Some of the NRLiSt BDB ligands were structurally dissimilar to all the other ligands of their respective dataset and stood alone in their own cluster. Since in the LB approach two or more molecules are necessary to create a pharmacophore, whenever possible, these singletons were attributed to other clusters.
3D ligand-based pharmacophores
For each cluster, a LB pharmacophore was generated and evaluated for its selectivity. Conversely to the outcomes obtained with the SB approach, pharmacophores that were selective for agonist ligands and pharmacophores that were selective for antagonist ligand were obtained for all of the 27 NRs of the NRLiSt BDB. Since the ligand-based approach does not depend on the availability and the diversity of PDB structures, the recalls and the MCC values associated with the 3D LB pharmacophores were superior to those achieved with the 3D SB pharmacophores for all datasets but the PPAR_gamma agonist dataset. Conversely to the SB approach again, using the 3D LB pharmacophores, higher recalls and MCC values were obtained with the antagonist datasets. In particular, a recall of 100 % and a MCC value of 1 were achieved for all antagonist datasets but 5 whereas only 8 agonist datasets were fully covered by the pharmacophores that were generated.
The LB pharmacophores included exclusion volume spheres that were automatically generated around the best alignment of the ligands used for the creation of the pharmacophore. In the SB pharmacophore modeling approach, the exclusion volume spheres actually represent the shape of the active site whereas in the LB approach, the exclusion volume spheres sterically limit the pharmacophore. For some LB pharmacophores, we manually added exclusion volume spheres to remove decoys compounds using the hypothesis that if inactive compounds map all the require features, their inactivity can be due to steric clashes with the binding site [35, 36].
Combination of structure-based and ligand-based pharmacophores
Virtual ligand screening protocols often combine LB and SB approaches since including all possible information for a target enhance the chance to find hits . Several studies succeeded in identifying active compounds using hierarchical or parallel association of pharmacophore modeling with either LB and SB approaches [28, 41–49] and a few studies associated SB and LB pharmacophores [50–53]. Combinatorial use of multiple pharmacophore models for a target allowed to cover a larger chemical space of actives compared to a single model [38, 54], and can be used to enhance the hit identification success rate [43, 55–57]. In this study, we decided to combine SB and LB approaches by using one VLS method, the LigandScout pharmacophore modeling tool. The SB and LB selective pharmacophores previously generated were gathered in two pharmacophore ensembles: “SBLB agonist selective pharmacophores” and “SBLB antagonist selective pharmacophores”.
3D SBLB pharmacophores performance
The performance obtained by using a combination of SB and LB pharmacophores was largely superior to the one obtained with the SB approach alone and slightly better than the performance obtained with the LB approach alone since the LB pharmacophores were associated with extremely high recalls and MCC values (close to 100 and 1 % respectively). The recalls and MCC values of SBLB selective pharmacophores didn’t reach 100 and 1 % respectively for only 8 datasets out of the 54 that were used for this study. For two of these datasets (the LXR_alpha agonist and LXR_beta agonist datasets), one ligand could only be represented by a pharmacophore formed by two non-independent pharmacophore features (Additional file 1: Figure S7). This pharmacophore could not be used for virtual screening since no defined alignment may be found for a pharmacophore presenting less than three independent features. Thus, the ligand could not be retrieved by any of the LXR_alpha and LXR_beta agonist selective pharmacophores, which explains the resulting recall of 99.6 %. For the 6 remaining datasets, some active ligands could not be separated from decoys because of their high structural similarity (Additional file 1: Figures S8−S13).
To discriminate agonist and antagonist ligands of the NRs and cover the chemical space of the NRLiSt BDB datasets, 718 pharmacophores were necessary, with a mean value of 13 pharmacophores per dataset. The number of pharmacophores per data set was significantly correlated with the number of ligands in the respective datasets (Kendall’s tau coefficient, p-value = 9.55e−15). Since the agonist datasets are in majority larger than the antagonist datasets, the average number of pharmacophores per dataset was higher for agonist datasets compared to antagonist datasets (mean values of 15 and 11 pharmacophores per dataset respectively). Consequently, in a large majority, the number of generated selective SBLB pharmacophores was superior for agonist datasets and 4 out the 5 NRs for which the antagonist SBLB pharmacophores outnumbered the agonist ones presented datasets with more antagonist ligands than agonist ligands.
These pharmacophores were composed of 3–16 features, their majority being composed of 4 or 5 features. This finding is of great interest since it shows that pharmacophores with a limited number of features can be used to discriminate agonist and antagonist ligands of the NRs. However, it is important to note that the count of pharmacophore features did not include the number of exclusion volume spheres. The most represented pharmacophore features were hydrophobic groups and hydrogen bond acceptors, and the less represented pharmacophore features were negative and positive ionizable areas. However, some significant differences between SBLB agonist selective pharmacophores and SBLB antagonist selective pharmacophores could be observed, the main trend being that for 9 NRs, the SBLB agonist selective pharmacophores included less HBA features than the SBLB antagonist selective pharmacophores. Interestingly, some overlap exists between SB and LB pharmacophores. For example, the ER_alpha agonist selective pharmacophores include 5 SB pharmacophores and 34 LB pharmacophores. The 5 SB pharmacophores could be aligned with 5, 19, 20, 23 and 25 LB pharmacophores, with a mean of respectively 90, 78, 64, 89 and 87 % of SB features overlapped by LB features (Additional file 1: Table S55; Additional file 1: Figure S14). Similarly, the 6 SB ER_alpha antagonist selective pharmacophores could be aligned with respectively 10, 10, 10, 12, 13 and 15 out of the 17 LB ER_alpha antagonist selective pharmacophores with a mean of respectively 55, 75, 75, 78, 72, 54 and 77 % of SB features overlapped by LB features (Additional file 1: Table S56). For 7 % of the pharmacophores, one or more pharmacophore features were set as optional, which means the ligands mapping all pharmacophore features and the ligands mapping all pharmacophore features but the optional one were considered as hits.
These pharmacophores are able to retrieve a wide range of ligands, from 1 ligand for the most stringent pharmacophores to 1299 for the most powerful one. We could not identify any significant correlation between the number of pharmacophore features and the number of ligands retrieved by each pharmacophore (Additional file 1: Figure S15). Hence, the pharmacophores presenting a small number of features didn’t necessarily retrieve more ligands than the pharmacophores with a larger number of features.
Each SBLB combination of agonist selective pharmacophores and antagonist selective pharmacophores was tested for its selectivity for their dedicated NRs ligands on all NRLiSt BDB datasets. We observed that the pharmacophores generated for this study were selective for the NRs activity for which they have been created, the “SBLB antagonist selective pharmacophores” being more selective than the “SBLB antagonist selective pharmacophores”. Particularly, 6 out of the 8 combinations of pharmacophores for which the average recall of the cross-screening study was above 50 % were “SBLB agonist selective pharmacophores”. In addition, 16 out of the 17 combinations of pharmacophores for which the average recall of the cross-screening study was below 5 %, were “SBLB antagonist selective pharmacophores”. This selectivity is associated to three features: the number of ligands in the dataset used to generate the pharmacophores (Kendall’s tau coefficient, p-value = 3.476e−8), the number of pharmacophores included in the combination for the considered dataset (Kendall’s tau coefficient, p-value = 5.915e−5) and the active ligands over decoys ratio (Kendall’s tau coefficient, p-value = 4.461e−11). Hence, the selectivity of the combination of pharmacophores against its dedicated NRs activity decreased for datasets with a large number of ligands (and consequently with a large number of pharmacophores included in the combination) or with a number of active ligands largely outnumbering the number of decoys. This last point was particularly true for the PPAR_alpha, PPAR_beta and PPAR_gamma agonist datasets that encompassed more than 800 active ligands, and for which not even 10 antagonist decoys were available. A larger number of PPAR antagonist ligands would be necessary to obtain a pharmacophore that would be more selective for their agonist ligands.
Finally, the selectivity of the combination of pharmacophores also depended on the selectivity of the NRs ligands. Indeed, NRs ligands presented cross-reactivity, which means that one ligand can bind to several NRs, and this cross-reactivity is particularly true between NRs isoforms. This cross-reactivity is reflected in the NRLiSt BDB (NRLiSt BDB Supplementary Information ) and also in the pharmacophores that were generated for this study. Indeed, the combinations of pharmacophores created for this study displayed a lack of selectivity between the different NRs isoforms (ER_alpha/ER_beta, LXR_alpha/LXR_beta, PPAR_alpha/PPAR_beta/PPAR_gamma, RAR_alpha/RAR_beta/RAR_gamma, ROR_alpha/ROR_gamma, RXR_alpha/RXR_beta/RXR_gamma, TR_alpha/TR_beta), which is clearly visible with the small checkerboard along the central diagonal in the Fig. 9.
In the present work, we aimed to create NRs agonist selective pharmacophores and NRs antagonist selective pharmacophores. Our main objective was to evaluate the use of a 3D pharmacophore modeling approach to discriminate agonist and antagonist ligands of the NRs. We generated with the LigandScout software 3D structure-based pharmacophores and 3D ligand-based pharmacophores using respectively the NRs PDB structures and the sets of NRs ligands included in the NRLiSt BDB. We evaluated and compared the performance of these two approaches by focusing on two features: (1) the ability to generate selective pharmacophores for agonist or antagonist ligands, (2) the recalls associated to the combinations of pharmacophores in order to study the coverage the chemical space of the NRLiSt BDB datasets. Using the structure-based approach, we obtained selective pharmacophores for both agonist and antagonist ligands of the NRs, but since the recalls were correlated with the availability and the diversity of PDB structures, no selective pharmacophore could be generated for some datasets. Using the ligand-based approach, we created pharmacophores selective for agonist ligands and pharmacophores selective for antagonist ligands for each NRLiSt BDB datasets that yielded high performances. However, the best performances were obtained by combining the structure-based and the ligand-based approaches. We identified that: (1) the number of pharmacophores necessary to cover each NRLiSt BDB dataset depended on the number of ligands included in the dataset and (2) a limited number of pharmacophore features, mostly 4 or 5, were sufficient to discriminate agonist and antagonist ligands of the NRs. We also study the selectivity of the combination of the pharmacophores obtained for each dataset against all other datasets and we demonstrated that our pharmacophores were selective for their dedicated NRs ligands, and that this selectivity was associated with three features: the number of ligands in the dataset, the number of pharmacophores in the combination and the active ligands over decoys ratio. In conclusion, we have been able to generate 3D agonist and antagonist selective pharmacophores that cover most of the NRLiSt BDB active ligands chemical space. These 3D pharmacophores can be used as a predictor of the pharmacological activity of NRs ligands.
constitutive androstane receptor
directory of useful decoys enhanced
estrogen receptor α
estrogen receptor β
estrogen related receptor α
farnesoid X receptor α
liver X receptor α
liver X receptor β
Matthew’s correlation coefficient
- NRLiSt BDB:
nuclear receptors ligands and structures benchmarking database
Protein Data Bank
peroxisome proliferator activated receptor α
peroxisome proliferator activated receptor β
peroxisome proliferator activated receptor γ
pregnane X receptor
retinoic acid receptor α
retinoic acid receptor β
retinoic acid receptor γ
retinoic acid receptor related orphan receptor α
retinoic acid receptor related orphan receptor γ
retinoid X receptor α
retinoid X receptor β
retinoid X receptor γ
steroidogenic factor 1
thyroid hormone receptor α
thyroid hormone receptor β
vitamin D receptor
Sladek FM (2003) Nuclear receptors as drug targets: new developments in coregulators, orphan receptors and major therapeutic areas. Expert Opin Ther Targets 7(5):679–684
Sladek FM (2011) What are nuclear receptor ligands? Mol Cell Endocrinol 334(1–2):3–13
Chen T (2008) Nuclear receptor drug discovery. Curr Opin Chem Biol 12(4):418–426
Schapira M, Abagyan R, Totrov M (2003) Nuclear hormone receptor targeted virtual screening. J Med Chem 46(14):3045–3059
Mestres J (2002) Virtual screening: a real screening complement to high-throughput screening. Biochem Soc Trans 30(4):797–799
Oprea TI, Matter H (2004) Integrating virtual screening in lead discovery. Curr Opin Chem Biol 8(4):349–358
Kumar A, Zhang KY (2015) Hierarchical virtual screening approaches in small molecule drug discovery. Methods 71:26–37
Schmieder P, Koleva Y, Mekenyan O (2002) A reactivity pattern for discrimination of ER agonism and antagonism based on 3-D molecular attributes. SAR QSAR Environ Res 13(2):353–364
Bisson WH, Cheltsov AV, Bruey-Sedano N, Lin B, Chen J, Goldberger N et al (2007) Discovery of antiandrogen activity of nonsteroidal scaffolds of marketed drugs. Proc Natl Acad Sci USA 104(29):11927–11932
Nose T, Tokunaga T, Shimohigashi Y (2009) Exploration of endocrine-disrupting chemicals on estrogen receptor alpha by the agonist/antagonist differential-docking screening (AADS) method: 4-(1-adamantyl)phenol as a potent endocrine disruptor candidate. Toxicol Lett 191(1):33–39
Carrieri A, Giudici M, Parente M, De Rosas M, Piemontese L, Fracchiolla G et al (2013) Molecular determinants for nuclear receptors selectivity: chemometric analysis, dockings and site-directed mutagenesis of dual peroxisome proliferator-activated receptors alpha/gamma agonists. Eur J Med Chem 63:321–332
Zhang L, Sedykh A, Tripathi A, Zhu H, Afantitis A, Mouchlis VD et al (2013) Identification of putative estrogen receptor-mediated endocrine disrupting chemicals using QSAR- and structure-based virtual screening approaches. Toxicol Appl Pharmacol 272(1):67–76
Ekins S, Goldsmith M-R, Simon A, Zsoldos Z, Ravitz O, Williams AJ (2013) LASSO-ing potential nuclear receptor agonists and antagonists: a new computational method for database screening. J Comput Med 2013:8. doi:10.1155/2013/513537
Kolsek K, Mavri J, Sollner Dolenc M, Gobec S, S Turk (2014) Endocrine disruptome—an open source prediction tool for assessing endocrine disruption potential through nuclear receptor binding. J Chem Inf Model 54(4):1254–1267
Politi R, Rusyn I, Tropsha A (2014) Prediction of binding affinity and efficacy of thyroid hormone receptor ligands using QSAR and structure-based modeling methods. Toxicol Appl Pharmacol 280(1):177–189
Lagarde N, Zagury JF, Montes M (2014) Importance of the pharmacological profile of the bound ligand in enrichment on nuclear receptors: toward the use of experimentally validated decoy ligands. J Chem Inf Model 54(10):2915–2944
Ng HW, Zhang W, Shu M, Luo H, Ge W, Perkins R et al (2014) Competitive molecular docking approach for predicting estrogen receptor subtype alpha agonists and antagonists. BMC Bioinform 15(Suppl 11):S4
Brzozowski AM, Pike AC, Dauter Z, Hubbard RE, Bonn T, Engstrom O et al (1997) Molecular basis of agonism and antagonism in the oestrogen receptor. Nature 389(6652):753–758
Moras D, Gronemeyer H (1998) The nuclear receptor ligand-binding domain: structure and function. Curr Opin Cell Biol 10(3):384–391
Weatherman RV, Fletterick RJ, Scanlan TS (1999) Nuclear-receptor ligands and ligand-binding domains. Annu Rev Biochem 68:559–581
Bourguet W, Germain P, Gronemeyer H (2000) Nuclear receptor ligand-binding domains: three-dimensional structures, molecular interactions and pharmacological implications. Trends Pharmacol Sci 21(10):381–388
Bachmair F, Hoffmann R, Daxenbichler G, Langer T (2000) Studies on structure–activity relationships of retinoic acid receptor ligands by means of molecular modeling. Vitam Horm 59:159–215
Spencer TA, Li D, Russel JS, Collins JL, Bledsoe RK, Consler TG et al (2001) Pharmacophore analysis of the nuclear oxysterol receptor LXRalpha. J Med Chem 44(6):886–897
Schuster D, Langer T (2005) The identification of ligand features essential for PXR activation by pharmacophore modeling. J Chem Inf Model 45(2):431–439
Zhao W, Gu Q, Wang L, Ge H, Li J, Xu J (2011) Three-dimensional pharmacophore modeling of liver-X receptor agonists. J Chem Inf Model 51(9):2147–2155
Grienke U, Mihaly-Bison J, Schuster D, Afonyushkin T, Binder M, Guan SH et al (2011) Pharmacophore-based discovery of FXR-agonists. Part II: identification of bioactive triterpenes from Ganoderma lucidum. Bioorg Med Chem 19(22):6779–6791
von Grafenstein S, Mihaly-Bison J, Wolber G, Bochkov VN, Liedl KR, Schuster D (2012) Identification of novel liver X receptor activators by structure-based modeling. J Chem Inf Model 52(5):1391–1400
Temml V, Voss CV, Dirsch VM, Schuster D (2014) Discovery of new liver X receptor agonists by pharmacophore modeling and shape-based virtual screening. J Chem Inf Model 54(2):367–371
Teske K, Nandhikonda P, Bogart JW, Feleke B, Sidhu P, Yuan N et al (2014) Identification of Vdr antagonists among nuclear receptor ligands using virtual screening. Nucl Recept Res. doi:10.11131/2014/101076
Lewis SN, Garcia Z, Hontecillas R, Bassaganya-Riera J, Bevan DR (2015) Pharmacophore modeling improves virtual screening for novel peroxisome proliferator-activated receptor-gamma ligands. J Comput Aided Mol Des 29(5):421–439
Langer T, Wolber G (2004) Pharmacophore definition and 3D searches. Drug Discov Today Technol 1(3):203–207
Seidel T, Ibis G, Bendix F, Wolber G (2010) Strategies for 3D pharmacophore-based virtual screening. Drug Discov Today Technol 7(4):e221–e228
Lagarde N, Ben Nasr N, Jeremie A, Guillemain H, Laville V, Labib T et al (2014) NRLiSt BDB, the manually curated nuclear receptors ligands and structures benchmarking database. J Med Chem 57(7):3117–3125
Wolber G, Langer T (2005) LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters. J Chem Inf Model 45(1):160–169
Manetti F, Botta M, Tafi A (2006) Application of pharmacophore models in medicinal chemistry. In: Langer T, Hoffmann RD (eds) Pharmacophores and pharmacophore searches, vol 32. WILEY-VCH, Weinheim, pp 253–282
Laggner C, Wolber G, Kirchmair J, Schuster D, Langer T (2008) Pharmacophore-based virtual screening in drug discovery. In: Varnek A, Tropsha A (eds) Chemoinformatics approaches to virtual screening. Royal Society of Chemistry, London, pp 76–119
Wolber G, Dornhofer AA, Langer T (2006) Efficient overlay of small organic molecules using 3D pharmacophores. J Comput Aided Mol Des 20(12):773–788
Vuorinen A, Schuster D (2015) Methods for generating and applying pharmacophore models as virtual screening filters and for bioactivity profiling. Methods 71:113–134
Vuorinen A, Nashev LG, Odermatt A, Rollinger JM, Schuster D (2014) Pharmacophore model refinement for 11β-hydroxysteroid dehydrogenase inhibitors: search for modulators of intracellular glucocorticoid concentrations. Mol Inform 33:15–25
Kaserer T, Beck KR, Akram M, Odermatt A, Schuster D (2015) Pharmacophore models and pharmacophore-based virtual screening: concepts and applications exemplified on hydroxysteroid dehydrogenases. Molecules 20(12):22799–22832
Kumar A, Chaturvedi V, Bhatnagar S, Sinha S, Siddiqi MI (2009) Knowledge based identification of potent antitubercular compounds using structure based virtual screening and structure interaction fingerprints. J Chem Inf Model 49(1):35–42
Liu X, Xie H, Luo C, Tong L, Wang Y, Peng T et al (2010) Discovery and SAR of thiazolidine-2,4-dione analogues as insulin-like growth factor-1 receptor (IGF-1R) inhibitors via hierarchical virtual screening. J Med Chem 53(6):2661–2665
Chen Z, Tian G, Wang Z, Jiang H, Shen J, Zhu W (2010) Multiple pharmacophore models combined with molecular docking: a reliable way for efficiently identifying novel PDE4 inhibitors with high structural diversity. J Chem Inf Model 50(4):615–625
Di-wu L, Li LL, Wang WJ, Xie HZ, Yang J, Zhang CH et al (2012) Identification of CK2 inhibitors with new scaffolds by a hybrid virtual screening approach based on Bayesian model; pharmacophore hypothesis and molecular docking. J Mol Graph Model 36:42–47
Gabrielsen M, Kurczab R, Siwek A, Wolak M, Ravna AW, Kristiansen K et al (2014) Identification of novel serotonin transporter compounds by virtual screening. J Chem Inf Model 54(3):933–943
Nicolaes GA, Kulharia M, Voorberg J, Kaijen PH, Wroblewska A, Wielders S et al (2014) Rational design of small molecules targeting the C2 domain of coagulation factor VIII. Blood 123(1):113–120
Niu MM, Qin JY, Tian CP, Yan XF, Dong FG, Cheng ZQ et al (2014) Tubulin inhibitors: pharmacophore modeling, virtual screening and molecular docking. Acta Pharmacol Sin 35(7):967–979
Wang Q, Park J, Devkota AK, Cho EJ, Dalby KN, Ren P (2014) Identification and validation of novel PERK inhibitors. J Chem Inf Model 54(5):1467–1475
Kaserer T, Rigo R, Schuster P, Alcaro S, Sissi C, Schuster D (2016) Optimized virtual screening workflow for the identification of novel G-quadruplex ligands. J Chem Inf Model 56(3):484–500
Vitale RM, Gatti M, Carbone M, Barbieri F, Felicita V, Gavagnin M et al (2013) Minimalist hybrid ligand/receptor-based pharmacophore model for CXCR4 applied to a small-library of marine natural products led to the identification of phidianidine a as a new CXCR4 ligand exhibiting antagonist activity. ACS Chem Biol 8(12):2762–2770
Gangwal RP, Das NR, Thanki K, Damre MV, Dhoke GV, Sharma SS et al (2014) Identification of p38alpha MAP kinase inhibitors by pharmacophore based virtual screening. J Mol Graph Model 49:18–24
Ekins S, Freundlich JS, Coffee M (2014) A common feature pharmacophore for FDA-approved drugs inhibiting the Ebola virus. F1000Research 3:277
Pogorelcnik B, Brvar M, Zajc I, Filipic M, Solmajer T, Perdih A (2014) Monocyclic 4-amino-6-(phenylamino)-1,3,5-triazines as inhibitors of human DNA topoisomerase IIalpha. Bioorg Med Chem Lett 24(24):5762–5768
Spitzer GM, Heiss M, Mangold M, Markt P, Kirchmair J, Wolber G et al (2010) One concept, three implementations of 3D pharmacophore-based virtual screening: distinct coverage of chemical search space. J Chem Inf Model 50(7):1241–1247
Sanders MP, Barbosa AJ, Zarzycka B, Nicolaes GA, Klomp JP, de Vlieg J et al (2012) Comparative analysis of pharmacophore screening tools. J Chem Inf Model 52(6):1607–1620
Temml V, Kaserer T, Kutil Z, Landa P, Vanek T, Schuster D (2014) Pharmacophore modeling for COX-1 and -2 inhibitors with LigandScout in comparison to Discovery Studio. Future Med Chem 6(17):1869–1881
Warszycki D, Mordalski S, Kristiansen K, Kafel R, Sylte I, Chilmonczyk Z et al (2013) A linear combination of pharmacophore hypotheses as a new tool in search of new active compounds—an application for 5-HT1A receptor ligands. PLoS One 8(12):e84510
NL designed the protocol, applied it to all datasets and analyzed the results; SD applied the methods to the ER_alpha dataset; JFZ participated in the scientific discussions; MM designed and directed the work; NL, JFZ and MM wrote the manuscript. All authors read and approved the final manuscript.
We thank Prof Gerhard Wolber, Prof Thierry Langer and inte:ligand company for generously providing the LigandScout software.
The authors declare that they have no competing interests.
All authors belong to the GBA laboratory of the Conservatoire National des Arts et Métiers, Paris, France. NL is postdoctoral researcher. SD is graduate student. JFZ is full professor, chair of bioinformatics. MM is assistant professor.