Molecular modelling of human 5-hydroxytryptamine receptor (5-HT2A) and virtual screening studies towards the identification of agonist and antagonist molecules

The serotonin receptors, also known as 5-hydroxytryptamine (5-HT) receptors, are a group of G protein-coupled receptors (GPCRs) and ligand-gated ion channels found in the central and peripheral nervous systems. GPCRs have a characteristic feature of activating different signalling pathways upon ligand binding and these ligands display several efficacy levels to differentially activate the receptor. GPCRs are primary drug targets due to their central role in several signal transduction pathways. Drug design for GPCRs is also most challenging due to their inherent promiscuity in ligand recognition, which gives rise to several side effects of existing drugs. Here, we have performed the ligand interaction study using the two prominent states of GPCR, namely the active and inactive state of the 5-HT2A receptor. Active state of 5-HT2A receptor model enhances the understanding of conformational difference which influences the ligand-binding site. A 5-HT2A receptor active state model was constructed by homology modelling using active state β2-adrenergic receptor (β2-AR). In addition, virtual screening and docking studies with both active and inactive state models reveal potential small molecule hits which could be considered as agonist-like and antagonist-like molecules. The results from the all-atom molecular dynamics simulations further confirmed that agonists and antagonists interact in different modes with the receptor.


Introduction
Serotonin (5-HT), a major neurotransmitter found in the central nervous system, is also present in many peripheral tissues. Its numerous biological functions are mediated by different types of serotonin receptors (Uphouse, 1997). These serotonin receptors, also known as 5-hydroxytryptamine receptors (5-HTR) or 5-HT receptors are a group of G protein-coupled receptors (GPCRs) and ligand-gated ion channels. They mediate both excitatory and inhibitory neurotransmission (Cowen, 1991). The serotonin receptors influence various biological and neurological processes such as aggression, anxiety, appetite, cognition, learning, memory, mood, nausea, sleep, and thermoregulation. During the past two decades, multiple 5-HT receptor subtypes have been characterized ranging from 5-HT 1 to 5-HT 7 (Hoyer et al., 1994). Of these, the 5-HT 2 receptor family currently comprises three receptor subtypes, 5-HT 2A , 5-HT 2B , and 5-HT 2C receptors, which are similar in terms of their molecular structure, pharmacology, and signal transduction pathways (Barnes & Sharp, 1999). In all these 5-HT 2 subtypes, 5-HT 2A receptor have shown extensive clinical implications since it is involved in psychological disorders ranging from schizophrenia, depression, obsessive compulsive disorder, attention deficit-hyperactivity disorder, to eating disorders such as anorexia nervosa and autism spectrum disorders (Nichols & Nichols, 2008). This particular receptor has been found to bind to a number of agonists, antagonists, and inverse agonists, which in turn regulate the signalling via G-proteins to treat specific disease conditions (Pytliak, Vargová, Mechírová, & Felšöci, 2011). For example, altered 5-HT receptor level is implicated in schizophrenia (De Angelis, 2002) and depression (Celada, Puig, Amargos-Bosch, Adell, & Artigas, 2004).
An important aspect of this GPCR signalling mechanism is that the ligand does not pass through the membrane. Instead, the signal is transferred to the inside of the cell by conformational changes in the receptor protein, which are coupled to ligand-binding events (De Lean, Stadel, & Lefkowitz, 1980). The receptor protein in the membrane is dynamic and can take on a number of conformations; two prominent ones are called the inactive and active conformations (Perez & Karnik, 2005). Conformational changes at the ligand-binding site and other places in active and inactive forms of β 2 -AR will be helpful for 5-HT 2A receptor drug discovery.
Recently, high-resolution structures were determined in class A GPCRs (Tate & Schertler, 2009). Many of these structures are determined in inactive conformations, bound to either inverse agonists that reduce basal activity or neutral antagonists that maintain basal activity. The four crystal structures of class A GPCRs: bovine rhodopsin (Choe et al., 2011;Deupi et al., 2012;Scheerer et al., 2008;Standfuss et al., 2011), human β 2 -AR (Rasmussen et al., 2011a(Rasmussen et al., , 2011b, human A 2A R (Lebon et al., 2011;Xu et al., 2011) and rat NTSR1 (White et al., 2012) were obtained in active (or intermediateactive) states. An important landmark in GPCR biology was the determination of the active-state ternary complex of β 2 -AR in complex with the heterotrimeric G protein (Rasmussen et al., 2011b). These structures provide unprecedented insights into the structural and functional diversity of this protein family (Venkatakrishnan et al., 2013).
Many computational studies were previously performed to model transmembrane proteins and to understand the ligand-binding modes (AbdulHameed, Hamza, & Zhan, 2006;Hamza, AbdulHameed, & Zhan, 2008). In a previous study, sequence analysis followed by modelling and docking studies of 5-HT 2A receptor was performed to reveal hotspot residues important for ligand binding (Kanagarajadurai, Malini, Bhattacharya, Panicker, & Sowdhamini, 2009). The structure of 5-HT 2A receptor was modelled using structure of the inactive state of β 2 -AR (PDBID:2RH1) (Cherezov et al., 2007) and the residues involved in ligand binding were reported. After the determination of the active state of β 2 -AR structure, several studies have been attempted to understand the difference in ligand binding of the active and inactive state of the receptor (Deupi & Standfuss, 2011;Huber, Menon, & Sakmar, 2008;Lebon et al., 2011). Modelling-and QSAR-based identification of antagonist for 5-HT 2 subtype was performed to achieve the discrimination of ligand interaction between the subtypes (Brea et al., 2002). There was also an attempt to model the active state of 5-HT 2A receptor using steered molecular dynamics (MD) studies (Ísberg, Balle, Sander, Jørgensen, & Gloriam, 2011). Recently, a similar attempt was reported to design biased agonist using active state 5-HT 2A receptor (Marti-Solano et al., 2015). Moreover, a study about the single amino acid difference (Ser239 in [TM5]) between mouse and human and their impact in ligand interaction (Canal et al., 2013) was reported.
In this paper, we report a systematic approach for the derivation of active state model, investigation of important ligand-binding residues by known molecule docking, virtual screening (VS) of small molecules, validating the hits by enrichment analysis, and followed by MD study to propose new agonist-and antagonist-type molecules. We also report possible differences between active and inactive state model in the context of their interactions with the ligand. We have also observed that certain residues are preferred for interaction between active and inactive state of 5-HT 2A receptor by means of docking and MD studies.
Different ligands can differentially activate signalling pathways of GPCRs via the same receptor. This differential activation arises independent of differences in the affinity of the ligand to the receptor. Ligands often display differences in efficacy and/or potency at one signalling pathway vs. another (Raote, Bhattacharya, & Panicker, 2007). In recent years, few studies have attempted to gain the knowledge on influence of functional residues on the ligand which leads to different efficacy levels (Goddard & Abrol, 2007;Kanagarajadurai et al., 2009;Rasmussen et al., 2011a;Strange, 2008). The current structural data about GPCRs provide enormous information about active and inactive conformational changes and residues involved in interaction.
The intention of this analysis was to investigate the difference between the active and inactive forms of 5-HT 2A receptor structure using available structural information. The endogenous ligand and some of the known agonist-and antagonist-type molecules were docked to decipher the binding patterns and important residues involved in agonist-and antagonist-type interactions. This will eventually help to understand binding characteristics of agonist and antagonist interactions with 5-HT 2A receptor. Further, using this knowledge, new agonist-and antagonist-type molecules were identified through VS. Difference in mode of interaction between the active and inactive state of the receptor was enabled by induced fit docking (IFD) experiments. Further, MD was performed to investigate the conformational changes and stability of the proposed model and receptor-ligand interactions. Finally, multiple sequence alignment of serotonin receptors subtypes (except 5-HT 3 ), especially of transmembrane region was pursued, which would provide insights about some important residues specific to the 5-HT 2A receptor subtype in order to avoid promiscuous binding.

Molecular modelling and validation
Until 2013, the structures of serotonin receptors were modelled using β 2 AR active and inactive state of the structure (Cherezov et al., 2007;Rasmussen et al., 2011a). In 2013, an important milestone in the research of serotonin receptor is the determination of crystal structure of 5-HT 1E  and 5-HT 2B . Since both the structures were from serotonin receptor subtype, it should be appropriate to use the structure information for modelling 5-HT 2A receptor. However, the differential signalling patterns are mirrored in the crystal structures, which show features of an intermediate active state for the 5-HT 1B receptor, and a β-arrestin-biased activation state for the 5-HT 2B receptor. It was also reported that the residues involved in biased signalling (β-arrestin mediated) and unbiased signalling (GPCR-mediated) are different. Hence, we used the active state β 2 AR structure as a template for modelling the 5-HT 2A receptor. As such, the active-and inactive state-mediated structural changes are similar in both β 2 -AR and serotonin receptors. In the previous study, the inactive state of the 5-HT 2A receptor was modelled (Kanagarajadurai et al., 2009) using inactive state of β 2 AR (PDBID: 2RH1) structure and the same modelled structure was utilized for this study for the inactive state.
The crystal structure of the human β 2 adrenergic receptor (β 2 -AR) (PDBID: 3SN6) (Rasmussen et al., 2011b) was used as a template for modelling the active state of 5-HT 2A receptor. The human 5-HT 2A receptor sequence (ID: P28223) and the sequence of human (β 2 -AR) (PDBID: 3SN6) were aligned using PRALINE server (Pirovano, Feenstra, & Heringa, 2008). PRALINE is a multiple sequence alignment tool, designed especially for aligning transmembrane proteins (Figure 1) which would help to improve the alignment in the TM regions. The structure of the template was obtained from RCSB (Berman et al., 2000) (http://www.rcsb.org/pdb). The coordinates corresponding to 1-30, 239-264 and 341-365 segments were not available in the template crystal structure due to poor electron density and hence the corresponding residues were removed from the query sequence before the alignment. The final alignment was used to construct the model using the software MODEL-LER (Šali & Blundell, 1993) (Version 9.1). A set of 40 structures were generated, from which the lowest energy structure was used for further processes. The initial lowenergy model, obtained from MODELLER, was validated using PROCHECK (Laskowski, MacArthur, Moss, & Thornton, 1993) server. The structure was further energy minimized using the SYBYL software package (Version 7.1) (Tripos Associates Inc.) and utilized to build the N-terminal (1-70 residues), C-terminal regions (402-471) and ICL3 (243-276 residues) regions using pre-installed PRODAT database tool in SYBYL (Tripos Associates). The structure was again energy minimized using the same parameters successively after the construction of the loop regions. The final structure was validated again using PROCHECK (Laskowski et al., 1993).

Docking studies on endogenous and known ligands
In order to examine the binding characteristics of endogenous ligand serotonin and other known molecules, docking studies were performed using the model structure. This study would enhance the understanding of the influence of functionally important residues on the agonist and antagonist molecules. A total of 17 known ligands (8 agonists and 9 antagonists; please see Results): agonists such as Serotonin, DOI (2,5-Dimethoxy-4-iodoamphetamine), Mescaline, LSD (Lysergide), 5-MeO-alpha-ET, Psilocin, Bufotenine, Dimethyltryptamine) and Antagonists such as Nefazodone, Aripiprazole, Haloperidol, Cyproheptadine, Trazodone, Clozapine, Ketanserin, Spiperone, Risperidone were collected from ZINC database (Irwin & Shoichet, 2005) and their experimental affinity (Ki) values were retrieved from ChEMBL database (Gaulton et al., 2011). Both active state model and previously proposed inactive state model (Kanagarajadurai et al., 2009) were used for docking. The docking studies were performed using Schrödinger Glide software (version 9.3, Schrödinger, Inc., New York, NY, 2009) with IFD module. The conformational search area was defined by a grid, which covers maximum region of all the TM helices. From the top 20 poses, the best pose was selected based on known residue interactions and GScore (Glide Score). The experimental values (Ki) of binding affinity for these molecules against 5-HT 2A are in nM range. They have been converted using free energy calculations (ΔG°= −RTlnK) where the results are in kcal/mol.

VS and IFD
Initially, endogenous ligand and known agonist and antagonist molecules were docked to understand the binding pattern, and the range of GScore was also obtained. This would also provide information about the site and the residues important for binding, which was utilized for selecting the best molecules from the VS experiments. In order to enrich the data-set of agonistand antagonist-type molecules available for 5-HT 2A receptor, VS followed by IFD is next performed to identify more potential small molecules which could function as agonist-or antagonist-type molecules.
VS was performed using small molecule library consisting of 1,40,809 molecules from the clean lead-like category of Zinc database (Irwin & Shoichet 2005). These molecules satisfy Lipinski's rule of five (Lipinski, 2000) and molecular weight is in the range of 150-350, which is appropriate for GPCR ligand-binding site and readily available for docking. The VS workflow has number of tasks prior to docking. Initially, the molecules were considered for Qikprop calculation and Lipinski's rule validation, removal of any reactive groups. Then, ligand preparation task was performed to correct the improper geometry, removal of duplicate entry, and finally generate possible conformations of the ligand. The centroid for VS was set in the ligand-binding site composed of TM3, TM5, TM6, and TM7 based on the known ligand-binding site of GPCRs. The screening was performed by three consecutive methods: namely high throughput virtual screening, special precision (Halgren et al., 2004), and XP (Extra precision) (Friesner et al., 2006) modes with a different sampling and scoring method. The final XP method performs extensive sampling with the stringent scoring calculation for the best binding pose.
The hits from VS were short-listed for further IFD docking on the basis of known interactions of 5-HT 2A with agonist-and antagonist-type molecules. IFD method was employed to find the more accurate complex structure and to rescue the false negatives (relatively poorly scoring true binders) in VS experiments which could be docked in an improved manner. The choice of an IFD method is quite important for GPCR structures, since the ligand-binding site is very small and introducing flexibility to that region would definitely improve interactions between the amino acids and the ligand. This was basically achieved by selection of a centroid point and ligands were docked by Glide docking protocol using a softened potential (Van der Waals radii scaling) and optional removal of side chains and application of constraints. The top poses were further refined using Prime side chain prediction for each protein/ligand complex, with primary focus on residues within a given distance of any ligand pose (default 5 Å), with optional inclusion or exclusion of other residues. The receptor structure in each pose was an induced fit to the ligand structure and conformation. Further, Glide re-docking of ligand structure was performed into the induced-fit receptor (Kalid, Toledo Warshaviak, Shechter, Sherman, & Shacham, 2012) with implicit membrane model for the receptor.

Scaffold-based clustering
The VS hits were further analysed for similarity/dissimilarity with the known molecules. All the hits, along with the known agonists and antagonists chosen for study, were clustered using CANVAS module (Sastry, Lowrie, Dixon, & Sherman, 2010) in Schrodinger software. All the molecules were clustered based on scaffold similarity using binary and hierarchical clustering. The diversity of ligands was computed as follows: the MACCs fingerprints of the active ligands were determined, and hierarchical clustering was performed on the matrix of pairwise similarities determined using the Tanimoto metric (Flower, 1998) using the program CANVAS.  (Park, Scheerer, Hofmann, Choe, & Ernst, 2008), covered by red box, are conserved in all the 5-HT GPCR class, except in 5-HT 6 receptor. The dashed red line indicates the disulphide bridge formed across the CYS residues in TM3 (C148) and ECL2 (C227).

Validation of VS hits using enrichment calculation
The key requirement for success in VS is that the docking methodology should rank actives very early in the larger set of compounds (Truchon & Bayly, 2007). For this purpose, a simple enrichment calculation was performed using known actives, VS hits along with the decoy set docked in both active and inactive state model. This analysis would help to understand the ability of the docking protocol to avoid false positive and false negatives. The known actives (62 compounds) were collected from ZINC database (Irwin & Shoichet, 2005) and decoy set (1000 compounds) were collected from GDD database (Gatica & Cavasotto, 2012). All the compounds (62 actives + 428 VS hits + 1000 decoys) were docked in both active and inactive form and ranked by GScore. The enrichment was calculated using ROC plot which includes information about sensitivity and specificity. Enrichment curves are a useful tool to characterize the ability of a docking program to select active compounds and discard inactive compounds. ROC scores vary between 0 and 1.0; closer the score is to 1.0, the better the classifier at distinguishing actives from decoys (Daga, Polgar, & Zaveri, 2014;Kalva, Vinod, & Saleena, 2014). Two sets of analysis, such as recognition of only actives and recognition of VS hits in both active and inactive form, were performed. Enrichment factors were calculated at 2, 5, and 10% of the database which reflects the ability of docking method to find true positives among the decoys in the database compared to a random selection.

MD simulations
MD simulations were performed to evaluate the stability and conformational changes of the proposed model and also to investigate the protein-ligand contacts, especially the residues involved in interaction with agonist-and antagonist-type molecules at different timescales. For this purpose, MD simulation was performed on the active state 5-HT 2A (Banks et al., 2005). The active state model was pre-equilibrated in the lipid layer composed of POPE (1-palmitoyl-2-oleoyl-Sn-glycero-3-phosphoethanolamine) using the system builder tool of Desmond in Maestro 9.6. Water molecules were added to the ligand-receptor complex with simple point charge water model. MD simulation was performed with the Desmond system builder where the protein and ligand were packed in the orthorhombic box (sides a = 90 Å, b = 90 Å, c = 90 Å), allowing for a 10 Å buffer region between protein atoms and box sides. The system was neutralized with the appropriate counter ions with an ion neutralizer. All the simulations were run in the NPT ensemble, and the temperature was kept constant at 300 K by coupling to the Nose-Hoover thermostat (Hoover, 1985). Once the system was built, the solvated system was subjected to energy minimization using steepest descent algorithm. The minimized protein was further relaxed before the actual simulation by a series of energy minimizations and short MD simulations. The energy and trajectory atomic coordinates for all simulations were recorded at every 250 ps and 1 ns. The simulations were carried out for 10 ns for active state model and best agonist and antagonist complexes. Energy fluctuations and RMSD of the complex in each trajectory were analysed with respect to simulation time. The backbone and side chain root mean square fluctuations (RMSF) of each residue of active state model were monitored for consistency. The intermolecular interactions of agonist and antagonist type interactions were assessed for stability of the docking complex. Such dynamics simulations were done merely to ascertain the robust interactions between ligand and receptor; hence, a relatively short time scale of 10 ns was considered.

Multiple sequence alignment of serotonin receptors
The multiple sequence alignment program PRALINE ™ (Pirovano et al., 2008) was used for the alignment of all the subtypes of human serotonin receptors, except 5-HT 3 . PRALINE uses progressive alignment algorithms to handle TM topology and alignment of TM regions using a state-of-the-art predictor.

Results and discussion
Modelling of active state 5-HT 2A receptor and validation We modelled the active state of the 5-HT 2A receptor using the active state β 2 -AR structure as a template.
The human 5-HT 2A receptor shares 34% sequence identity with human β 2 -AR. The sequence alignment was analysed to confirm the preservation of conserved residues throughout the alignment and motifs, such as DR [Y/F] at the C-terminal end of TM3 and NPXXY at TM7 (Figure 1). The length of TM helices were also equivalent and TM helices were found to be highly conserved between the two receptor sequences, but the loop regions were more variable.
The cysteine residues at TM3 and ECL2 (extracellular Loop 2), thought to be involved in the formation of a disulphide bond, were conserved between the template and the query (Figure 1). The model structure, after loop building and energy minimization through SYBYL, was validated using PROCHECK. The PROCHECK results (Supplementary Figure S1) for the model, excluding the loop regions, show more than 99% of the residues are in allowed regions (97.3% in the strictly allowed region). The full-length structure (model with loops) showed more than 96% of the residues in the allowed region of the Ramachandran plot (Supplementary Figure S1). The conformation of the final model (Figure 2) was also compatible to retain the disulphide bond formed between the Cys residues of TM3 and ECL2.
Comparison of 5-HT 2A receptor model with β 2 -AR active and 5-HT subtype (5-HT 1E and 5-HT 2B ) structures First, we compared the important ligand-binding residues of the modelled structure with the template structure β 2 AR (PDBID: 3SN6). Though sequence identity between 5-HT 2A receptor and β 2 AR is only 34%, the ligand-binding site is highly conserved (Table 1). Structure superposition of the ligand-binding site reveals that all the important residues superimpose well (Supplementary Figure S2).
Next, we compared the 5-HT 2A receptor model with the recent crystal structures from the serotonin subtype 5-HT 1B and 5-HT 2B receptors. The equivalent residue for Ser239 is Gly221 in 5-HT 1B receptor and Ser212 in 5-HT 2B receptor. In addition, position of hydroxyl group in Ser239 is toward the binding pocket. But, in the case of 5-HT 1B structure which has equivalent residue, Ser212 hydroxyl group points towards the extracellular region (Supplementary Figure S3), suggesting that 5-HT 1B structure is in intermediate active state. In addition, the conformation of ECL2 region, which plays an important role in ligand-binding, is different in 5-HT 1B and 5-HT 2B receptor structure. Table 2 shows the residue difference at the ligand-binding site. Ser239 was earlier reported to be involved in agonist interaction (Shapiro et al., 2000). This is consistent with the finding that when the agonist binds, there is an inward movement of extracellular region of TM5, which results in displacement of intracellular part of TM5 away from the TM helices. Moreover, the conformation of ECL2, where Asp 231, is prior to TM5 in sequence, is quite different in 5-HT 2B structure.
Comparison of active and inactive state models of 5-HT 2A receptor We have also compared our active state model with the previously modelled inactive state model to observe the structural differences. In 5-HT 2A receptor, the TM helices TM3, TM5, TM6, TM7, and ECL2 (extracellular loop 2) are involved in ligand binding (Kanagarajadurai et al., 2009). In the active state of the GPCRs, TM5 and TM6 undergo major conformational changes which are influenced by changes in the binding site. The same was also observed in the superposition of active and inactive state model (Supplementary Figure S4). Major structural deviations were observed in the cytoplasmic part of the active state model and less structural differences were seen in the extracellular part of the model. The difference of side   TM3  D155, T160  D155, T116  TM5  S239  S203  TM6  N343, F339  N293, F289  TM7  Y370  Y316  ECL2 D231 -chain positioning of the ligand-binding residues such as Asp155, Ser 239, Asp231, and Asn 343 was also observed. This is due to the structural difference between active and inactive state models.

Docking of endogenous and known ligands
The docking modes and affinities of few known agonists (like DOI and Butotenine) and antagonists (like Seriptone) were quite different between active and inactive state models (Figure 3). The affinity and GScore details are shown in Figure 3(A) and (B  (Kristiansen et al., 2000b). The other agonist molecules interact in a similar manner as serotonin (Supplementary Figure S5), from our docking studies, except that the -OH group on the indole group of the ligand interacts with either Ser 239 [TM5] or Asn 343 [TM6] based on its position. In case of antagonists, from the docking studies, we observed a variety of interactions (with less conserved patterns) which could be a reason for different functional selectivity. However, the halogen, attached to the aromatic ring in antagonists, was usually anchored between TM3 and TM5 where Ser159 and Tyr160 are positioned. However, there are several residues that contribute to hydrophobic interactions.

Analysis of VS hits
The VS experiment resulted in 428 molecules which are considered as hits. The earlier reports and studies about the interactions of agonist and antagonist molecules were helpful in selecting the best lead molecules for further analysis. For example, the interactions of important residues such as Asp155, Ser159, Ser239, Asn363, and Tyr360 with the lead molecule were checked and selected.
To further analyse the interactions of these molecules in active and inactive state of the receptor, the molecules were also docked with models of both states of 5-HT 2A receptor. From the 428 molecules, the best 64 molecules were selected for further IFD analysis (Figure 4). These 64 molecules showed GScore in the range of −10.956 to −7.797 kcal/mol. The IFD poses of 64 molecules with active and inactive state of the 5-HT 2A receptor were further examined for their interactions. This analysis gave rise to 15 best molecules which could act as better lead molecules and these molecules were investigated for specific agonist-and antagonist-type interactions.
The scaffold-based clustering of 428 VS hits, along with known molecules, gave useful information about the similarity/dissimilarity of the identified hits with the known molecules. The similar scaffold molecules were grouped into clusters and cluster size describes the number of molecules in that particular cluster.
There are total of 108 clusters and the biggest cluster contains 117 molecules. Table 3 provides the clustering details of known molecules and similar VS hits and about the similarity of the best 15 molecules with known molecules. Four of these molecules (ZINC01673166, ZINC03872345, ZINC35111670, and ZINC35111667) have similar scaffold-like mescaline (agonist) and four molecules (ZINC12188847, ZINC34467228, ZINC19851433 and ZINC71896677) have similar scaffold-like known antagonists such as haloperidol, ketanserin, spiperone, and aripiprazole.

Validation of VS results
The VS hits and known ligands, along with 1000 decoy set, were docked with both active and inactive 5-HT 2A model. By using ROC calculation, enrichment studies were prepared on two data-sets such as ranking of known molecules (62 actives) and ranking of VS hits (428 molecules in top N% of the results (Table 4). The known ligands and hits selected by VS were retrieved with a good ROC value, as shown in Figure 5. In the Table 2. Comparison of ligand-binding residues in 5-HT 2A , 5-HT 2B and 5-HT 1B receptors.

IFD studies with both active and inactive conformation
A total of 15 molecules were selected based on interaction and GScore values and these molecules were docked with both modelled active and inactive conformation of 5-HT 2A receptor to predict the difference in binding pattern and interactions. Among these 15 molecules, 10 molecules interact in a manner similar to agonists in the active state. It was also found that these molecules interact with better affinity in active state rather than inactive state (Table 5). All the molecules have protonated amine  (428)   Notes: The plot shows known actives (left) shows nearly 90% of the molecules were recovered in active and lesser molecules recovered in inactive state docking. All the VS hits were recovered (right) with the good ROC value in both active and inactive state model.

which forms salt or H-bond interactions with Asp155
[TM3] and terminal functional group interacts with Asn363 [TM7]. The -OH or -OCH3 functional group in the aromatic head always forms H-bond interaction with Ser239 and Asn343, which were proved experimentally to be involved in agonist interactions. The GScore comparison of all the selected molecules with the known agonist-and antagonist-type molecule also showed that these proposed molecules have equal or better GScore (Supplementary Figure S6). The aromatic rings in the molecules are stabilized by residues Phe339, Phe340 [TM6], and Phe243 [TM5] which contribute toward π-π interactions. While these molecules were docked in the inactive state structure of 5-HT 2A receptor, the abovementioned interactions were not observed, except the interaction with Asp155. Instead of Ser239 [TM5], the backbone atom of Ser207 [TM4] forms H-bond interaction with ligand functional groups.
The conformation of the molecule itself was slightly different in such a way that the head group gets positioned in the cleft between TM3, TM4, and TM5 and none of aromatic residues contribute to π-π stacking (Supplementary Figure S7). The interaction difference was also reflected in the GScore, such that the docking of ligands to the active state showed better score then inactive state docking of this agonist-type molecules.
From the selected 15 molecules, a total of five molecules showed interaction pattern similar to antagonist or inverse agonists. These molecules were also docked with both active and inactive state conformation of 5-HT 2A  receptor structure and affinity details were recorded in the form of GScore. All the molecules have aromatic head groups with halogen attached and a protonated amine group. The protonated amine forms either salt bridge or H-bond interaction with Asp155, as observed in agonist-type interactions. None of the aromatic residues in TM5 and TM6 contribute to π-π interactions; instead, Trp151 [TM3] contributes towards interaction with the tail aromatic group. The halogen attached in the head aromatic group was anchored by the hydrophobic residues and it was also observed that the head group was docked in the binding cleft between TM3, TM4, and TM5. In addition, Asp231 from ECL2 also contributed towards H-bond interaction (Supplementary Figure S8).

Agonist-type molecule interactions
The availability of agonist-and antagonist-bound GPCR crystal structures provided insights about the interactions followed by molecules with different efficacy level. To behave like an agonist, the molecule should have a set of functional groups and specific interactions with certain amino acids in the site. Two major anchor interactions are necessary for full agonists. A common protonated amino group in the ligand formed a salt bridge with carboxyl group of Asp155 in TM3. The polar groups of catechol or similar moiety interacted with side chains of Ser 239 in TM5 (Figure 6(A)). Moreover, β-OH of agonist-formed H-bond network with Asp155 & Asn363, polar interaction prevailed between Asn343 and β-OH group in agonists and the rotameric positions of Trp336, Phe339, and Phe340 remain as in the inactive state, contributing towards π-stacking (Figure 3, serotonin interactions as an example). There were total of 10 molecules which showed similar types of interactions with the receptor. The two-dimensional interactions of agonist-type interactions were shown (Supplementary Figure S7).

Antagonist-type molecule interactions
A conserved scaffold is the hydroxyamine motif of the ligands with inverse agonist, antagonist, or full/partial agonist activities, which retains a conserved interaction with Asp113 and Asn363. The head groups are mainly aromatic moieties with a halogen group attached. These head groups, with aromatic moieties of all compounds, are anchored by strong hydrophobic interactions in the binding cleft (Figure 6(B)). The important residue, Ser239, is not accessible for antagonists to interact in the inactive form. This could be an important point that agonist interaction with Ser239 is important for initiating the activation process. From the 15 molecules, five molecules have the specified functional group and interactions with the receptor. Further, investigation of common and variant interactions between the known molecules and the final 15 selected molecules was performed with both active and inactive state model of the receptor. Table 6 shows the details of the interacting residues between known (8 agonists and 9 antagonists) and VS selected final molecules. The comparison was between known agonists and selected final set agonist-type molecules from the final set (Table 6) and the known antagonists with the selected antagonist-type molecules. We checked their interactions in both active and inactive state. Here, common residues were those which interact with the known molecules as well as the final 15 molecules and the variant residues were those interact either with only known molecules or only with 15 molecules. It was interesting to note that Asp155 [TM3] is invariably in the interacting mode with

Analysis of MD simulations
The changes in conformations and the stability of proposed model and receptor-ligand complexes were investigated through 10-ns MD simulations. The conformations obtained during 10-ns simulations for both apo form of active state and the best scored agonist and antagonist bound complexes were analysed using the simulation interaction diagram option of the Desmond program. Generally, three properties are analysed for proposed apo form of the model (1) RMSD (2) RMSF (3) protein secondary structure stability. The final RMSD variation from the initial model of C-alpha and backbone atoms shows that system is stable beyond 10 ns (Supplementary Figure S9). As the GPCR protein has its intrinsic conformational switches in the whole protein, it is quite possible that the protein model has acquired the changes in the order of 3-5 Å which are stabilized around a fixed value at the end of simulations. The RMSF calculations show that the large fluctuation is observed in ICL2 loop that connects TM5 and TM6 helices. It is also important to note that TM5 and TM6 undergo large conformational changes upon change of active state to inactive state (Audet & Bouvier, 2012).
The 10 snapshots, obtained at each 1 ns, were investigated for further changes in interactions with the receptor. The important conserved residues E/DRY [TM3] and NPXXY[TM7] which were located at the cytoplasmic ends of the transmembrane segments were analysed carefully in all the 10 structures. These conserved sequences have been proposed to be important for activation or for maintaining the receptor in the inactive state (Hofmann et al., 2009). All these residues occupy very similar positions as it was in the reference structure, except Arg 172 [TM3] which acquires a slightly different conformation at the end of simulation as seen in β 2 AR-Nb80 structure (Rasmussen et al., 2011a). The NPXXY motif Tyr 380 [TM7] moved more inside to the helix bundle (Figure 7(A)).
It was observed that TM6 moved outward by 11 Å as measured at the α-carbons of Glu 218 (red dotted lines) between the reference and the final structure (10th  Notes: The transmembrane regions are highlighted in green colour. The residues that are fully conserved across 5-HT subtypes are highlighted in yellow boxes. The 5-HT 2 subtype-specific and 5-HT 2A -specific residues are highlighted in red boxes. structure). The TM5 helix conformation was not changed very much during the simulations. These structural observations showed that the active state model is more stabilized by the above-mentioned changes.

MD analysis on receptor-ligand complexes
From the selected set of agonist-and antagonist-type molecules, the best scored two molecules were considered for 10-ns molecular simulations to study their stability of the docked complexes of both active and inactive models. For the receptor-ligand complex, four properties were considered for investigations (a) RMSD, (b) RMSF, (c) Secondary structure changes, and (d) Ligand-protein contacts. The first three features were mainly to ensure that there are no abnormal structural changes over the time. The important feature with respect to receptor-protein complexes are protein-ligand contacts. This would help to analyse the presence and absence of specific interactions with binding site residues. In case of agonist-type molecule, the 10-ns simulations show that

Sequence alignment of human 5-HT subtypes
The recent progress in designing or identifying lead molecules for a specific disorder caused by a particular subtype of serotonin receptor is quite challenging. The conserved and subtype-specific residues in the ligand-binding site would be of great help to design such lead molecules. Apart from the residues, the positions of those residues in TM helices also play an important role. The sequence alignment of all the human 5-HT subtypes explored the conserved residues across subtype and subtype-specific residues of ligand-binding site, which comprises TM3, TM5, TM6, TM7, and ECL2.This information would be helpful to design molecules for the particular subtype and in identifying the 5-HT 2A receptor-specific residues through the alignment of TM regions. From the docking studies, it was observed that certain residues from TM3, TM5, TM6, TM7, and ECL2 are important for ligand interactions. Some of those residues, namely Asp155

Conclusion
The modelling of the active state of 5-HT 2A receptor explains the conformational difference between active and inactive state recorded from the β 2 -AR crystal structures. Through VS experiments, a large number of molecules were screened and further studies on the hits lead to the selection of 15 molecules which could behave like agonist, antagonist, and inverse agonists. The docking studies of these molecules on active and inactive conformation of 5-HT 2A receptor revealed a set of residues important for ligand interactions. Enrichment studies and clustering analysis showed encouraging results for VS hits as well as the final selected molecules since 85 and 79% of known small molecule agonists/antagonists could be recognized within top 20% when docked against active and inactive state model of the receptor, respectively. Due to the inward movement of cytoplasmic part of TM5, residue Ser239 in TM5 is involved in interactions with the agonist-type molecules when the receptor model is viewed in the active state. However, in the inactive state, the opposite trend is observed, where TM5 moves back to its original position and small binding cleft is created between TM3, TM4, and TM5. All the ligand-head groups are anchored by that cleft. It was also found that certain residues are preferred for agonist and antagonist interactions. is involved in interaction. The proposed active state model and the observed protein-ligand interactions were further confirmed by MD simulations. Sequence analysis of human 5-HT subtypes helps to identify residues which are conserved across all 5-HT subtypes and residues specific to 5-HT 2 subtypes. These subtype specific ligand-interacting residue information could be useful for design or identify more ligands for specific subtypes. The current study supported the earlier findings that different pattern of interaction is attained by molecules with different efficacy level.