Structure-based discovery of potent and selective melatonin receptor agonists

Melatonin receptors MT1 and MT2 are involved in synchronizing circadian rhythms and are important targets for treating sleep and mood disorders, type-2 diabetes and cancer. Here, we performed large scale structure-based virtual screening for new ligand chemotypes using recently solved high-resolution 3D crystal structures of agonist-bound MT receptors. Experimental testing of 62 screening candidates yielded the discovery of 10 new agonist chemotypes with sub-micromolar potency at MT receptors, with compound 21 reaching EC50 of 0.36 nM. Six of these molecules displayed selectivity for MT2 over MT1. Moreover, two most potent agonists, including 21 and a close derivative of melatonin, 28, had dramatically reduced arrestin recruitment at MT2, while compound 37 was devoid of Gi signaling at MT1, implying biased signaling. This study validates the suitability of the agonist-bound orthosteric pocket in the MT receptor structures for the structure-based discovery of selective agonists.


Introduction
The type 1A and 1B melatonin receptors (MT 1 and MT 2 ) are G protein-coupled receptors (GPCRs) that respond to the neurohormone melatonin (N-acetyl-5-methoxytryptamine) (Pévet, 2016;Reppert et al., 1994). Melatonin is found in all mammals, including humans, where it regulates sleep and helps to synchronize the circadian rhythm with natural light-dark cycles (Brzezinski, 1997;Xie et al., 2017). Chemically, melatonin is synthesized from serotonin in the pineal gland of the brain during darkness (Ganguly et al., 2002). Both MT 1 and MT 2 share canonical helical 7-transmembrane (7-TM) topology Stauch et al., 2019), although they are differentially expressed and implicated in diverse biological functions and pathologies (Dubocovich and Markowska, 2005). While exogenous melatonin has been commonly used for the treatment of insomnia and jetlag, more effective and long-lasting MT agonists such as ramelteon have been approved for primary chronic insomnia treatment, because of their low side-effect profile as compared to other sleeping aids such as benzodiazepines (Hardeland et al., 2011;Erman et al., 2006). Other MT agonists such as tasimelteon and agomelatine, are used for non-24-hour sleepwake disorders in blind individuals and as an atypical anti-depressant for major depressive disorders, respectively (Lavedan et al., 2015;de Bodinat et al., 2010). Recent studies also suggest MT receptors play an essential role in learning, memory, and neuroprotection (Liu et al., 2016) and illustrate the potential utility of partial and selective MT 2 receptor agonists as antinociceptive drugs (Ló pez-

Benchmarking receptor models
To evaluate the ability of the structure-based receptor models to recognize high-affinity melatonin receptor ligands, we performed extensive docking of a subset of known ligands of MT 1 and MT 2 receptors (Figure 1-figure supplement 1) into (i) the unmodified 3D structures obtained from X-ray crystallography (MT 1 _Xtal, MT 2 _Xtal), as well as (ii) into the receptor models where the ligandbinding pocket was optimized by conformational modeling (MT 1 _Opt, MT 2 _Opt). Analysis of the docking poses for the known MT ligands in both crystal structures and optimized MT receptor models showed favorable binding scores with docking poses consistent with the orientation and binding modes of crystallized ligands (Figure 1a-d). The major interactions include aromatic stacking of the heterocyclic core with ECL2 hydrophobic residue F179/192 ECL2 (the residue numbers for MT 1 and MT 2 listed for UniProt (Bateman et al., 2017) ids: P48039 and P49286, respectively, followed by superscripted Ballesteros -Weinstein numbering scheme Ballesteros and Weinstein, 1995), as well as hydrogen bonding interactions with N162/175 4.60 and Q181/194 ECL2 Stauch et al., 2019) . The performance of each model was then evaluated as the area under the corresponding receiver operator characteristic (ROC) curve (AUC), benchmarking the ability of these models to correctly detect ligands among decoys. The AUC values for the optimized models of MT receptors showed substantial improvement over AUC values for MT crystal structures (MT 1 _Opt = 87 vs. MT 1 _Xtal = 69; and MT 2 _Opt = 82 vs. MT 2 _Xtal = 70) (see Figure 1e). Overall, these results validated the improved VLS performance of the optimized models of MT 1 and MT 2 receptors, which were then used for large-scale prospective VLS.

