Drug repositioning to propose alternative modulators for glucocorticoid receptor through structure-based virtual screening

Abstract Drug repositioning has recently become one of the widely used drug design approaches in proposing alternative compounds with potentially fewer side effects. In this study, structure-based pharmacophore modelling and docking was used to screen existing drug molecules to bring forward potential modulators for ligand-binding domain of human glucocorticoid receptor (hGR). There exist several drug molecules targeting hGR, yet their apparent side effects still persist. Our goal was to disclose new compounds via screening existing drug compounds to bring forward fast and explicit solutions. The so-called shared pharmacophore model was created using the most persistent pharmacophore features shared by several crystal structures of the receptor. The shared model was first used to screen a small database of 75 agonists and 300 antagonists/decoys, and exhibited a successful outcome in its ability to distinguish agonists from antagonists/decoys. Then, it was used to screen a database of over 5000 molecules composed of FDA-approved, worldwide used and investigational drug compounds. A total of 110 compounds satisfying the pharmacophore requirements were subjected to different docking experiments for further assessment of their binding ability. In the final hit list of 54 compounds which fulfilled all scoring criteria, 19 of them were nonsteroidal and when further investigated, each presented a unique scaffold with little structural resemblance to any known nonsteroidal GR modulators. Independent 100 ns long MD simulations conducted on three selected drug candidates in complex with hGR displayed stable conformations incorporating several hydrogen bonds common to all three compounds and the reference molecule dexamethasone. Communicated by Ramaswamy H. Sarma


