Molecular Modeling of the 3D Structure of 5‑HT1AR_Discovery of Novel 5‑HT1AR Agonists via Dynamic Pharmacophore-Based Virtual Screening

JOURNAL OF CHEMICAL INFORMATION AND MODELING pubs. acs. org/jcim Molecular Modeling of the 3D Structure of 5-HT1AR: Discovery of Novel 5-HT1AR Agonists via Dynamic Pharmacophore-Based virtual Screening Lili xut Shanglin Zhon Kunqian Yu. t s Bo Gao, A\ Hualiang jiang *f Xuechu Zhen, and Wei Fu*,t 'Department of Medicinal Chemistry Key Laboratory of Smart Drug Delivery, Ministry of Education, School of Pharmacy,Fudan University, Shanghai 201203, China State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, Department of Pharmacology Soochow University College of Pharmaceuticals Sciences, Suzhou 215123, China S Supporting Information ABSTRACT: The serotonin receptor subtype 1A(5-HTIAR) has een implicated in several neurological conditions, and potent 5- Best Hit: Fwo1 HT1AR agonists have therapeutic potential for the treatment of Ki=51.4 nm EC50=7 nm de Ixiety, schizophrenia, and Parkinsons disease. In the present study, a homology model of 5-HTIaR was built based on the latest released high-resolution crystal structure of the B2AR in its active state(PDB: 3SN6). a dynamic pharmacophore model, which takes the receptor flexibility into account, was constructed, validated, and applied to our dynamic pharmacophore-based virtual y screening approach with the aim to identify potential S-HTIAR agonists. The obtained hits were subjected to S-HTIAR binding and functional assays, and 10 compounds with medium or high K and ECso values were identified. Among them, FWOl(K:=51.9 nM, ECso=7 nM)was evaluated as the strongest agonist for 5 HTLAR The active S-HTIAR model and dynamic pharmacophore model obtained from this study can be used for future discovery and design of novel S-HTIAR agonists. Also, by integrating all computational and available experimental data,a stepwise 5-HTIAR signal transduction model induced by agonist Fwol was proposed. ■ INTRODUCTION antipsychotic), F-13640(highly selective and effective full Serotonin(5-hydroxytrptamine, S-HT), one of the most agonist; undergoing clinical trials for the treatment of pain) important neurotransmitters in the brain, was involved (Table 1). Given the rising interest on 5-HT1AR agonists,new most of the human behavioral and neuropsychological processes, modulating mood, cognition and memory, feeding, Table 1. Pharmacophore Model Validation Using Statistical sexuality, sleep, and pain, etc. Abnormal S-HT transmission isParameters associated with pathogenesis of psychiatric disorders and neurodegenerative disorders. To date, at least 16 serotonin value receptor subtypes have been cloned which are grouped into 1000 seven subfamilies(5-HT1-7)based on their different signaling total number of actives in database(A) echanisms. Among them 5-HTIaR was the first one to be total hits(Ht) fully sequenced and widely studied. Its involvement in anxiety and depression has long been documented -GI evidence has revealed new therapeutic applications of 5- 9 ratio of actives in the hit list HTIAR in the treatment of schizophrenia, Parkinsons disease, enrichment factor or enhancement(E) Accordingly, intensive efforts have been made to explore the fals herapeutic application of 5-HTIAR agonists in the past three H score(goodness of hit list) decades. Representative 5-HTIAR agonists include R-8- E=(Ha/Ht)/(A/D);GH =[((Ha/4HtA)(3A +Ht))(1((HtHa)/ OH-DPAT(full agonist; widely used agonistic tool drug), (DAD))I, GH score of 0.6-0 7 indicates a good model. Buspirone(partial agonist; anxiolytic), Vilazodone(SSRI/S- HTIAR partial agonist; antidepressant),Aripirazole(mixed Received: August 14, 2013 eceptor profile with 5-HTIAR agonistic activity; atypical Published: November 18, ACS Publications o2013 American Chemical Society 3202 dxdoLor/10. 1021/c 00481plJ Chem Inf Model. 2013, 53, 3202-3211