Prospective VLS and candidate selection
A library of 8.4 million commercially available compounds was docked into the optimized MT 1 _Opt and MT 2 _Opt structural models (see Materials and methods), and for every compound, docking scores and binding interactions were predicted. The top 5000 scoring compounds were selected from these VLS runs for each receptor, which were further evaluated by redocking into both MT 1 and MT 2 receptors with increased computational sampling. The initial hit list contained 700 compounds predicted to bind to both receptors. To evaluate these top docking solutions, we created additional models of MT receptors by restoring thermostabilizing mutations (Figure 1-figure supplement 2) in the proximity of the orthosteric site to wild-type residues (A104 3.29 G and W251 6.48 F in MT 1 ; W264 6.48 F in MT 2 ), and performed further conformational optimizations. We determined that the impact of these mutations on the docking results was negligible. The dock scores for selected MT ligands were better than the standard ICM VLS cutoff À32.0 kJ/mol, which is better than or comparable to the docking score of melatonin (À29.3) and other high affinity MT receptor ligands .
To capture chemotype diversity, we selected the top 500 compounds for each receptor using chemical clustering in combination with docking scores. A final set of 62 compounds (23 from only MT 1 VLS; 25 from only MT 2 VLS; 14 from both MT 1 and MT 2 VLS) were selected for purchase based on a multidimensional composite criterion accounting for compound novelty, chemical diversity, well-defined interaction patterns with binding site residues N162/175 4.60 and/or Q181/194 ECL2 , and interaction similarities with ligands observed in the crystal structures (See Figure 2; Supplementary file 1 Table S1).
Most of the compounds represented new chemotypes with Tanimoto chemical distance values >0.22 (Abagyan et al., 2016), separating them from known high-affinity MT ligands available in CHEMBL24 (Gaulton et al., 2017). We also chose a close analog of melatonin -compound 28 (Tanimoto distance = 0.05), which to our knowledge, has not yet been characterized as a ligand for MT receptors (Gaulton et al., 2017;Kim et al., 2019). Compound 28 served as an additional positive control, which also helped us to evaluate the effect of a single chemical group substitution in melatonin on the binding and function at the MT receptors.

Experimental hit identification and validation
The selected 62 candidate compounds from VLS were tested experimentally for binding and functional profiles in both MT receptors. Eleven compounds (Table 1; Figure 3) demonstrated sub-   micromolar potencies in G i/o mediated signaling assays (18% hit rate). Ten of these eleven compounds also showed binding affinities K i <10 mM in a competition binding assay (Figure 3-figure supplement 1). The melatonin derivative 28 identified by VLS was as potent as melatonin itself in MT 2 (EC 50 = 0.04 nM) and had the same potency (EC 50 = 0.04 nM) at MT 1 . The most potent new chemotype, 21, displayed an EC 50 = 0.36 nM for MT 2 , with a 30-fold selectivity over MT 1 receptor (MT 1 EC 50 = 12 nM). Overall, seven hits had EC 50 <100 nM for at least one of the melatonin receptors. Similar to other low molecular weight MT ligands, most of the hits identified belong to a library of fragment-like compounds with molecular weights less than 250 Da, and have exceptionally high ligand efficiency (LE), far exceeding the~0.3 value considered as a standard for a promising lead. For example, compound 21 (Mol. Wt. = 224 Da) had the highest LE values of 0.83 kcal/mol per nonhydrogen atom for MT 2 and 0.69 kcal/mol per non-hydrogen atom for MT 1 receptors (Hopkins et al., 2004;Hopkins et al., 2014). The excellent LE of these molecules allows the potential for further optimization of their drug-like properties.