Introduction
As a steroid hormone activated transcription factor, glucocorticoid receptor (GR) belongs to the superfamily of ligandregulated nuclear receptors (NRs) comprising the estrogen receptor (ER), progesterone receptor (PR), androgen receptor (AR), mineralocorticoid receptor (MR), peroxisome proliferation receptor, vitamin D receptor and the thyroid hormone receptor. It has an important role for inflammation response, and is involved in the control of cell growth, apoptosis, differentiation and metabolism (Bereshchenko et al., 2018;Kadmiel & Cidlowski, 2013). It is a ubiquitously expressed intracellular protein that regulates gene transcription. Following the first X-ray crystal structure resolved by Bledsoe et al. in 2002(Bledsoe et al., 2002, the GR ligand binding domain has become the target of all rational drug-design studies conducted in the last twenty years towards finding selective agonists that possess the desired anti-inflammatory and immunosuppressive properties with fewer side effects. Expressed by NR3C1 gene which is located on chromosome 5 (5q31), GRs have a modular structure which consists of several domains listed in sequential order; N-terminal transactivation domain (NTD, 1 to 421 residues), DNA binding domain (DBD, 422 to 486 residues), a hinge region (486 to 528 residues), ligand binding domain (LBD, 528 to 777 residues) and a C-terminal domain (Oakley & Cidlowski, 2011). The receptor consists mainly of alpha helices and a few beta sheets that fold into a three-layered alpha-helical domain which creates a fully occluded ligand-binding pocket surrounded by four helices. As a result of a ligand binding which triggers a conformational rearrangement in LBD, the receptor is released from the heteromeric complex which consists of several chaperone proteins including heat-shock protein HSP90, and enters the nucleus where it regulates gene expression directly or indirectly, through interacting with site-specific DNA sequences and several other coregulatory proteins and transcription factors.
Glocucorticoids (GCs) are the natural steroid ligands of GRs. Since 1950s, they have been widely used in the treatment of inflammatory and autoimmune diseases. However, their desired anti-inflammatory and immunosuppressive effects were often accompanied by adverse side effects such as diabetes mellitus, glaucoma, hypertension, Cushing's syndrome, osteoporosis and many others (Schake et al., 2002). GCs in complex with GRs regulate gene transcription by both transactivation and transrepression mechanisms. Many detrimental steroidal side effects are caused by transactivation whereas transrepression is considered as the key mechanism for the antiinflamatory activity. Over the years, there have been several attempts to identify ligands which preferentially induce transrepression with little or no transactivating mechanism. One such study by Schacke et al. more than a decade ago introduced a nonsteroidal selective GR agonist (SEGRA) called ZK216348 which showed a markedly superior side effect profile in addition to an improved therapeutic index comparable to classical GCs (Sch€ acke et al., 2004). However, ZK216348 displayed a poor selectivity, binding to progesterone receptor (PR), androgen receptor (AR) and mineralocorticoid receptor (MR) in addition to GR. In later years, the same research group redesigned ZK216348 to increase the selectivity towards GR by replacing the methylbenzoxazine with a quinoline moiety and came up with a prototype with reduced activity in a transactivation assay (Jaroch et al., 2010). In fact, the quinoline core was presented as a novel scaffold in an early study by Coghlan et al. which reported the first nonsteroidal GR-selective ligands with functional profiles similar to those of natural GCs (Coghlan et al., 2001). In the last two decades, several other nonsteroidal ligands such as ZK245186, CpdA, LGD5552, Org214007-0, Mapracorat, Dagrocorat and indazole ethers were presented as promising GR modulators which strongly binds GR and represses inflammatory gene expression both in vitro and in vivo studies (Baiula & Spampinato, 2014;De Bosscher et al., 2005;Hemmerling et al., 2017;L opez et al., 2008;Ripp et al., 2018;Sch€ acke et al., 2009;van Lierop et al., 2012). They all displayed a lower side effect compared to classical GCs.
Despite the shift towards nonsteroidal ligands due to their dissociative profile in transactivation versus transrepression, a successful attempt made by Uings et al. identified a novel steroidal GR ligand GW870086 which preserved its strong antiinflammatory effects while activating only a subset of genes that are normally activated by classical GCs (Uings et al., 2013). Thus, it was represented as a new topical steroid with a different safety profile to existing therapies. In addition, two triterpenoids, protopanaxadiol and protopanaxatriol, which share structural similarities with GCs were recently identified as potent and selective GR modulators which induced transrepressional activities with no or limited transactivation (Karra et al., 2019).
Despite several published X-ray crystal structures of GR complexed with either agonists or antagonists Biggadike et al., 2008;Carson et al., 2014;Edman et al., 2014;He et al., 2014;Hemmerling et al., 2017;Madauss et al., 2008;Ripa et al., 2018;Suino-Powell et al., 2008) there is a limited amount of virtual screening studies for the discovery of GR agonists. A decade ago, an application of ligand-based virtual screening was performed for the identification of novel nonsteroidal GR modulators where a total of 264,000 commercially available compounds was screened against a pharmacophore model (Onnis et al., 2010). More recently, Potatimis et al. used a structured-based virtual screening protocol using a pharmacophore model derived from a pyrrolidinone amide analogue (Potamitis et al., 2019). In this study, we present another structure-based screening protocol which used both pharmacophore and docking evaluations using a database which incorporated FDA-approved and World-approved drugs in addition to investigational compounds. Drug repositioning or "repurposing" of existing drug molecules has become a promising approach for identifying effective compounds in the last decade (Ashburn & Thor, 2004). The discovery of a new drug embraces many challenges which are encountered in different stages of its timeline, such as in vitro and in vivo screening, toxicology tests, clinical development, and many more, each requiring several years and a substantial amount of budget. Several success stories exist in drug repositioning and thus became a well-established approach in recent years. This work is the first attempt that combines pharmacophore modelling and several docking techniques for identifying potential GR agonist candidates among the existing drug compounds.
Our pharmacophore model was generated using three X-ray crystal structures of LBD domain of hGR complexed with deacylcortivazol, dexamethasone and dibC (desisobutyrylciclesonide) in the expanded binding pocket. After an extensive elimination of compounds in the database, the remaining hit molecules were subjected to a series of docking and scoring protocols. Final hit list comprised 35 steroidal compounds which all displayed GR agonist activity, and 19 nonsteroidal compounds with unique structures which presented no marked similarity with any experimentally reported nonsteroidal GR modulators.

