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-3211chemical 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