Chemical and conformational diversity of hits
Most of the hit compounds, as shown in Chemical structure 1, are novel and display diverse chemotypes distinct from known high-affinity MT ligands (ChEMBL, pAct >6.0), with Tanimoto distance exceeding 0.4 for all but two ligands (28 and 23). While the majority of known MT agonists reported in ChEMBL have either substituted indene or naphthalene core, only two of the eleven hits reported here have fused heterocycles and several others have two substituted aromatic rings connected by a flexible chain. Most compounds have diverse substitutions at positions topologically equivalent to the 5-methoxy, acetylamido and C2 position of melatonin ( Figure 2). Two compounds-21 and 29-  have a pyrimidine scaffold, whereas four compounds-23, 37, 44, and 57have a methoxyphenyl group in place of the 5-methoxy indole scaffold in melatonin. Another interesting core is the cyclopentyl-fused thienopyridine of compound 45. Only 2 compounds, 28 and 54, have substituted indoles similar to melatonin. The predicted binding poses of the selected hit compounds in their docking models of MT receptors are shown in Figure 4. Nine out of eleven hits have methoxy or a similar group predicted to make hydrogen bonding interactions with N162/175 4.60 , which was found to be a critical residue for receptor activation, despite playing a limited role in ligand affinity or structural stability of the receptor .
Furthermore, seven of the hits were predicted to form hydrogen bonding interactions with Q181/ 194 ECL2 similar to alkylamide tail of melatonin. Five hits were predicted to occupy a significant space in the pocket flanked by TMs-II, III, and VII forming hydrophobic interactions, especially with residues Y281/294 7.38 and Y285/298 7.43 , as previously found in the MT receptor structures Stauch et al., 2019). These hydrophobic interactions are similar to those formed by the phenyl moiety of 2-phenylmelatonin. Both types of hydrogen bonding and hydrophobic interactions were found to be critical for a ligand's steric fit into the MT receptor binding pocket and are the primary determinants of ligand affinity.

Structural basis of subtype selectivity of the hits
Six of the identified hits were found to be at least 30 fold more potent at MT 2 compared to MT 1 in the G i/o-mediated cAMP inhibition assays. Among the hits reported here with novel scaffolds, compound 21 has the highest potency for both MT 2 (EC 50 = 0.36 nM) and MT 1 (EC 50 = 12 nM). Compound 21 was predicted to bind both the MT receptors in a similar orientation by forming hydrogen bonding interactions with N162/175 4.60 and Q181/194 ECL2 with its methoxy anchor and acetylamido tail, respectively. These interactions had been reported to be critical for ligand affinity and potency at the MT receptors Stauch et al., 2019).
Other compounds also possess remarkable MT 2 selectivity. For example, compound 47 is 187fold selective towards MT 2 (EC 50 = 10 nM for MT 2 , and 2.34 mM for MT 1 , respectively). The pyrrole ring mimics the indole ring of melatonin, the amide group forms hydrogen bonding with N162/ Chemical structure 1. Chemical structures of hit compounds with EC 50 <1 mM at the MT receptors.
175 4.60 and the chlorophenyl group forms hydrophobic interactions with ECL2 and TM-II, III, and VII residues ( Figure 4). Despite the lack of polar interactions with Q181/194 ECL2 , the compound displays sub-micromolar potency for MT 2 . Similarly, compound 45 also lacks a substitution topologically equivalent to acetylamido tail of melatonin (R2 feature) and yet has a sub-micromolar potency and 17-fold selectivity for MT 2 (EC 50 = 427 nM). In contrast, compound 44 was predicted to form interactions with Q181/194 ECL2 , but it lacks an R3 equivalent substitution, which still makes it 267-fold selective for MT 2 (EC 50 = 263 nM). These findings suggest that either R2 or R3 could be sufficient in maintaining the potency and selectivity at MT 2 .