Pharmacophore screening via LigandScout
A structure-based pharmacophore model was created for each seventeen complexes using the widely used LigandScout software tool (Wolber & Langer, 2005). Based on their similarities, seventeen pharmacophore models were clustered into groups of either two, three or four. Then, a so-called "shared" pharmacophore model that solely consists of the common pharmacophore features was generated for each group separately. In all shared models, excluded volumes representing the space which was sterically occupied by the receptor were simply discarded to prevent any bias in the results, as the receptor did not occupy a fixed space around the ligand. Furthermore, 75 agonists and 300 decoys  3cld,3e7c,3k22,3k23,3csj,4lsj,4p6w,4p6x,4udd,4udc,5g3j,5g5w,5nft,6el6,6el7,6el9. Sequence similarity indicated with red-white-blue color gradient where blue indicates the highest similarity and red the lowest similarity, b) 2 D representation of steroidal (with green box) and nonsteroidal GR (with orange box) agonists found in each complex structure.
(antagonists/unknown) extracted from ZINC (Sterling & Irwin, 2015), PubChem (Kim et al., 2019) and DUD-E decoy (Mysinger et al., 2012) databases were used to construct a small database which was then screened against each of the shared models in order to test the discriminatory power. The maximum number of pharmacophore features that can be omitted during screening was set to 1, i.e. a hit molecule was allowed to skip one feature at most. The hit compound was further evaluated using the software tool's in-house scoring function (Wolber & Langer, 2005). A receiver operation characteristic (ROC) curve was obtained to demonstrate the model's ability to distinguish agonists from antagonists. The pharmacophore model which had the highest curve above the diagonal and also the highest AUC (area under the ROC curve) value was selected to be used for screening the large database. The large database comprised drug molecules which were either approved by major jurisdictions or still under investigation for potential use. Extracted from ZINC15 database, a total of 1394 FDA-approved, 2724 Worldapproved (not FDA-approved) and 1417 investigational molecules were used for drug repurposing.

Docking experiments
Hit compounds of the pharmacophore screening with a score value above a defined threshold were selected and subjected to docking via AutoDock (Morris et al., 2009) and GOLD (Jones et al., 1997) docking tools. AutoDock tool uses a semi-empirical scoring function which combines molecular mechanics-based forcefield terms (e.g. vdW and electrostatics) and empirical terms such as solvation and ligand entropy. ChemPLP and GoldScore were the two in-house scoring functions selected for GOLD docking. In order to evaluate the accuracy of scoring functions, each of the seventeen agonists was redocked to its own receptor structure in the complex. The native state of each agonist was successfully predicted in all docking runs with satisfactorily low RMSD values below 2.0 Å.
The target protein used for the large database screening was the apo form of the crystal structure with PDB id 3bqd after removal of its agonist. It was arbitrarily chosen among seventeen receptors as they share a high degree of structural similarity ($1-2 Å of RMSD). It should be noted that 3bqd is the deacylcortivazol bound form of the receptor and the expansion of the ligand binding site as a result of deacylcortivazol allowed docking of nonsteroidal ligands with various sizes that would not be fitted into an unbound form of the receptor and thus provided an ideal template for our screening studies in the search of both steroidal and nonsteroidal GR modulators (Suino-Powell et al., 2008).

Molecular dynamics (MD) simulations of hit compounds
Three hit compounds selected from the docking experiments and the steroid dexamethasone as a reference were each subjected to 100 ns long independent MD simulations in complex with the target receptor that was used for docking (PDB id: 3bqd). The drug-receptor complex was solvated in a box of TIP3P water molecules, and ionized with Na þ2 and Clto neutralize the total charge of the system. The whole system incorporated a total of 41,006 atoms of which 36,780 belonged to water molecules. OPLS-AA/M forcefield (Jorgensen et al., 1996) was employed for describing the interaction potential of both receptor and drug molecules. The forcefield parameters for drug molecules were generated via the Web-based service LigParGen (Jorgensen & Tirado-Rives, 2005). Using NAMD v2.13 software tool (Phillips et al., 2005), a 1000 steps of energy minimization was followed by 100 ns MD simulation for which the first 5 ns was taken as equilibration phase on the basis of the root mean square deviation (RMSD) profile. All simulations were conducted at constant NPT at 310 K using Langevin dynamics for all non-hydrogen atoms, with a Langevin damping coefficient of 5 ps À1 . The system was kept at a constant pressure of 1 atm by using a Node-Hoover Langevin piston (Feller et al., 1995) with a period of 100 fs and damping timescale of 50 ps. Long-range electrostatic interactions were treated by particle mesh Ewald (PME) method, with a grid point density of 1 A˚. A cutoff of 12 A˚was used for van der Waals and short-range electrostatics interactions with a switching function. Time step was set to 2 fs by using SHAKE algorithm for bonds involving hydrogens (Ryckaert et al., 1977) and the data was recorded at every 20 ps to generate a trajectory of 5000 conformations.