Molecular Modeling of the 3D Structure of 5‑HT1AR: Discovery of Novel 5‑HT1AR Agonists via Dynamic Pharmacophore-Based Virtual Screening Lili Xu,†,§ Shanglin Zhou,‡,§ Kunqian Yu,‡,§ Bo Gao,∥ Hualiang Jiang,*,‡ Xuechu Zhen,*,∥ and Wei Fu*,† † Department of Medicinal Chemistry & Key Laboratory of Smart Drug Delivery, Ministry of Education, School of Pharmacy, Fudan University, Shanghai 201203, China ‡ State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China ∥ Department of Pharmacology, Soochow University College of Pharmaceuticals Sciences, Suzhou 215123, China *S Supporting Information ABSTRACT: The serotonin receptor subtype 1A (5-HT1AR) has been implicated in several neurological conditions, and potent 5- HT1AR agonists have therapeutic potential for the treatment of depression, anxiety, schizophrenia, and Parkinson’s disease. In the present study, a homology model of 5-HT1AR was built based on the latest released high-resolution crystal structure of the β2AR in its active state (PDB: 3SN6). A dynamic pharmacophore model, which takes the receptor flexibility into account, was constructed, validated, and applied to our dynamic pharmacophore-based virtual screening approach with the aim to identify potential 5-HT1AR agonists. The obtained hits were subjected to 5-HT1AR binding and functional assays, and 10 compounds with medium or high Ki and EC50 values were identified. Among them, FW01 (Ki = 51.9 nM, EC50 = 7 nM) was evaluated as the strongest agonist for 5- HT1AR. The active 5-HT1AR model and dynamic pharmacophore model obtained from this study can be used for future discovery and design of novel 5-HT1AR agonists. Also, by integrating all computational and available experimental data, a stepwise 5-HT1AR signal transduction model induced by agonist FW01 was proposed. ■ INTRODUCTION Serotonin (5-hydroxytrptamine, 5-HT), one of the most important neurotransmitters in the brain, was involved in most of the human behavioral and neuropsychological processes, modulating mood, cognition and memory, feeding, sexuality, sleep, and pain, etc.1 Abnormal 5-HT transmission is associated with pathogenesis of psychiatric disorders and neurodegenerative disorders. To date, at least 16 serotonin receptor subtypes have been cloned which are grouped into seven subfamilies (5-HT1−7) based on their different signaling mechanisms.2 Among them 5-HT1AR was the first one to be fully sequenced and widely studied. Its involvement in anxiety and depression has long been documented.3−5 Growing evidence has revealed new therapeutic applications of 5- HT1AR in the treatment of schizophrenia,6 Parkinson’s disease,7 and neural damage.8 Accordingly, intensive efforts have been made to explore the therapeutic application of 5-HT1AR agonists in the past three decades.9−11 Representative 5-HT1AR agonists include R-8- OH-DPAT (full agonist; widely used agonistic tool drug),12 Buspirone (partial agonist; anxiolytic),13 Vilazodone (SSRI/5- HT1AR partial agonist; antidepressant),14 Aripirazole (mixed receptor profile with 5-HT1AR agonistic activity; atypical antipsychotic),15 F-13640 (highly selective and effective full agonist; undergoing clinical trials for the treatment of pain)16 (Table 1). Given the rising interest on 5-HT1AR agonists, new Received: August 14, 2013 Published: November 18, 2013 Table 1. Pharmacophore Model Validation Using Statistical Parameters parameter value total compounds in database (D) 1000 total number of actives in database (A) 25 total hits (Ht) 33 active hits (Ha) 21 % yield of actives 63.6 % ratio of actives in the hit list 84 enrichment factor or enhancement (E) a 25.45 false negatives 4 false positives 12 GH score (goodness of hit list)a 0.68 a E = (Ha/Ht)/(A/D); GH = [((Ha/4HtA)(3A +Ht))(1((HtHa)/ (DA)))], GH score of 0.6−0.7 indicates a good model. Article pubs.acs.org/jcim © 2013 American Chemical Society 3202 dx.doi.org/10.1021/ci400481p | J. Chem. Inf. Model. 2013, 53, 3202−3211

Journal of chemical Information and modelin Article chemical entities with high efficacy and good pharmacokinet favorable interactions with the key residues demonstrated by profiles are still limited. Most reported molecules were available mutagenesis studies. The most reasonable complex originated from several widely studied structural classes such was then submitted to QM/MM minimization encoded in DS as indolylalkylamines and arylpiperazines. Though ligand-based to eliminate bad contacts. The ligand, together with the side approach can lead to molecules with high and even super high chain atoms of D3.32 and S5. 42, was included in the QM s-HTIAR affinity more easily, the diversification of 5-HT1ar region for a quantum calculation using Dmol, and the rest was agonist scaffolds is impaired. Also, the lack of accurate 5-HT1ar handled by CHARMm force field in the MM region. structure is another significant impediment. Molecular Dynamics Simulation. The MD simulations Thanks to the great advances in X-ray crystallography, the were performed using the GROMACS 4.5.1 package. The R- advent of high-resolution crystal structure of the B2AR in its 8-OH-DPAT-S-HTIAR complex was embedded in an explicit active state (PDB: 3SN6) provides us with a good hydrated POPC membrane bilayer. Protein was inserted opportunity to model the 3D structure of 5-HTIAR accurately. according to the Inflate Gro methodology described by Based on the constructed model of 5-HTIAR, we further Kandt, reaching an area per lipid of N75 A. The syster adopted dynamic pharmacophore-based virtual screening was then solvated with SPC waters in a 80 x 80 x 86 A box, approach for the discovery of novel S-HTIAR agonists In this and bad waters were removed. A neutralized system with an approach, a dynamic pharmacophore model was developed and ionic concentration of 154 mmol/L was reached by randomly based virtual screening(DPB-VS)method, first developed by Cl. The resulting system for R-8-OH-DPAT-5-HTIAR Heather A Carlson, has been successfully applied to systems contains 35156 atoms. The Berger lipid parameter was used like HIV-1 integrase, 9 fatty acid amide hydrolase(FAAH),2 for the popc molecules30 in combination with GROMOS96 and Histone Deacetylase 8(HDAC8). Practices have proved 53A6 force field for the protein. The molecular topology ofr that dynamic pharmacophore models perform better than static 8-OH-DPAT was prepared with PRODRG, and the partial ones for they take the flexibility of active site into account. In harge was calculated by using the ChelpG method our study a 100 ns molecular dynamic simulation of the R-8- implemented in the Gaussian 092 with the DFT/B3LYP/6- DH-DPAT-S-HT1AR complex was conducted in order to 311g** basis set. The other two systems of the FwO1-5-HTIAR generate a collection of representative agonistic conformations, complex and the ligand-free S-HTIAR in its in and the active site of 5-HTIaR was then mapped by using five set up in a similar way types of chemical probes. Followed by cluster analysis of Prior to MD simulation, energy minimizations were features, a dynamic pharmacophore model for 5-HT1AR was performed to eliminate poor contacts. After 1000 steps of subjected to docking based-virtual screening(DB-VS) Finally, minimization, the maximum force was converged to less than a set of compounds displaying agonistic activity at 5-HT1AR 10.00 kcal /mol/A. After each system was heated gradually from were revealed by biological tests. Among them, Fwol displays 0 to 310 K by v-scale thermostat, a 1-ns NPT equilibration was the bindig t aode of f wor an d wh r was in vestgerte bey peorer ede wit pr te keen ie and restrained sing the Nd the means of molecular docking and dynamics simulations study. Parrinello-Rhaman method to maintain a constant pressure of 1 Finally, a stepwise 5-HTIAR signal transduction model hereof bar. 500-ps unrestrained equilibration ran afterward. Periodic boundary conditions were applied. a time step of 2.0 fs was employed. All bonds were constrained by the LINCS algorithm. MATERIALS AND METHODS Electronic interactions were calculated using the Partical-Mesh Homology Modeling. The e of 5- Ewald(PME) algorithm. A 100-ns production run was carried TIAR was downloaded from the UniProtKB database(Entry out for the R-8-OH-DPAT-S-HT1AR complex system and 100 code:PO8908),and sequence similarity search was performed ns for the FWO1-5-HTIAR complex and the ligand-free 5. using the NCBI BLAST server. The lately disclosed active- HTIAR systems with coordinates saved every 2 ps for later state structure(PDB code: 3SN6) 7 of B2AR was selected as the template to construct the agonistic conformation of 5- Cluster Analysis. The 5000 protein conformations HTIAR Sequence alignment of s-HTIAR and B,AR was carried extracted every 20 ps from the trajectory of the R-8-OH out using the ClustalX 2.0 program. Homology modeling was DPAT-S-HTIAR complex simulation system were clustered performed with Discovery Studio 3. 5(hereafter abbreviated to based on root-mean-square deviation (RMSD) of the DS). 4 Fifty models were generated after loop refinement, and conformations using the GROMOS conformational cluster the one with the lowest Discrete Optimized Protein Energy analysis method as implemented in the GROMACS. a cutoff (DOPE) score was submitted to energy minimization(1000 value of 1.2 A was employed as the criteria to assign a cluster, teps steepest descent with backbone constrained ). The generating a total of 36 clusters. The representative structures PROCHECK program> was used to evaluate the stereo- from the top 5 popular clusters(-I,- ll, Ill,- IV, -v)were used chemical quality of S-HTIAR to build the dynamical pharmacophore model. Molecular Docking. GoldSuite 5.0 was employed to Active Site Mapping and Pharmacophore Model onduct flexible docking. Briefly, the binding pocket was Generation. Th was used to map the defined to include all residues within 10.0 A of Cy carbon atom active sites of the five representative structures of S-HtIAR to in conserved D3. 32(superscripts refer to Ballesteros-Weinstein detect energetically favorable interactions with the following numbering). Full flexibility was allowed for ligands.The probes: negative ionizable(Co0-), positive ionizable(N1+), number of genetic algorithm(GA)runs was set hydrogen-bond acceptor(O), hydrogen-bond donor(N1),and GoldScore was selected as the scoring function. The top- hydrophobic probes(DRY). The output from the GRID ranking solutions were visually inspected by considerin calculations was visualized and superimposed using VMD 32 dx dolor/10.1021/c400481plJ Chem Inf Model. 2013, 53, 3202-3211