Functional selectivity of the hit ligands
All the discovered hits show activity as agonists in G i/o -protein signaling assays at both MT 1 and MT 2 receptors. At the same time, some compounds show functional profiles notably distinct from full and balanced agonism, especially at MT 2 . Thus, four of the hits, 28, 29, 57, and 62 had their efficacy (E max ) reduced to less than 70% in MT 2 , and are therefore considered partial agonists (Audinot et al., 2003). The identified hits were also evaluated for their b-arrestin recruitment (Figure 3-figure supplements 2 and 3), with the comparative analysis of G-protein and b-arrestin activity shown in Figure 5. In the case of the MT 1 receptor, there are no significant deviations from the overall balanced G-protein/Arrestin signaling profiles for most compounds (Figure 5a). One exception is compound 37, which completely lacks G-protein signaling, though it still binds to MT 1 and displays substantial b-arrestin activity. In the case of MT 2 , however, there are several compounds that show a marked reduction in b-arrestin signaling compared to G-protein, especially compounds 21 and 28, which show bias factors of 15.5 and 33.9, respectively (Figure 5b). These results suggest that MT ligands may show rather distinct functional bias profiles in G-protein and b-arrestin signaling, as observed in many other GPCRs (Kenakin, 2019;Roth, 2019), though the biological importance of this bias in MT remains to be investigated (Cecon et al., 2018).
To gain more insights into these variations in ligand activity at MT receptors, we analyzed conformational differences among ligand-receptor pairs for these compounds. As the hit compounds are fragment-like and may attain multiple energetically-favorable poses at the orthosteric site upon docking, the specific conformational features driving partial agonism remain unclear. However, analysis of compounds 28 and 37 with the most pronounced bias to G-protein in MT 2 and b-arrestin in MT 1 , respectively, suggests possible explanations for these phenomena.
Compound 28 is very similar to melatonin, except the amide is replaced by a urea. This substitution renders compound 28 as a partial agonist at MT 2 (E max = 69.4%) while largely maintaining full agonism at MT 1 (E max = 95.3%). Docking predictions suggest that compound 28 assumes an orientation in the binding pocket similar to 2-phenylmelatonin with subtle differences, as shown in Figure 6. In MT 2 orthosteric site, the urea of compound 28 forms hydrogen bonding interactions with the side chains of polar residues Q194 ECL2 and N268 6.52 owing to its additional nitrogen. Such interactions, however, are energetically unfavorable in the case of MT 1 with possible steric clashes Stauch et al., 2019). Instead, the interactions of acetylamido group of melatonin with Q181 ECL2 are replaced by a hydrogen bond between Y281 7.38 and oxygen from the urea in compound 28. These interactions become favorable in MT 1 as the residue Y281 7.38 is rotated towards TM-VI placing it 4 Å away from T178 ECL2 . In the case of MT 2 , however, residues Y294 7.38 -T191 ECL2 are 3 Å apart forming an intermolecular hydrogen bond with Y281 7.38 oriented away from TM-VI allowing favorable orientation of Q194 ECL2 to form a hydrogen bond with compound 28.
A similar pattern of ligand-receptor interaction is observed from the docking of the most selective compound 37 into MT receptors. Compound 37 has a distinct and much bulkier substitution with a 3-cyclopropyl-1,2,4-oxadiazol group (R2 feature). In MT 1 , this oxadiazol group is predicted to form hydrogen bonds with Q181 ECL2 , T178 ECL2 and Y281 7.38 , while the cyclopropyl group is predicted to fit in the sub pocket formed by ECL2, TM-V, and VII. These interaction pattern changes in MT 2 , where the prominent hydrogen bonding between side chains of T191 ECL2 and Y294 7.38 is formed, precluding hydrogen bonding of these residues to compound 37. Moreover, the methoxy group (R1 feature) of the compound, which maintains a hydrogen bond with N175 4.60 in MT 2 is lost with N162 4.60 side chain in MT1, due to a subtle shift of the compound. Indeed, the methoxy -N162 4.60 interaction is found to be critical in receptor activation , and loss of this interaction is likely to explain lack of activity of compound 37 in G i -mediated signaling at MT 1 . This interaction difference, however, does not seem to affect the b-arrestin mediated signaling by compound 37 in MT 1 . This peculiar feature of 37 is supported by our mutational studies, where N162Q mutation in MT 1 actually increased potency of 37 slightly (2-fold), while Y281F, as expected, reduced potency by over 10x (Figure 6-figure supplement 1). For comparison, 28 drastically (>100 fold) reduced potency in both MT 1 mutants N162Q and Y281F. Taken together, these results support a key role of N162/175 4.60 anchoring interactions in G i -mediated receptor activation, and also suggest a distinct role of residues Y281/294 7.38 in governing ligand bias at MT receptors. Further analysis, including mutation and SAR studies of compound 37 derivatives are needed for comprehensive validation of this hypothetical mechanism in future studies.