Results and discussion
Generation of "shared" pharmacophore models As mentioned in Methods section, a pharmacophore model was created for each of seventeen complexes, two of which were illustrated in Figure 2 as the steroidal deacylcortivasol bound complex (PDB id: 3bqd) and the nonsteroidal indazole amide bound complex (PDB id: 3k23). The models of the remaining fifteen complexes were provided in Supplementary Figure 1. The hydrophobic features were represented with yellow spheres surrounding the atomic groups on the ligand, whereas hydrogen bond donor and acceptor groups were both represented with a sphere and an arrow. The red sphere corresponds to hydrogen bond acceptor group on the ligand and the red arrow emanating from the receptor simply points towards the acceptor group. Similarly, the green sphere represents the hydrogen bond donor group of the ligand and the green arrow emanates from the ligand towards its counterpart in the receptor. Table 1 gives the full list of the number of pharmacophore features for each complex along with ligand ids and PDB ids. Based on the similarity of their features, some of these models were grouped together to create a "shared" pharmacophore which incorporated the features shared by all its group members. The groups consisted of either two, three or four pharmacophore models as listed in Table 2 and also depicted in Figure 3. This led to a total of seven distinct shared pharmacophore models which were further used for screening a small database in order to test each model's discriminatory power.
Critical assessment of "shared" pharmacophore models for discriminatory power A small database composed of 75 agonists and 300 decoys (antagonists/unknown) was screened against each shared pharmacophore model to evaluate its capability to distinguish agonists from antagonists. Figure 4 illustrates ROC curves for each model along with its total number of hits and AUC values. Model #1 (DAY-CV71782-DEX) displayed the highest number of hits of 284 (73 agonists versus 211 decoys) with the highest enrichment factor of 5.0 and 4.8 in the top scoring 5% and top 10%, respectively (See last column in Table 2). Despite the fact that model #2 (HCY-MOF-GW6) had the highest AUC value of 0.826, none of the compounds in the hit list passed the Pharmacophore Fit score threshold of 45 which was selected as the "model exhaustion" point. On the other hand, 22 hits in model #1 had scores above 45 which yielded 100% true positives. For screening larger databases, the same threshold of 45 was used along with the selected pharmacophore model #1.

Establishing thresholds for score values in dockingbased virtual screening
Prior to screening the large database, 22 hit agonists having pharmacophore score value above 45 were subjected to docking. The goal here was to determine a threshold value for each score. For AutoDock, a threshold value of À10.0 kcal/mol was adopted such that the majority of the agonists in the list would be retained. As a result, only 2 out of 22 agonists were eliminated. Likewise, the threshold value was taken as 40 for both GoldScore and ChemPLP score values which led to a final list of 21 and 20 agonists respectively. As illustrated for the training database on the left side of the flowchart in Figure 5, all three docking results were later merged into a consensus list which consisted of compounds satisfying all three docking score thresholds. Accordingly, 18 agonists were left in the final hit list.