chemical entities with high efficacy and good pharmacokinetic profiles are still limited. Most reported molecules were originated from several widely studied structural classes such as indolylalkylamines and arylpiperazines. Though ligand-based approach can lead to molecules with high and even super high 5-HT1AR affinity more easily, the diversification of 5-HT1AR agonist scaffolds is impaired. Also, the lack of accurate 5-HT1AR structure is another significant impediment. Thanks to the great advances in X-ray crystallography, the advent of high-resolution crystal structure of the β2AR in its active state (PDB: 3SN6)17 provides us with a good opportunity to model the 3D structure of 5-HT1AR accurately. Based on the constructed model of 5-HT1AR, we further adopted dynamic pharmacophore-based virtual screening approach for the discovery of novel 5-HT1AR agonists. In this approach, a dynamic pharmacophore model was developed and applied to search compounds. This dynamic pharmacophore￾based virtual screening (DPB-VS) method, first developed by Heather A. Carlson,18 has been successfully applied to systems like HIV-1 integrase,19 fatty acid amide hydrolase (FAAH),20 and Histone Deacetylase 8 (HDAC8).21 Practices have proved that dynamic pharmacophore models perform better than static ones for they take the flexibility of active site into account. In our study, a 100 ns molecular dynamic simulation of the R-8- OH-DPAT-5-HT1AR complex was conducted in order to generate a collection of representative agonistic conformations, and the active site of 5-HT1AR was then mapped by using five types of chemical probes. Followed by cluster analysis of features, a dynamic pharmacophore model for 5-HT1AR was built. Compounds retrieved from DPB-VS were further subjected to docking based-virtual screening (DB-VS). Finally, a set of compounds displaying agonistic activity at 5-HT1AR were revealed by biological tests. Among them, FW01 displays the strongest agonistic activity toward 5-HT1AR. Furthermore, the binding mode of FW01 and 5-HT1AR was investigated by means of molecular docking and dynamics simulations study. Finally, a stepwise 5-HT1AR signal transduction model hereof was proposed. ■ MATERIALS AND METHODS Homology Modeling. The amino acid sequence of 5- HT1AR was downloaded from the UniProtKB database (Entry code: P08908), and sequence similarity search was performed using the NCBI BLAST server.22 The lately disclosed active￾state structure (PDB code: 3SN6)17 of β2AR was selected as the template to construct the agonistic conformation of 5- HT1AR. Sequence alignment of 5-HT1AR and β2AR was carried out using the ClustalX 2.0 program.23 Homology modeling was performed with Discovery Studio 3.5 (hereafter abbreviated to DS).24 Fifty models were generated after loop refinement, and the one with the lowest Discrete Optimized Protein Energy (DOPE) score was submitted to energy minimization (1000 steps steepest descent with backbone constrained). The PROCHECK program25 was used to evaluate the stereo￾chemical quality of 5-HT1AR. Molecular Docking. GoldSuite 5.026 was employed to conduct flexible docking. Briefly, the binding pocket was defined to include all residues within 10.0 Å of Cγ carbon atom in conserved D3.32 (superscripts refer to Ballesteros-Weinstein numbering27). Full flexibility was allowed for ligands. The number of genetic algorithm (GA) runs was set to 10, and GoldScore was selected as the scoring function. The top￾ranking solutions were visually inspected by considering favorable interactions with the key residues demonstrated by available mutagenesis studies. The most reasonable complex was then submitted to QM/MM minimization encoded in DS to eliminate bad contacts. The ligand, together with the side chain atoms of D3.32 and S5.42, was included in the QM region for a quantum calculation using Dmol3, and the rest was handled by CHARMm force field in the MM region. Molecular Dynamics Simulation. The MD simulations were performed using the GROMACS 4.5.1 package.28 The R- 8-OH-DPAT-5-HT1AR complex was embedded in an explicit hydrated POPC membrane bilayer. Protein was inserted according to the InflateGRO methodology described by Kandt,29 reaching an area per lipid of ∼75 Å. The system was then solvated with SPC waters in a 80 × 80 × 86 Å box, and bad waters were removed. A neutralized system with an ionic concentration of 154 mmol/L was reached by randomly replacing water molecules with the proper number of Na+ and Cl−. The resulting system for R-8-OH-DPAT-5-HT1AR contains 35156 atoms. The Berger lipid parameter was used for the POPC molecules30 in combination with GROMOS96 53A6 force field for the protein. The molecular topology of R- 8-OH-DPAT was prepared with PRODRG,31 and the partial charge was calculated by using the ChelpG method implemented in the Gaussian 0932 with the DFT/B3LYP/6- 311g** basis set. The other two systems of the FW01-5-HT1AR complex and the ligand-free 5-HT1AR in its inactive state were set up in a similar way. Prior to MD simulation, energy minimizations were performed to eliminate poor contacts. After 1000 steps of steepest decent and 200 steps of conjugate gradient energy minimization, the maximum force was converged to less than 10.00 kcal/mol/Å. After each system was heated gradually from 0 to 310 K by v-scale thermostat, a 1-ns NPT equilibration was performed with protein and ligand restrained, using the Nose￾Hoover thermostat to keep the temperature at 310 K, and the Parrinello-Rhaman method to maintain a constant pressure of 1 bar. 500-ps unrestrained equilibration ran afterward. Periodic boundary conditions were applied. A time step of 2.0 fs was employed. All bonds were constrained by the LINCS algorithm. Electronic interactions were calculated using the Partical-Mesh Ewald (PME) algorithm. A 100-ns production run was carried out for the R-8-OH-DPAT-5-HT1AR complex system and 100 ns for the FW01-5-HT1AR complex and the ligand-free 5- HT1AR systems with coordinates saved every 2 ps for later analysis. Cluster Analysis. The 5000 protein conformations extracted every 20 ps from the trajectory of the R-8-OH￾DPAT-5-HT1AR complex simulation system were clustered based on root-mean-square deviation (RMSD) of the conformations using the GROMOS conformational cluster analysis method as implemented in the GROMACS. A cutoff value of 1.2 Å was employed as the criteria to assign a cluster, generating a total of 36 clusters. The representative structures from the top 5 popular clusters (-I, -II, -III, -IV, -V) were used to build the dynamical pharmacophore model. Active Site Mapping and PharmacophoreModel Generation. The GRID 22 program33 was used to map the active sites of the five representative structures of 5-HT1AR to detect energetically favorable interactions with the following probes: negative ionizable (COO−), positive ionizable (N1+), hydrogen-bond acceptor (O), hydrogen-bond donor (N1), and hydrophobic probes (DRY). The output from the GRID calculations was visualized and superimposed using VMD.34 Journal of Chemical Information and Modeling Article 3203 dx.doi.org/10.1021/ci400481p | J. Chem. Inf. Model. 2013, 53, 3202−3211