Off-target profiling
To verify the ligand selectivity, the lead compounds 21, 28, and 37, were subjected to binding profiling at a panel of 47 common drug targets (including many GPCRs and neurotransmitter transporters). At a final concentration of 10 mM, they did not show any substantial binding at these targets, except for compound 28 that displayed over 50% inhibition at three 5-HT receptors (Figure 6-figure supplement 2), with binding affinity of 851 nM (pK i = 6.07 ± 0.11) for 5-HT 1A , 1525 nM  (pK i = 5.82 ± 0.02) for 5-HT 1D , and 286 nM (pK i = 6.54 ± 0.03) for 5-HT 7A . Assessment of functional activity of 28 in G s -mediated cAMP production shows that it is a weak partial agonist with EC 50 = 1819 nM, as compared to 0.04 nM at both MT 1 and MT 2 receptors. These results demonstrated that the lead compounds 21, 28, and 37 act specifically at MT receptors and show high selectivity for MT 1 or MT 2 over many common drug targets.

Discussion
The discovery of potent and selective MT ligands with novel chemotypes holds promise for the development of next-generation drugs to treat circadian rhythm and mood disorders, pain, insomnia, type-2 diabetes, and cancer (Karamitri and Jockers, 2019;Liu et al., 2016). Herein, we utilized the recently solved 3D structures of the melatonin receptors, in complex with the agonist 2-phenylmelatonin, Stauch et al., 2019) to perform prospective virtual ligand screening of large fragment-like compound libraries. This approach resulted in the discovery of ten new chemotypes of potent agonists, both full and partial, for MT 1 and MT 2 . The number of submicromolar hits and potency of the best among them is one of the highest reported for a VLS campaign in class A GPCRs (Lyu et al., 2019) and in-line with another VLS screen for MT receptors, published while this study was in revision (Stein et al., 2020). This is remarkable, considering that most GPCR structures have a limited capacity to distinguish agonists vs. antagonists (Costanzi and Vilar, 2012;Weiss et al., 2018) and prospective VLS campaigns often result in antagonists even when an agonist-bound VLS model is used (Lyu et al., 2019;Roth, 2019).
There are several factors, related to both the VLS procedure and the intrinsic properties of the MT receptors, that likely contributed to the high hit rate and agonistic potency of the hits in our study. Thus, the high quality of the initial crystal structure, further improved by ligand-guided optimization of the pocket for VLS, has been critical for the success of our previous VLS campaigns, and likely played a similar role here (Katritch et al., 2010;Lane et al., 2013;Zheng et al., 2017). At the same time, some intrinsic properties of MT receptors also likely facilitated successful VLS for agonists. As we mentioned above, endogenous ligand melatonin itself has unusually high picomolar potency at MT receptors (~4 pM at MT 1 and 50 pM at MT 2 , see Table 1). Melatonin and most other high-potency ligands are small (<250 DA) and yet they still fully occupy the very small, enclosed MT pocket. Chemical space of such size-limited fragment-like libraries is much smaller than the usual drug-like space, and can be more exhaustively searched, likely resulting in higher hit rates. Moreover, most known MT receptor ligands show agonist activity, while antagonists of similar potency are notoriously hard to find, suggesting that agonists may be intrinsically preferred ligands for MT receptors (Jockers et al., 2016).
Two of the hit compounds, 37 and 62, are MT 2 -selective partial agonists, which may have a desirable profile for eliciting antinociceptive effects mediated by melatonin receptors, and may be potentially useful in developing novel analgesics for pain management with reduced side effects (Ló pez-Canul et al., 2015). Of note, while some of the newly discovered hits are selective for MT 2, none of the hits in this study had substantial MT 1 selectivity. This may be explained by the lack of a bulky chemical group at the R1 position, which are known to confer strong MT 1 selectivity, e.g. in bitopic CTL 01-05-B-A05 ligand that stretches out of the pocket via narrow side channel . Design of such bitopic ligands with MT 1 selective chemotypes would need to explore larger compounds (MW >500), which were not considered in the current VLS screen.
This study represents a successful application of structure-based VLS to identify agonists with novel chemotypes, sub-nanomolar potencies, and a high degree of receptor subtype selectivity for a class A GPCR (Lyu et al., 2019;Wang et al., 2017). This study also represents a successful implementation of molecular modeling and structure-based virtual screening, aimed at the melatonin receptors, facilitated by the availability of high-quality structures capturing vital ligand-receptor interactions (Alkozi et al., 2018;Johansson et al., 2019;Stauch et al., 2019). Prevalence of agonists in the hit set suggests the importance of activated, agonist-bound conformations of the orthosteric pocket models for successful agonist screening. Note that, even though the receptors were thermostabilized by 9 and 8 point mutations (MT 1 /MT 2 , respectively) to aid crystallization rendering the receptor conformations inactive on the intracellular side, the agonist-bound orthosteric pockets remain relevant for structure-based drug discovery applications. Our benchmarking also corroborated the important role of ligand guided receptor optimization (LiBERO) (Katritch et al., 2012) in improving the outcomes of a structure-based VLS, similar to some of our previous VLS campaigns (Lane et al., 2013;Zheng et al., 2017). Another critical aspect of this successful VLS is the discovery of novel chemotypes with reliable docking poses. With our screening library assembled to be fragment-like with regards to molecular weights, our hits are diverse and amenable to chemical optimization to improve their pharmacological profiles. Thus, our results also illustrate the utility of fragment-like compounds in the early stages of drug discovery.
The chemical diversity, selectivity, high potency and agonist activities of the identified hits serve as a valuable starting point for the development of tool compounds to explore the biology of melatonin receptors. With the potential for selective modulation over the melatonin receptor subtypemediated biology, these novel chemotypes could provide new leads for the development of nextgeneration treatments for insomnia, pain, sleep and mood-related disorders, type 2 diabetes, and cancer. Receptor model preparation and optimization X-ray crystal structures of MT 1  and MT 2  receptors in complex with 2-phenylmelatonin (PDB IDs Berman et al., 2000: 6ME3, 6ME6) were used to generate virtual screening models. Both structures were converted from PDB coordinates to ICM objects using the object conversion protocol implemented in ICM-Pro (Abagyan et al., 2016). This process includes the addition of hydrogens and assignments of secondary structures, the energetically favorable protonation states to His, Asn and Gln side chains, and of formal charges to the ligand in a complex with the receptor, followed by local minimization of polar hydrogens using energy minimization protocols in ICM-Pro. The orthosteric ligand-binding pocket was further optimized with energy-based global optimization in ICM using Biased Probability Monte-Carlo (BPMC), where the orthosteric ligand and amino acid side chains within 5 Å radius were kept flexible and co-optimized , as described in LiBERO protocol (Katritch et al., 2012) and its previous applications (Lane et al., 2013;Zheng et al., 2017). To validate the models, a set of 20 known MT receptor ligands were selected from ChEMBL database (Gaulton et al., 2017) along with 780 MT receptor decoys selected for each MT 1 and MT 2 receptor from GPCR decoy database (GDD) (Gatica and Cavasotto, 2012) and docked into crystal structures, and optimized ligand models of MT receptors. Following the previously described ligand guided receptor binding pocket optimization protocol, the Receiver Operator Characteristic curves (ROC) were plotted based on the True Positive Rates (TPR) and False Positive Rates (FPR) (Katritch et al., 2012) to evaluate the model optimization. The AUC values were calculated as the areas under these ROC curves and used as a model selection criteria for prospective VLS runs. The RMSD values of ligand binding pocket side chain heavy atoms for MT 1 and MT 2 were 0.51 Å and 0.76 Å , respectively, compared to their corresponding crystal structures.