Final evaluation of the docked compounds in the binding pocket
In the final stage of filtering, three best poses obtained from three docking results, AutoDock, ChemPLP and GoldScore, were aligned to each other in the binding pocket in order to determine the consistency of the results. If either one of the best poses of one docking was different from another best pose by more than 4-5 Å, that compound was simply discarded from the hit list. Accordingly, none of 22 compounds showed inconsistency in their docked poses, thus were kept in the final hit list. High-throughout virtual screening through pharmacophore screening and consensus docking of three databases Three subsets of compounds extracted from ZINC database were subjected to pharmacophore screening using the selected pharmacophore model with ID #1. As indicated on the right side of the flowchart in Figure 5, nearly two thirds of the compounds in World-approved subsets were discarded at the first stage, simply because they didn't hold the required pharmacophore features. Also, both FDA-approved and investigational subsets diminished half in size at the end of the first stage of filtering. All the hit compounds satisfying the feature requirements were then evaluated based on their pharmacophore fit score value. A significant portion of compounds (> 90%) in each subset was discarded as they displayed a score value lower than 45. Once more, the highest elimination was observed for World-approved subset which lost nearly 98% of its content. As a result, 109 compounds in total from all three databases passed the second stage of filtering.
Each of 109 molecules was docked using different docking tools and scoring functions. The best pose extracted from each docking experiment was evaluated based on the corresponding predefined threshold value. Among three scoring evaluations, AutoDock yielded the least amount of compounds (75 out of 109) which satisfied the docking threshold, whereas the majority of compounds (>93%) passed the docking threshold for GoldScore and ChemPLP. Then, all the docked poses were merged to create a consensus pool which comprised a total of 67 hits satisfying all three docking threshold requirements. Among them, 27 were FDA-approved, 13 were World-approved and 27 were investigational compounds. Finally, any compound with its three best poses (one from each docking experiment) displaying a discordance in their relative orientation in the binding pocket was simply discarded. Accordingly, 13 compounds were discarded and thus in the final hit list of 54 molecules, there were 22 FDA-approved, 12 World-approved and 20 investigational compounds of which their best pose was illustrated in Figure 6.
The corresponding ZINC database IDs, compound names (trade names if available) and all three score values were provided in two separate tables in Supplementary Tables 2  and 3 for each of 54 compounds. The order of the list in both tables was based on ChemPLP score from highest to lowest. Furthermore, AutoDock score was plotted against GoldScore value for each of 54 best poses as depicted in Figure 7 where a considerable amount of correlation between the two score values was observed especially for FDA-approved and investigational compounds. Same plot was created for ChemPLP versus AutoDock scores and a similar distribution profile was observed.
Additionally, all hit compounds were separately aligned to the pharmacophore model as illustrated in Figure 8a for steroids and Figure 8c for nonsteroids. It is clear that the hit compounds incorporated the most essential features of the model with satisfactory pharmacophore fit score values. Moreover, the neighbouring residues interacting with the hit compound was provided in Supplementary Table 1. Accordingly, five residues which interacted with more than 40 hit compounds were recognized as Asn564, Gly567, Met604, Phe623 and Gln642 (See grey shaded cells in the table). Among them, two polar uncharged residues, Asn564 and Gln642 interacted via hydrogen bonds with the majority of the compounds, whereas the hydrophobic residues Gly, Met and Phe were in close contact making short-range hydrophobic interactions. Furthermore, Asn564 and Gln642 were the two critical residues that interacted with all 17 ligands in the resolved complex structures which were initially used to build the pharmacophore models. The interaction diagram of one selected steroidal and one nonsteroidal compound manifesting the highest score value in its category was also provided in Figure 8b and d. The 3 D version of the interacting residues nearby the compounds was also provided in Suppl Figure 2.