Journal of Chemical Information and Modeling Article (a) Phi(degrees) Figure 1. Homology model of S-HTIAR a, Homology model of S-HTIAR in cartoon representation. b, Ramachadran plot calculated for the S-HTIAR For each of the five probes, the grid points were superimposed the ability at a concentration of 10 uM to inhibit the binding of to identify clusters. Corresponding pharmacophore feature was a tritiated radioligand to the corresponding receptor was tested. generated in situ at the geometric center of each cluster with Compounds that inhibited binding by more than 90% were the tolerance volume determined by the radius of gyration of further assayed at six or more concentrations, ranging above corresponding cluster. For the pharmacophore model devel oped in this study, clusters were selected based on the binding and below ICso. The K, values were calculated using the mode of R-8-OH-DPAT in 5-HTR. Clusters of the N1+ following equation: K= ICso/(1+C/Kd). [H]5-HT was used as the standard radioligand for S-HTIA receptor. Duplicate residue D3. 32 were selected; clusters of both O and Ni probes tubes were incubated at 30C for 50 min with increasing within a hydrogen-bonding distance from the conserved SS.42 concentrations(1 nM to 100 uM)of each compound and with were selected; two hydrophobic clusters of Dry probe within a 0.7 nM CHS-HT in a final volume of 200 uL of binding buffer distance of 5 A from conserved F6.51, F6.52, and w7. 40 were containing 50 mM Tris and 4 mM MgCl2 (pH 7.4) elected. Thus,a four-feature pharmacophore model of 5- Nonspecific binding was assessed by parallel incubations with HTIAR was generated. Moreover, to reduce the false positive 10 uM S-HT. The reaction was started by addition of membranes(15 pg/tube) and stopped by rapid Itration position of side chains of D3.32, F6.51, and Y7.43 were added. through a Whatman GF/B glass fber filter and subsequent Dynamic Pharmacoph ore-Based Virtual Screening In washing with cold buffer [50 mM Tris and 5 mM ethy our virtual screening approach, the Ligand Pharmacophore potential molecules from Maybridge and Specs chemical 24-well cell harvester. Scintillation cocktail was added, and the Mapping protocol embedded in DS was employed to retrieve enediaminetetraacetic acid(EDTA)(pH 7.4) using a Brandel databases with the dynamic pharmacophore model as a 3D radioactivity was determined in a Micro Beta liquid scintillation query. The two databases, each containing 61623 and 166778 counter. The ICso and K values were calculated by nonlinear compounds,were filtered by Lipinski's"Rule of Five"36 to regression(PRISM, Graphpad, San Diego, CA)using a create druglike databases in advance using Prepare Ligands sigmoidal function. protocol. For each molecule in the database, a maximum of 1 SH-GTPyS Assays. The [S]GTPyS binding assay was conformers with an energy threshold of 20 kcal/mol were performed at 30C for 30 min with 10 ug of membrane protein generated using FASt algorithm. The Best and Flexible in a final volume of 100 uL with various concentrations of the mapping option was adopted. Compounds were required to compounds. The binding buffer contained 50 mM Tris (pH match at least three features of input pharmacophore(Fit Value 7.5),5 mM MgCl, 1 mM ethylenediaminetetraacetic acid 22.95). Hits were then subjected to GOLD docking. Specially, (EDTA), 100 mM NaCl, 1 mM DL-dithiothreitol (DTT), and e automatic genetic algor h the search efficiency was set at 30% which is recommended for 40 uM guanosine triphosphate(GDP). The reaction was virtual eening. Other parameters were the same as described initiated by the addition of [S]GTPys(final concentration of before in the molecular docking section. Outcome solutions 0.1 nM). Nonspecific binding was measured in the presence of were sorted by GoldScore. Compounds with GoldScore >45 uM S-guanylimidodiphosphate(GpP(NH)P). The re- were retained for further visual inspection action was terminated by the addition of 1 mL of ice-cold Binding Assays. All the selected compounds wer washing buffer(50 mM Tris, 5 mM Mg Cl, I mM EDTA, 100 subjected to competitive binding assays for 5-HTIA receptor, mM Naci)and was rapidly filtered with GE/C glass fber filters sing membrane preparation obtained from stable transfected (Whatman)and washed three times. Radioactivity was CHO cells as previously described by our laboratory. First, determined by liquid scintillation counting dxdoLor/10. 1021/c400481p. Chem. Inf Model. 2013, 53, 3202-3211

For each of the five probes, the grid points were superimposed to identify clusters. Corresponding pharmacophore feature was generated in situ at the geometric center of each cluster with the tolerance volume determined by the radius of gyration of corresponding cluster. For the pharmacophore model devel￾oped in this study, clusters were selected based on the binding mode of R-8-OH-DPAT in 5-HT1AR. Clusters of the N1+ probe that were within a distance of 3 Å from the conservative residue D3.32 were selected; clusters of both O and N1 probes within a hydrogen-bonding distance from the conserved S5.42 were selected; two hydrophobic clusters of DRY probe within a distance of 5 Å from conserved F6.51, F6.52, and W7.40 were selected. Thus, a four-feature pharmacophore model of 5- HT1AR was generated. Moreover, to reduce the false positive rate of virtual screening, three excluded volumes right at the position of side chains of D3.32, F6.51, and Y7.43 were added. Dynamic Pharmacophore-Based Virtual Screening. In our virtual screening approach, the Ligand Pharmacophore Mapping protocol embedded in DS was employed to retrieve potential molecules from Maybridge and Specs chemical databases35 with the dynamic pharmacophore model as a 3D query. The two databases, each containing 61623 and 166778 compounds, were filtered by Lipinski’s “Rule of Five” 36 to create druglike databases in advance using Prepare Ligands protocol. For each molecule in the database, a maximum of 100 conformers with an energy threshold of 20 kcal/mol were generated using FAST algorithm. The Best and Flexible mapping option was adopted. Compounds were required to match at least three features of input pharmacophore (Fit Value ≥2.95). Hits were then subjected to GOLD docking. Specially, the automatic genetic algorithm search option was used, and the search efficiency was set at 30% which is recommended for virtual screening. Other parameters were the same as described before in the molecular docking section. Outcome solutions were sorted by GoldScore. Compounds with GoldScore >45 were retained for further visual inspection. Binding Assays. All the selected compounds were subjected to competitive binding assays for 5-HT1A receptor, using membrane preparation obtained from stable transfected CHO cells as previously described by our laboratory.37,38 First, the ability at a concentration of 10 μM to inhibit the binding of a tritiated radioligand to the corresponding receptor was tested. Compounds that inhibited binding by more than 90% were further assayed at six or more concentrations, ranging above and below IC50. The Ki values were calculated using the following equation: Ki = IC50/(1+C/Kd). [3 H]5-HT was used as the standard radioligand for 5-HT1A receptor. Duplicate tubes were incubated at 30 °C for 50 min with increasing concentrations (1 nM to 100 μM) of each compound and with 0.7 nM [3 H]5-HT in a final volume of 200 μL of binding buffer containing 50 mM Tris and 4 mM MgCl2 (pH 7.4). Nonspecific binding was assessed by parallel incubations with 10 μM 5-HT. The reaction was started by addition of membranes (15 μg/tube) and stopped by rapid filtration through a Whatman GF/B glass fiber filter and subsequent washing with cold buffer [50 mM Tris and 5 mM ethyl￾enediaminetetraacetic acid (EDTA) (pH 7.4)] using a Brandel 24-well cell harvester. Scintillation cocktail was added, and the radioactivity was determined in a MicroBeta liquid scintillation counter. The IC50 and Ki values were calculated by nonlinear regression (PRISM, Graphpad, San Diego, CA) using a sigmoidal function. [ 35S]-GTPγS Assays. The [35S]GTPγS binding assay was performed at 30 °C for 30 min with 10 μg of membrane protein in a final volume of 100 μL with various concentrations of the compounds. The binding buffer contained 50 mM Tris (pH 7.5), 5 mM MgCl2, 1 mM ethylenediaminetetraacetic acid (EDTA), 100 mM NaCl, 1 mM DL-dithiothreitol (DTT), and 40 μM guanosine triphosphate (GDP). The reaction was initiated by the addition of [35S]GTPγS (final concentration of 0.1 nM). Nonspecific binding was measured in the presence of 100 μM 5′-guanylimidodiphosphate (Gpp(NH)p). The re￾action was terminated by the addition of 1 mL of ice-cold washing buffer (50 mM Tris, 5 mM MgCl2, 1 mM EDTA, 100 mM NaCl) and was rapidly filtered with GF/C glass fiber filters (Whatman) and washed three times. Radioactivity was determined by liquid scintillation counting. Figure 1. Homology model of 5-HT1AR. a, Homology model of 5-HT1AR in cartoon representation. b, Ramachadran plot calculated for the 5-HT1AR model. Journal of Chemical Information and Modeling Article 3204 dx.doi.org/10.1021/ci400481p | J. Chem. Inf. Model. 2013, 53, 3202−3211

Journal of chemical Information and modeling Article Chart 1. Representative S-HTIAR Agonists Buspirone F13640 Y7 X. Y-80A Figure 2. Side view of the R-8-OH-DPAT-S-HTIAR complex embedded in a hydrated POPC lipid bilayer(left). s-HTIAR is shown in cartoon representation with R-8-OH-DPAT in the upper center. Water molecules and lipid molecules are shown in sticks; and choline N(blue),P atom (orange)and Nat(gray), CL(green) ions are represented in spheres. The detail of the binding mode of R-8-OH-DPAT with S-HTIAR is amplified right).S-HTIAR is shown in gray ribbons and R-8-OH-DPAT in green ball and sticks. Hydrogen bonds are depicted in dotted lines. RESULTS AND DISCUSSION structure of P2 aR bound to the partial inverse agonist carazolol 3D Structure of the 5-HT1AR Model. Sequence alignment (PDB Code: 2RH1) F88), indicated that the sequence identity in the transmembrane our model was built taking the active state conformation of the region is 45% between S-HTIAR and B,AR. The most P2AR-Gs complex as template and represents the active recently disclosed crystal structure of the P,AR-Gs complex in conformation binding with agonists. Thus it is suitable f active state(PDB Code: 3SN6)was the optimal template to me investigation of discovery of 5-HTIAR agonists. construct the S-HTIAR model. Therefore, it was expected tha 8-OH-DPAT(Chart 1), a prototypical 5-HTIAR agonist, has een extensively used for pharmacological studies. R-8-OH- our model represents the active state of the receptor. Figure SI DPaT displays higher efficacy than its S-counterpart and acts as shows the final sequence alignment of 5-HT1ar to the a full S-HTIAR agonist, thus it was chosen herein to induce template. The PROCHECK statistics showed that 99.2% of active receptor conformations. The predicted binding mode the residues in the 5-HT,AR model (Figure la) were either in the most favored or in the additionally allowed regions of the ( Figure 2, right)is quite consistent with all theavailable experimental studies: the protonated nitrogen forms a salt Ramachandran plot; only one residue in the loop is located in negatively charged D3. 32, and the hydroxyl the disallowed region( Figure 1b), suggesting that the overall group is involved in the hydrogen bonds with Ser5. 42 and main chain and side chain conformations are reasonable Lys191(ECL2). The importance of conserved D3. 32 and S5.42 Besides, several structural features present in the 5-HT1aR for 5-HT binding to S-HTIAR was supported by a site-directed model are characteristic of active conformations of family a mutagenesis study. Additionally, F6.51 and F6.52are GPCRs: the broken Arg3.50-Glu630 ionic lock due to a large participated in the aromatic stacking with R-8-OH-DPAT separation of TM3-TM6 in the cytoplasmic ends and Tyr7.53 The two phenylalanine residues in TM6 contribute to of the NPxxY motif (so-called"Tyrosine Toggle Switch") hydrophobic interaction and were proved to be crucial for extending into the protein interior. Compared to previously ligand binding in the case of the dopamine D2 receptor. The eported S-HTIAR models which were built based on the crystal di-n-propyl substituent was found stretched in the hydrophobic dxdoLor/10. 1021/c400481p. Chem. Inf Model. 2013, 53, 3202-3211