Materials and methods
To perform additional evaluation of screening results with the thermostabilizing mutants in the proximity of the orthosteric site, as displayed in Figure 1-figure supplement 2, were restored to wild-type (WT) residues. The Phe residue at F251/264 6.48 located 4.3 Å from the ligand was mutated to Trp, followed by local minimization of side-chain conformations using energy-based sampling and minimization protocols . Similarly, A104 3.29 , located 5.2 Å was also restored to Gly in the MT 1 receptor model. Docking to this model suggests that these stabilizing mutations do not substantially impact the binding of known ligands and selected hit candidates into the orthosteric pocket.

Screening library
We selected a subset of commercially available (in-stock and on-demand) fragment-like compounds from the ZINC database with physicochemical properties similar to already reported melatonin receptor ligands (Gaulton et al., 2017;Sterling and Irwin, 2015). We considered compounds with molecular weight 250 Da and logP values ranging 1 to 5 based on the logP data of high-affinity MT ligands (Figure 2-figure supplement 1). The initial dataset comprised of~10 million compounds was converted from SMILES to 3D format, and formal charges were assigned. This set was further reduced to 8.4 million compounds after applying additional filters for net charges (between À1 to 1) and removing compounds with highly reactive functional groups and promiscuous PAINS chemotypes ('molPAINS' and 'bad groups' in ICM-Pro v.3.8-6) (Baell and Holloway, 2010).