Critical assessment of three selected hit compounds via MD simulations
Two steroidal investigational (IDs:12a,13a) and one nonsteroidal FDA approved (ID:17c) compound were selected from the list of 54 hits. Their best poses with the receptor were used as the initial conformation for MD simulations. In addition, the well-known FDA approved steroid dexamethasone was used as a reference molecule for the critical assessment of our results. As illustrated in Figure 9a, in all four simulations, RMSD value of the protein with respect to its initial state (t ¼ 0 ns) indicated a stabilization of the structure within the first 5 ns. During 95 ns of production phase, there was no major structural deviation of the receptor from its initial state. RMSD profile of each drug compound in the complex was depicted in Figure 9b. Accordingly, dexamethasone and one steroidal compound (12a) preserved their initial conformations to a great extent, whereas compounds 13a and 17c rapidly changed their conformations by 2-3 Å in the early stages of the trajectory, yet preserved their states throughout the simulation.
Root mean square fluctuation profile of each system was illustrated in Figure 9c alongside the binding site residues indicated with red squares. No major difference was observed in the mobility of the receptor when bound to different compounds, except a few deviations observed in the loop regions. Binding site residues well coincided with the least mobile regions of the receptor which provided a stable interaction between the receptor and the drug molecule. The change in the radius of gyration of the complex was monitored throughout the trajectory as illustrated in Figure 9d. For all four complexes, the size of system was well preserved at its initial value, roughly 18.5 Å. Finally, the hydrogen bond network between the ligand and the receptor was investigated in detail for every snapshot in the trajectory. As illustrated in Figure 9e, the total number of hydrogen bonds between the drug compound and the receptor was monitored throughout the simulation. For the reference molecule dexamethasone, half of the snapshots incorporated at least one hydrogen bond with the receptor. For the hit compounds 13a and 17c, one third of the trajectory displayed at least one hydrogen bond. On the other hand, the hit steroid compound 12a presented a total of five hydrogen bonds with the receptor in eight conformations, which was the highest number observed in all simulations. Nearly 70% of the time, compound 12a incorporated at least one hydrogen bond with the receptor, which indicated its strong binding affinity compared to the other compounds.
Moreover, neighbouring residues close to the bound ligand by less than 5 Å were determined and their frequency of occurrence over the entire trajectory was calculated. The list in Figure 9f shows a total of 24 residues with 100% occurrence in at least one of the drug-receptor simulations. 16 out of 24 residues highlighted in green were present in more than 97% of the collected 5000 snapshots. Finally, residues involved in hydrogen bonding with the ligand were indicated in red colour. Met560, Leu563, Asn564, Gln570, Gln642, Cys736 were the most critical residues interacting via hydrogen bonds with the drug compound at all times in the majority of the four complexes.  part of two of seventeen X-ray crystal complexes initially selected for constructing our screening protocol, and both appeared in our hit list with relatively high score values (See Supplementary Table 2). Moreover, GW870086 which is a potent anti-inflammatory compound identified by Uings et al. (Uings et al., 2013) shared the same steroid scaffold of model A with a unique set of side groups and listed in Table 3 for comparison. GW870086 was known with its distinctive ability to regulate only a subset of genes that are normally affected by classical GRs, thus presented itself as the most specific steroids available. Under model A, there were a total of six steroidal compounds from investigational subset (IDs: 11a, 12a, 13a, 14a, 15a and 16a), two of which were isomers (15a and 16a). Thus, our results highlighted these compounds to be further investigated for potential GR modulators with fewer side effects.
Despite the fact that the shared pharmacophore model used for screening was derived from three steroidal compounds (dexamethasone (DEX), deacylcortivasol (DAY) and desisobuytyryl ciclesonide (CV71782)), 19 out of 54 hits were found to be nonsteroidal. Among 19, two pairs were found to be isomers, thus each of 17 compounds displayed unique structures as depicted in Table 4. In fact, it was not possible to cluster them based on their similarity due to their structural distinctiveness. Moreover, their structural similarities to experimentally reported nonsteroidal glucocorticoid modulators (  2018; van den Heuvel et al., 2016;van Lierop et al., 2012) were determined using Tanimoto coefficient (Tc) similarity index via ChemMine tool (Backman et al., 2011). For each nonsteroidal hit, the GR modulator with the highest similarity index was listed alongside in the table. Overall, the highest Tanimoto index was observed as 0.35 which was well below the accepted cutoff value of 0.5 for an acceptable degree of similarity.
Sorted by ChemPLP score values in Supplementary Table 2, the majority of compounds in the highest score range was nonsteroidal and belonged to investigational subset. As listed in Supplementary Table 3, anlotinib which is a nonsteroidal tyrosine kinase inhibitor displayed the highest ChemPLP (92.87) and the lowest negative AutoDock score values (-12.90) and it is still under investigation for lung and colorectal cancer treatment. Furthermore, as illustrated in Table 4, it shares the highest similarity (Tc ¼ 0.339) with Compound-42 which was identified as a novel, potent and selective GR partial agonist by Takahashi et al. in 2007. In the high score range, there is another tyrosine kinase inhibitor which is brivanib alaninate and also shares the highest similarity (Tc ¼ 0.285) with Compound-42. However, Tc values are well below 0.5 to Table 3. A total of 35 steroidal hit compounds categorized under two groups based on the scaffolds represented as model A and B.
Ã Ligand in one of 17 X-ray structures used: PDB id 3cld. ÃÃ Ligand in one of 17 X-ray structures used: PDB id 4udc. ÃÃÃ GW870086: GR agonist acting as a potent anti-inflammatory agent. ÃÃÃÃ 9 b and 10 b are isomers. Table 4. A total of 17 nonsteroidal hit compounds and the experimentally identified GR modulator to which it shares the highest similarity indicated by Tanimoto coefficient (Tc).
suggest any striking structural similarity. Unlike steroidal compounds which were either GR agonist or metabolite, nonsteroidal compounds in the hit list displayed different activities such as G-protein coupled receptor (GPCR) agonists, retinoic acid metabolism inhibitor, topoisomerase inhibitor, androgen receptor modulator, and an antibiotic tetracycline. The fact that none of the nonsteroidal compounds were GR agonists propose them as alternative candidates of which the majority was still under investigational stages (14 investigational versus 3 FDA-approved and 2 World-approved). Thus, each of these compounds used for the treatment of various diseases can be repurposed as GR modulators.