■ RESULTS AND DISCUSSION 3D Structure of the 5-HT1AR Model. Sequence alignment indicated that the sequence identity in the transmembrane region is ∼45% between 5-HT1AR and β2AR. The most recently disclosed crystal structure of the β2AR-Gs complex in active state (PDB Code: 3SN6) was the optimal template to construct the 5-HT1AR model. Therefore, it was expected that our model represents the active state of the receptor. Figure S1 shows the final sequence alignment of 5-HT1AR to the template. The PROCHECK statistics showed that 99.2% of the residues in the 5-HT1AR model (Figure 1a) were either in the most favored or in the additionally allowed regions of the Ramachandran plot; only one residue in the loop is located in the disallowed region (Figure 1b), suggesting that the overall main chain and side chain conformations are reasonable. Besides, several structural features present in the 5-HT1AR model are characteristic of active conformations of family A GPCRs:39 the broken Arg3.50-Glu6.30 ionic lock due to a large separation of TM3-TM6 in the cytoplasmic ends and Tyr7.53 of the NPxxY motif (so-called “Tyrosine Toggle Switch”) extending into the protein interior. Compared to previously reported 5-HT1AR models which were built based on the crystal structure of β2AR bound to the partial inverse agonist carazolol (PDB Code: 2RH1) or bovine rhodopsin (PDB Code: 1F88), our model was built taking the active state conformation of the β2AR-Gs complex as template and represents the active conformation binding with agonists. Thus it is suitable for the investigation of discovery of 5-HT1AR agonists. 8-OH-DPAT (Chart 1), a prototypical 5-HT1AR agonist, has been extensively used for pharmacological studies.40 R-8-OH￾DPAT displays higher efficacy than its S-counterpart and acts as a full 5-HT1AR agonist, thus it was chosen herein to induce active receptor conformations. The predicted binding mode (Figure 2, right) is quite consistent with all the available experimental studies: the protonated nitrogen forms a salt bridge with the negatively charged D3.32, and the hydroxyl group is involved in the hydrogen bonds with Ser5.42 and Lys191 (ECL2). The importance of conserved D3.32 and S5.42 for 5-HT binding to 5-HT1AR was supported by a site-directed mutagenesis study.41 Additionally, F6.51 and F6.52 are participated in the aromatic stacking with R-8-OH-DPAT. The two phenylalanine residues in TM6 contribute to hydrophobic interaction and were proved to be crucial for ligand binding in the case of the dopamine D2 receptor.42 The di-n-propyl substituent was found stretched in the hydrophobic Chart 1. Representative 5-HT1AR Agonists Figure 2. Side view of the R-8-OH-DPAT-5-HT1AR complex embedded in a hydrated POPC lipid bilayer (left). 5-HT1AR is shown in cartoon representation with R-8-OH-DPAT in the upper center. Water molecules and lipid molecules are shown in sticks; and choline N (blue), P atoms (orange) and Na+ (gray), CL− (green) ions are represented in spheres. The detail of the binding mode of R-8-OH-DPAT with 5-HT1AR is amplified (right). 5-HT1AR is shown in gray ribbons and R-8-OH-DPAT in green ball and sticks. Hydrogen bonds are depicted in dotted lines. Journal of Chemical Information and Modeling Article 3205 dx.doi.org/10.1021/ci400481p | J. Chem. Inf. Model. 2013, 53, 3202−3211