Virtual ligand screening
The VLS of 8.4 million compounds library for MT 1 and MT 2 models were performed using the VLS protocol in ICM-Pro . The receptor energy potential maps were calculated using a fine potential grid (0.5 Å ). Several energy terms, including van der Waals, hydrophobic, electrostatic and hydrogen bonding interactions were considered for map calculations. Full torsional flexibility of ligands was allowed, and their internal conformational strain was considered while the receptor atoms were assigned rigid for docking. The docking was performed using BPMC conformational sampling and energy minimization protocol implemented in ICM-Pro for scoring and finding the best docking solutions at the default effort level 1. These top compounds were further docked into corresponding MT receptor models with an increased sampling effort value of 3. Each VLS run for the 8.4 million compound library used 32,000 CPU core hours on 3 Linux workstations with a total of 192 CPU cores. The chemical similarity of selected compounds was calculated using Tanimoto chemical distance function 'Distance(chem1 chem2)', available in Molsoft's ICM-Pro (Totrov, 2008). The fingerprints in this function use a combination of ECFP and linear fingerprints as described in ICM-Pro manual (http://www.molsoft.com/icm/fingerprints.html).

Binding and functional assays Radioligand binding assays
All compounds for in vitro testing were purchased from Enamine, Molport, and Chembridge in stock libraries, with verified identity and guaranteed purity of >95% (37 compounds) or >90% (25 compounds), see Supplementary file 2 for compound QC data).
The Radioligand binding assays were conducted by the NIMH Psychoactive Drug Screening Program (PDSP). The NIMH PDSP is directed by Bryan L. Roth, MD, PhD, at the University of North Carolina at Chapel Hill, North Carolina, and Program Officer Jamie Driscoll at NIMH, Bethesda, MD. Binding assay procedures and protocols are also available online at http://pdspdb.unc.edu/ pdspWeb/?site=assays. All the radioligand binding assays were performed using membranes prepared from transiently transfected HEK293T cells (purchased from ATCC, CRL-11268, authenticated by the supplier using morphology, growth characteristics, and STR profiling and certified mycoplasma-free) and in standard binding buffer (50 mM Tris, 10 mM MgCl 2 , 0.1 mM EDTA, 0.1% BSA, 0.01% ascorbic acid, pH 7.4) as recently reported . [ 3 H]melatonin (PerkinElmer, specific activity = 77.4-84.7 Ci/mmol) is used as the radioligand. Competitive binding assays were performed with various concentrations of test compounds (100 fM to 10 mM), [ 3 H]melatonin (0.2-1.7 nM), and MT 1 or MT 2 membranes in a total volume of 150 mL. Assay plates were sealed and incubated for 4 hr at 37˚C in a humidified incubator until harvesting. Plates were harvested using vacuum filtration onto 0.3% polyethyleneimine pre-soaked 96-well Filtermat A (PerkinElmer) and washed three times with cold wash buffer (50 mM Tris, pH 7.4). Filters were dried and melted with a scintillation cocktail (Meltilex, PerkinElmer). Radioactivity was counted using a Wallac TriLux Microbeta counter (PerkinElmer). Results were analyzed using GraphPad Prism 7.0.