Conclusions
Thousands of compounds from FDA-approved, Worldapproved and investigational databases were screened using a shared pharmacophore model holding the common pharmacophore features of three known crystal structures of GR receptor in complex with their steroidal agonists. Nearly half of the compounds in each database satisfied the shared pharmacophore features, yet a significant portion of them (> 90%) was discarded after filtering with a threshold score value of 45. This threshold value was determined after a small database composed of 75 agonists and 300 decoys (antagonists/unknown) was screened. None of the antagonists with the shared pharmacophore features displayed a score value greater than 45, thus discarded, whereas the majority ($30%, 22 out of 75) of screened agonists had scores well above that threshold.
Out of 5535 compounds, a total of 110 satisfying all the pharmacophore requirements were subjected to three different docking experiments for further evaluation as potential binders. Similar to pharmacophore screening, different thresholds were determined for each docking evaluation using the small database, prior to screening the large database. As a result, the majority of the compounds satisfied the predefined thresholds, especially in Goldscore and ChemPLP scoring evaluations. Upon merging the results by extracting the compounds that displayed desired affinities towards GR in all three docking experiments, a pool of 67 compounds was obtained. Finally, three best poses of each compound collected from three docking were aligned and compared to each other, and the compounds which displayed different conformations/orientations in their best pose were simply discarded. In the final evaluation step, 100 ns long MD simulations were conducted on complex systems with three selected compounds from the hit list along with the known agonist dexamethasone used as a reference. Throughout the simulations, all three drug compounds formed stable interactions with the receptor via multiple hydrogen bonds with several residues commonly observed in dexamethasone-GR complex.
In the final hit list of 54 compounds, 35 were steroidal GR agonists. The remaining 19 compounds were nonsteroidal despite the fact that the initial shared pharmacophore model was constructed from receptor structures in complex with steroidal compounds. Each of 19 compounds displayed unique structures with a low degree of similarity to experimentally reported nonsteroidal GR modulators. Consequently, even though the experimental validation lacks, the distinctiveness of each of these compounds might pave the way for their repurposing to achieve potential agonist activities with fewer side effects, when used in appropriate dosages.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Availability of data and material
Available if required.

Code availability
Available if required.

Author contributions
The manuscript was written through the contributions of all authors. All authors have given approval to the final version of the manuscript.