Journal of Chemical Information and Modeling Article (b) 22● Figure 3. a, Superimposition of the five major S-HTIAR conformations obtained from 100 ns MD simulation. Cluster-I representative structure is shown in white surface. Conserved residues in the binding pocket of Cluster -L, -Il,Ill,IV, and -V representative conformations are eparately in white, blue, green, yellow, and pink. b, Top view of the pharmacophore model mapped into the active pocket of s-HTIAR(only I representative structure was shown for clarity). Two hydrophobic elements(Hydro-l, -2), a positive ionizable element(POS), a hydrog acceptor element, and three excluded volumes(Excl.1,2,3)are represented in cyan, red, green, and black spheres, respectively. pocket(enclosed by ECL2 and the extracellular segments of Validation of a Pharmacophore Model. The Guner TM-3,-7)which is also typical among dopamine receptors. MD Henry (GH) scoring method"was used to assess the ability of simulation was started with the obtained R-8-OH-DPAT-5- our pharmacophore model to discriminate a small number of HTIAR complex embedded in a hydrated POPC lipid bilayer known active molecules against a greater number of decoy (Figure 2, left). During the simulation 8-OH-DPAT was molecules in the database. The decoy set contains a total of constantly trapped in the binding pocket of 5-HTIAR Taken 1000 compounds, including 25 highly active 5-HTIAR agonists together, our pi predicted 5-HT1AR model is credible, and taken from ChEMBL database, and 975 matched decoys conformations generated from MD simulation shall stand for generated from online DUD-E(Directory of Useful Decoys active states of the recept Enhanced)tool. The pharmacophore model was used to The Dynamic Pharmacophore Model Generation. Five screen the decoy set employing the BEST flexible searching representative structures obtained from clustering were used to method. A set of statistical parameters were set to analyze the represent the flexibility of s-HTIAR and to build the dynamic result (Table 1). 33 compounds were retrieved as hits with a hit pharmacophore model. GRID Probed the active site with five rate of 63.6%. In addition, the pharmacophore model showed types of functional groups(COO-, NI+, O, Nl, and DRY) to an enrichment factor of 25.45 and a GH score of 0.68, find the most favorable regions for them to interact with the indicating the good quality of our model protein, generating a total of 25 interaction maps(5 for each Database Searching and Bioassay Results. The protein conformation), and they were overlaid based on the validated pharmacophore model was then used as a 3D query superimposition of the five major S-HTIAR conformations to screen Maybridge and Specs databases, of which nondruglike (Figure 3a). Backbone atoms RMSD of representative compounds were rejected. A total of 18, 976 compounds respectively. On the basis of site-directed mutagenesis data and with GoldScore >45 were extracted for visual inspection. omputational studies available in the literature, #+344ie Compounds that form favorable interactions with key residues hydrogen bond contacts with D3. 32, Ser5. 42, T5.43, N7.39, in the active site such as D3. 32, SS. 42, F6.51, and F6.52 Y7.43 and T-t stacking interactions with Y2.64, w6.48, F6.51, chosen, and those with unreasonable binding modes F6.52, W7.40, Y7. 43 are crucial for S-HTIAR agonist binding discarded. At last, a total of 45 compounds, 16 from Specs only the features complementary to these key residues were 29 from Maybridge, were selected to purchase from onsidered. Generally, a positive ionizable feature(POS)was suppliers. The flowchart of hybrid searching approach was chosen to target the negatively charged carboxylic side chain of shown in Figure 4 D3. 32, the hydrogen-bond acceptor(HBA)feature located ge clusters found in Specs, May Bridge hydrophobic regions, which are surrounded by TMS/TM6 and 228401 Lipinski's Ro5& Veber's Rules TM2/TM7/ECL2, were combined into two hydrophobic 210,483 features(HYDRO-1,-2), respectively. Whereas, clusters of COO- probes were discarded due to the electrostatic repulsion between COO-and negatively charged residues in the binding D8-VS site of S-HTIAR(Figure 7), and the hydrogen-bond donor Visual inspection to SS.42/T5.43 was also omitted for it partially Bioassay overlapped with HYDRO-1. Thus, a four-feature pharmaco- phore model was produced( Figure 3b). Three exclusion volumes(Excl-1,-2, -3)were in situ generated from the side Figure 4. Workflow of dynamic pharmacophore-based virtual chains of the conserved residues D3. 32, F6.51, and Y7.43. searching approach. dx dolor/10.1021/c400481plJ Chem Inf Model. 2013, 53, 3202-3211