Signaling assays G s and G i/o -cAMP assays
GloSensor cAMP assays were conducted according to the recently published procedure  with minor modifications. Briefly, HEK293 T cells (as above) were transiently cotransfected with receptor (MT 1 or MT 2 ) and GloSensor cAMP (Promega) plasmids overnight, plated in Poly-L-Lysine coated 384-well white clear bottom plates in DMEM + 1% dialyzed FBS. Cells were used for assays at a minimum of 6 hr after plating. Culture medium was first removed and cells were stimulated with drugs in assay buffer (1x HBSS, 20 mM HEPES, 1 mg/ml BSA, 0.1 mg/ml ascorbic acid, pH 7.4) for 15 min at room temperature (this and all the following steps), followed by addition of a mixture of isoproterenol (final of 100 nM) and luciferin (final of 1 mM) for G i/o -cAMP production inhibition assays and luciferin (final of 1 mM) for G s -cAMP production assays. The plates were counted for luminescence after 25 min in a luminescence plate reader. Results were analyzed using GraphPad Prism 7.0.

Tango assays
Tango arrestin recruitment assays were carried out according to the previously published procedure Kroeze et al. (2015). In brief, HTLA cells were transiently transfected with receptor TANGO DNA constructs overnight in DMEM with 10% FBS. Transfected cells were then plated into poly-L-Lys coated 384-well plates using DMEM supplemented with 1% dialyzed FBS. After 6 hr incubation, drug dilutions, prepared in DMEM with 1% dFBS, were added for incubation overnight (16-20 hr). Medium and drug solutions were removed, Bright-Glo reagent (20 uL/well) was added for luminescence counting 20 min later. Results were analyzed in GraphPad Prism 7.0.
Bias factors were estimated according to the published procedure Kenakin et al. (2012) with modifications. Briefly, normalized and pooled results were analyzed by fitting the Black and Leff operational model in Prism 7.0 to obtain Log(t/K A ) values for each pathway (Tango and G i -cAMP). Within each signaling pathway, a difference of Log(t/K A ), DLog(t/K A ), between a test compound and selected reference (melatonin in this case) was calculated. For a testing compound, the difference of DLog(t/K A ), DDLog(t/K A ), was then obtained between two pathways. The bias factor is 10 DDLog(t/KA) .