pocket (enclosed by ECL2 and the extracellular segments of TM-3, -7) which is also typical among dopamine receptors. MD simulation was started with the obtained R-8-OH-DPAT-5- HT1AR complex embedded in a hydrated POPC lipid bilayer (Figure 2, left). During the simulation 8-OH-DPAT was constantly trapped in the binding pocket of 5-HT1AR. Taken together, our predicted 5-HT1AR model is credible, and conformations generated from MD simulation shall stand for active states of the receptor. The Dynamic Pharmacophore Model Generation. Five representative structures obtained from clustering were used to represent the flexibility of 5-HT1AR and to build the dynamic pharmacophore model. GRID probed the active site with five types of functional groups (COO-, N1+, O, N1, and DRY) to find the most favorable regions for them to interact with the protein, generating a total of 25 interaction maps (5 for each protein conformation), and they were overlaid based on the superimposition of the five major 5-HT1AR conformations (Figure 3a). Backbone atoms RMSD of representative structures in each Cluster (II, III, IV, V) to the major conformer (Cluster I) is 1.48 Å, 1.38 Å, 1.74 Å, 1.39 Å, respectively. On the basis of site-directed mutagenesis data and computational studies available in the literature,41,43,44 i.e. hydrogen bond contacts with D3.32, Ser5.42, T5.43, N7.39, Y7.43 and π−π stacking interactions with Y2.64, W6.48, F6.51, F6.52, W7.40, Y7.43 are crucial for 5-HT1AR agonist binding, only the features complementary to these key residues were considered. Generally, a positive ionizable feature (POS) was chosen to target the negatively charged carboxylic side chain of D3.32, the hydrogen-bond acceptor (HBA) feature located close to Y7.43 was retained, and large clusters found in two hydrophobic regions, which are surrounded by TM5/TM6 and TM2/TM7/ECL2, were combined into two hydrophobic features (HYDRO-1, -2), respectively. Whereas, clusters of COO- probes were discarded due to the electrostatic repulsion between COO- and negatively charged residues in the binding site of 5-HT1AR (Figure 7), and the hydrogen-bond donor feature adjacent to S5.42/T5.43 was also omitted for it partially overlapped with HYDRO-1. Thus, a four-feature pharmaco￾phore model was produced (Figure 3b). Three exclusion volumes (Excl-1, -2, -3) were in situ generated from the side chains of the conserved residues D3.32, F6.51, and Y7.43. Validation of a Pharmacophore Model. The Gü ner￾Henry (GH) scoring method45 was used to assess the ability of our pharmacophore model to discriminate a small number of known active molecules against a greater number of decoy molecules in the database. The decoy set contains a total of 1000 compounds, including 25 highly active 5-HT1AR agonists taken from ChEMBL database,46 and 975 matched decoys generated from online DUD-E (Directory of Useful Decoys, Enhanced) tool.47 The pharmacophore model was used to screen the decoy set employing the BEST flexible searching method. A set of statistical parameters were set to analyze the result (Table 1). 33 compounds were retrieved as hits with a hit rate of 63.6%. In addition, the pharmacophore model showed an enrichment factor of 25.45 and a GH score of 0.68, indicating the good quality of our model. Database Searching and Bioassay Results. The validated pharmacophore model was then used as a 3D query to screen Maybridge and Specs databases, of which nondruglike compounds were rejected. A total of 18,976 compounds obtained from PB-VS were docked into the representative structure of 5-HT1AR in the biggest cluster. 1500 compounds with GoldScore >45 were extracted for visual inspection. Compounds that form favorable interactions with key residues in the active site such as D3.32, S5.42, F6.51, and F6.52 were chosen, and those with unreasonable binding modes were discarded. At last, a total of 45 compounds, 16 from Specs and 29 from Maybridge, were selected to purchase from their suppliers. The flowchart of hybrid searching approach was shown in Figure 4. Figure 3. a, Superimposition of the five major 5-HT1AR conformations obtained from 100 ns MD simulation. Cluster-I representative structure is shown in white surface. Conserved residues in the binding pocket of Cluster -I, -II, -III, -IV, and -V representative conformations are colored separately in white, blue, green, yellow, and pink. b, Top view of the pharmacophore model mapped into the active pocket of 5-HT1AR (only Cluster￾I representative structure was shown for clarity). Two hydrophobic elements (Hydro-1, -2), a positive ionizable element (POS), a hydrogen bond acceptor element, and three excluded volumes (Excl-1, -2, -3) are represented in cyan, red, green, and black spheres, respectively. Figure 4. Workflow of dynamic pharmacophore-based virtual searching approach. Journal of Chemical Information and Modeling Article 3206 dx.doi.org/10.1021/ci400481p | J. Chem. Inf. Model. 2013, 53, 3202−3211