Homology Modelling and Docking Studies of Human a 2-Adrenergic Receptor Subtypes

Alpha2-adrenergic receptors (α2-ARs) belong to adrenergic receptor family. There are nine representative members is adrenergic family namely, α1a, α1b, α1d, α2a, α2b, α2c, β1, β2, and β3. The rationale for classification of adrenergic receptors initially was by agonist and antagonist binding characteristics, but later on they were identified as separate gene products, as reviewed by Blyund et al. [1]. α2a, α2b, α2c are expressed in various organs. The relative distribution of these receptors is as follows: Brain>spleen>kidney>aorta=lung=sk eletal muscle>heart=liver (IUPHAR database ) [2]. Specific functions are associated with each subtype: α2a-ARs are involved in regulation of cardiovascular function, blood pressure and plasma noradrenaline (norepinephrine). Peripheral α2b-ARs are responsible for causing sodium retention and vasoconstriction, thus are involved in mediating hypertension through stimulation by agonists and salt, vasoconstrictor response and may play a role in development and reproduction [3,4]. Peripheral α2c-ARs causes cold-induced vasoconstriction and are involved in responses to stress, locomotion, startle reflex and it also regulates the catecholamine release from the adrenal gland, through a feedback mechanism [5,6]. α2aand α2c-ARs are also involved in regulating norepinephrine release from the peripheral sympathetic neurons and thus act as presynaptic inhibitors [5]. Thus, α2-ARs serve as important drug targets against hypertension, cardiovascular dysfunction, regulating vasoconstriction, opiate withdrawal, alcohol addiction and as adjuvants for anesthesia during surgery [7,8]. Subtype specific drugs for β1 and α1 receptors, are available, but for α2 adrenergic receptors the efforts are on, yet not much progress has been achieved. Still, the treatment with nonspecific drugs like Clonidine, Medetomidine, and Brimonidine are in use for the treatment of hypertension, glaucoma, tumor pain, postoperative pain etc. [9]. Most of the time the treatment with α2-receptor agonists is withdrawn due to side effects such as sedation in hypertension. Thus, the availability of subtype specific ligands would be beneficial for the treatment of diseases related to α2-receptors. But the success of drug design efforts is achieved only in the presence of high resolution crystal structure. In the absence of high resolution crystal structure, Homology modeling is an attractive method to obtain the structure, and this method has been proved to be a suitable option to get atomic level resolution of protein structures [10,11]. A recent review by Costanzi et al. [12] elegantly describes the various issues related to homology modeling of GPCRs. Recent studies on homology models have demonstrated that their accuracy is comparable to the structures obtained through X-ray crystallography [13] and the ligand binding affinity is also on par with that of crystal structures [14]. But, the accuracy of homology model is solely dependent on template and when the sequence identity is less than 30% the reliability of model decreases. However, the earlier attempts of modeling GPCRs were on bacteriorhodopsin [15]. Later on, it was Bovine rhodopsin and for many years bovine Rhodopsin has been the only GPCR with experimental structural information available, and all the homology modelling efforts were focused on this structure [16,17]. Homology modelling and docking studies were earlier based on bovine Rhodopsin even though it showed lower sequence identity (21%) and lower transmembrane identity (26%) as the availability of high resolution GPCR structures was the limitation [18-25]. This was followed by the availability of other members of GPCR family, Human β2-adrenergic receptor [26], turkey β1-adrenergic receptor [27], squid Rhodopsin [28], Human adenosine A2a receptor [29], chemokine CXCR4 [30], Human Dopamine D3 receptor in complex with a D2/ D3 selective antagonist Eticlopride [31], and most recently histamine H1 (H1R) [32], M3 muscarinic acetylcholine receptor [33], Mu-opioid receptor [34], a lipid G protein-coupled receptor [35], M2 muscarinic receptor [36], Kappa opioid [37], Delta opioid [38], Neurotensin receptor 1 [39], chemokine CXCR1 [40], Protease activated receptor 1 [41], 5-hydroxytryptamine 1b [42], and 5-hydroxytryptamine 2b


Introduction
Alpha2-adrenergic receptors (α2-ARs) belong to adrenergic receptor family. There are nine representative members is adrenergic family namely, α1a, α1b, α1d, α2a, α2b, α2c, β1, β2, and β3. The rationale for classification of adrenergic receptors initially was by agonist and antagonist binding characteristics, but later on they were identified as separate gene products, as reviewed by Blyund et al. [1]. α2a, α2b, α2c are expressed in various organs. The relative distribution of these receptors is as follows: Brain>spleen>kidney>aorta=lung=sk eletal muscle>heart=liver (IUPHAR database ) [2]. Specific functions are associated with each subtype: α2a-ARs are involved in regulation of cardiovascular function, blood pressure and plasma noradrenaline (norepinephrine). Peripheral α2b-ARs are responsible for causing sodium retention and vasoconstriction, thus are involved in mediating hypertension through stimulation by agonists and salt, vasoconstrictor response and may play a role in development and reproduction [3,4]. Peripheral α2c-ARs causes cold-induced vasoconstriction and are involved in responses to stress, locomotion, startle reflex and it also regulates the catecholamine release from the adrenal gland, through a feedback mechanism [5,6]. α2a-and α2c-ARs are also involved in regulating norepinephrine release from the peripheral sympathetic neurons and thus act as presynaptic inhibitors [5]. Thus, α2-ARs serve as important drug targets against hypertension, cardiovascular dysfunction, regulating vasoconstriction, opiate withdrawal, alcohol addiction and as adjuvants for anesthesia during surgery [7,8]. Subtype specific drugs for β1 and α1 receptors, are available, but for α2 adrenergic receptors the efforts are on, yet not much progress has been achieved. Still, the treatment with nonspecific drugs like Clonidine, Medetomidine, and Brimonidine are in use for the treatment of hypertension, glaucoma, tumor pain, postoperative pain etc. [9]. Most of the time the treatment with α2-receptor agonists is withdrawn due to side effects such as sedation in hypertension. Thus, the availability of subtype specific ligands would be beneficial for the treatment of diseases related to α2-receptors. But the success of drug design efforts is achieved only in the presence of high resolution crystal structure. In the absence of high resolution crystal structure, Homology modeling is an attractive method to obtain the structure, and this method has been proved to be a suitable option to get atomic level resolution of protein structures [10,11]. A recent review by Costanzi et al. [12] elegantly describes the various issues related to homology modeling of GPCRs. Recent studies on homology models have demonstrated that their accuracy is comparable to the structures obtained through X-ray crystallography [13] and the ligand binding affinity is also on par with that of crystal structures [14]. But, the accuracy of homology model is solely dependent on template and when the sequence identity is less than 30% the reliability of model decreases. However, the earlier attempts of modeling GPCRs were on bacteriorhodopsin [15]. Later on, it was Bovine rhodopsin and for many years bovine Rhodopsin has been the only GPCR with experimental structural information available, and all the homology modelling efforts were focused on this structure [16,17]. Homology modelling and docking studies were earlier based on bovine Rhodopsin even though it showed lower sequence identity (21%) and lower transmembrane identity (26%) as the availability of high resolution GPCR structures was the limitation [18][19][20][21][22][23][24][25]. This was followed by the availability of other members of GPCR family, Human β2-adrenergic receptor [26], turkey β1-adrenergic receptor [27], squid Rhodopsin [28], Human adenosine A2a receptor [29], chemokine CXCR4 [30], Human Dopamine D3 receptor in complex with a D2/ D3 selective antagonist Eticlopride [31], and most recently histamine H1 (H1R) [32], M3 muscarinic acetylcholine receptor [33], Mu-opioid receptor [34], a lipid G protein-coupled receptor [35], M2 muscarinic receptor [36], Kappa opioid [37], Delta opioid [38], Neurotensin receptor 1 [39], chemokine CXCR1 [40], Protease activated receptor 1 [41], 5-hydroxytryptamine 1b [42], and 5-hydroxytryptamine 2b [43]. Recently, the modelling groups have used β2-adrenergic receptor as a template to model subtypes of α-adrenergic receptors, as it shares higher sequence identity (29-31%) and higher transmembrane identity (37-43%) with α2-ARs [44][45][46].

Homology modelling
The sequences of Human α2A (PO8913), α2B (P18089), and α2C (P18825) adrenergic receptors were retrieved from Swiss-Prot database (http://www.uniprot.org/). The template structure and sequence of Human Dopamine D3 receptor X-ray structure complexed with a selective antagonist Eticlopride (3PBL, 2.89 Å) was downloaded from PDB (www.rcsb.org). The sequence alignment was generated with Clustal W multiple sequence alignment tool ( Figure  S1) [47] and manually adjusted to avoid insertions and deletions in the transmembrane regions. Modelling was done using Modeller 9v8 [10].
The modeling parameters were set to output a single best structure with the least DOPE score. The water molecules could not be included for modelling as it is homology modelling method. The final models were energy minimized using the Steepest Descent minimization method until convergence is attained using the CHARMm Force field and submitted to PSVS software suite [48] for stereochemical validation [49]. The analysis output includes constraint analysis, goodness-of-fit and structure quality scores using information from prior knowledge and verifies the normality of torsion angles, bond angles, bond lengths and distances between unbounded neighbor atoms. The analysis measures are both global and site-and specific. Ballesteros-Weinstein convention was used to assign residue positions throughout our analysis for transmembrane helices [50] as well as loop regions [22].

Preparation of receptor and molecular volume calculations
The molecular volume of models of α2-ARs was determined using Discovery Studio 3.1 (Molecular Attributes, Accelrys Inc). The molecular volumes of the binding site cavity were obtained using calculate binding site from the receptor Cavities tool in DS 3.1. We selected the eraser algorithm, provided by the software, to detect probable binding site residues in the receptor. This algorithm detects binding site cavities based on the receptor shape Venkatachalam et al. [51].

Ligand selection
Twenty four ligands (agonist ( Figure 1) and antagonist (Figure 2)), were selected from the literature [22,46,52,53]. These ligands have been reported by the research studies to bind to and interact with the α2-AR. A set of eight agonists including endogenous ligands such as Dopamine, Adrenaline, known α2 agonists such as Clonidine, Dexmedetomidine and subtype selective agonist Guanfacine were selected from literature for docking studies [23,39,45]. Sixteen antagonists including subtype specific ligands such as BRL-44408, JP-1302, OPC-2836 and others such as ARC239, Clozapine, WB4101 that have been reported to bind to α2 adrenergic receptors were chosen for the docking studies to analyse their interactions with the α2-AR models based on Dopamine receptor (PDB ID:3PBL) [22,46,52]. The 3D structures were obtained from small molecule database Pubchem, present in NCBI server (http://pubchem.ncbi.nlm.nih.gov/).

Docking studies
Ligand preparation: All the agonist and antagonist structures were prepared using LigPrep (Schrödinger, LLC, New York, USA, 2012). The ligands were checked for the three-dimensional (3D) structure, realistic bond lengths and bond angles, covalent bonds, accompanying fragments, hydrogens, protonation states and were prepared using Ligprep. The LigPrep process consists of a series of steps that perform conversions, apply corrections to the structures, generate variations on the structures, eliminate unwanted structures, and optimize the structures. To initiate the process, the ligand structures in sdf format were used as input for the LigPrep module. The agonist and antagonist structures were prepared separately. We used the default OPLS_2005 as the Force field for minimization and the ionization states to be generated were set to the default p H of 7.0 ± 2.5. Ionization states were generated using Epik which generates these states in aqueous solution. We also selected the option to retain the original state of the input molecule. To ensure that other molecules such as water and counter ions are excluded from the ligand structure, we used the Desalt option in LigPrep. We selected the Generate Tautomer option to generate tautomeric forms of the input ligands. Finally, we used the Retain specified chiralities option to use the chiral information from the original ligand file after the module fixes any irregularities with regard to atom numbering and bond directions. The prepared and minimized ligands are generated as .mae outfile which is used for docking.
Protein structure preparation: The α2-adrenergic receptor models were prepared using protein preparation wizard of Schrödinger Maestro 9.3 (Schrödinger, LLC, New York, USA, 2012). Schrödinger offers a comprehensive protein preparation facility in the Protein Preparation Wizard, which is designed to ensure chemical correctness and to optimize protein. The protein structure pdb file was used as input for this module. In the process tab, we selected the options Assign bond orders, Add hydrogens, Create zero-order bonds to metals, create disulphide bonds, Fill in missing side chains using Prime and Delete waters beyond 5Å from hetero groups options to fix any irregularities in the structure. Finally we refined through restrained minimization using the OPLS 2005 Force field. This minimization is performed by the impref utility which uses Impact to perform minimization wherein a harmonic potential of 25 kcal mol -1 Å -2 is used to restrain the heavy atoms.
Ligand docking: The ligands were docked using the Glide application which uses Impact version v18007 program to perform the docking (Schrodinger, LLC, New York, USA, 2012) [54][55][56][57]. It performs grid-based ligand docking with energetics and searches for favorable interactions between one or more typical small ligand molecules and a typically larger receptor molecule, usually a protein [54]. The receptor grid generated and the prepared ligands were used as input. The ligand sampling was set to flexible which generates alternate conformations of the input ligands internally while docking and torsional sampling bias for nitrogen inversions, ring conformations and amide bonds. Since Epik was used for ionization during ligand preparation, Epik penalties were imposed during docking to compensate for higher The ligands were docked with the active site using the 'extra precision' Glide algorithm. The default van der Waals scaling factor of 0.8 with a partial charge cutoff of 0.15 were used to reduce penalty for close contacts for non-polar ligand atoms in the ligand tab of the glide panel. Glide generates conformations internally and passes these through a series of filters. It first places the ligand centre at various grid positions of a 1Å grid and rotates it around the three Euler angles. At this stage, crude score values and geometrical filters weed out unlikely binding modes. The next filter stage involves a grid-based force field evaluation and refinement of docking solutions including torsional and rigid body movements of the ligand. The final energy evaluation is done with GlideScore and a single best pose is generated as the output for a particular ligand. The generated poses were analyzed to observe their orientations and interactions.

Re-docking of co-crystallized ligands:
To check the validity of the method, Eticlopride complexed with 3PBL was used and re-docked using the protocol described above for Glide docking. Each re-docked pose was evaluated by considering the RMSDs and docking scores. The selected re-docked pose was further evaluated by its interactions and energetic analysis to investigate the efficiency of the docking search algorithm and scoring function by comparing its values with the bound conformation.

Homology modelling
Comparative modelling was used to model α2-adrenergic receptors of Human based on the X-Ray crystal structure of Human Dopamine D3 receptor (PDB ID: 3PBL) as they share highest sequence identity and transmembrane identity (α2a: 34%, 49%; α2b: 32%, 49%; α2c: 34%, 49%), in comparison to available crystal structures ( Table 1). The sequences of α2-ARs and the template 3PBL were aligned using Clustal W and then were manually adjusted to avoid insertions and deletions in the transmembrane domains. The transmembrane regions were aligned according to Baldwin's model for the alpha carbon positions in the transmembrane helices of the GPCR family [55]. In all three α2-ARs, the N, C-terminal part and the third intracellular loop (IL3), which is over 100 residues long, were not modelled because the available loop modelling algorithms are limited to up to 13 residues long loops [45].
The homology modelling package Modeller (version 9v8) [10] was used to generate the three-dimensional models for α2-ARs based on the X-ray structure of Dopamine Receptor using the sequence alignment presented in Figure S1. The resulting models were first geometrically refined in order to reduce the side chain steric clashes. Later on, the entire receptor was energetically minimized using steepest descent method and conjugate gradient using CHARMm Force field. The final structures were validated using PSVS validation suite [48]. The final models have over 90% of residues in the favorable regions of the Ramachandran map ( Figure S2) and all the main-chain parameters, like peptide bond planarity, bad nonbonded interactions, C-α distortion, overall G-factor, bond length distribution and side-chain parameters, are in the normal range.
Ramachandran Plot for α2a-AR showed 85.7% residues the most favored regions and 14.3% in the additionally allowed region. Verify 3D analysis showed a score of 0.21, a Zscore of-4.01 and the MolProbity Clash Zscore was-2.10.
Ramachandran Plot for α2b-AR showed 85.2% residues the most favoured region, 14% in the additionally allowed region, 0.4% in the generously allowed region and 0.4% in the disallowed region. The Verify 3D score was 0.21, the Zscore for this model was-4.01 and the MolProbity Clash Zscore was-2.70.
Ramachandran Plot for α2c-AR showed 80.6% residues the most favoured region, 18.6% in the additionally allowed region and 0.8% in the generously allowed region. The Verify 3D scaore was 0.24, the Zscore for this model was-3.53 and the MolProbity Clash Zscore was-1.77.
The Zscores for all three model structures are greater than-4, indicating that the structures are properly refined and can be used for further analysis.

Binding site cavity analysis
Binding site cavity volume was obtained by binding site cavity detection tool of DS 3.1. Our analysis showed largest binding cavity volume in α2b-AR (503 Å 3 ), followed by α2c-AR (471 Å 3 ) and the smallest cavity was observed in α2a-AR (403 Å 3 ). The results are in agreement with previous studies, where the binding cavity predicted was largest for α2b-AR [25].
Docking results α2a-AR: Agonists: Of the docked agonists, the Dopamine was observed to have higher glide score (-7.979107), followed closely by Adrenaline (-7.876917). The other agonists had a glide scores in the range-6.862784 to-4.244504. Dopamine was observed to interact through hydrogen bonding with the residues Asp 3.32 (TM3), Ser 5.46 (TM5) ( Figure 3) and was also observed to interact through other electrostatic interactions with the residues Ser 5.42 (TM5), Tyr 7.43 (TM7). Some of the agonists such as Adrenaline and 2-amino-1phenylethanol were also involved in π-π interactions with Phe 6.51 and Tyr 6.55. The details of the interactions of the other ligands are given in Supplementary table S1.

α2b-AR:
Agonists: Guanfacine was observed to have the highest binding affinity with a glide score of-7.533895. The endogenous ligands Adrenaline and Dopamine followed closely with glide scores of-6.899865 and-6.858532 respectively. Most of the agonists were observed to interact with the α2b receptor through hydrogen bonding with Asp 3.32 (Supplementary table S1). Phe 6.51 and Tyr 6.55 were involved in π-π interactions with some of the agonists (Supplementary table S1).
Antagonists: WB-4101 was found to have the highest glide score (-8.641979) and was involved in interactions with Asp 3.32 (hydrogen bond), Lys xl2.51 (π-cation) and Phe 6.51 (π -π). The other antagonists were observed to have a glide score in the range-8.592746 to-2.15852. Many of the antagonists were observed to be involved in interaction with the residue Asp 3.32 through hydrogen bond and with Phe 6.51 through π -π interactions (Supplementary table S1).
α2c-AR: Agonists: Dopamine was found to have the highest binding affinity of the agonists used for docking in our study with a glide score of-9.078919 (Supplementary table S1). It was found to interact with the residue Asp 3.32 through Hydrogen bond formation (Figure 3). It was also involved in Hydrogen bonding with Ser 5.41 (Supplementary table S1). Most of the other agonists such as Adrenaline, 2-amino-1phenylethanol were also observed to interact with Asp D3.32 through hydrogen bonding and π-π interactions with the residues Phe 6.51, Tyr 6.55, Phe 7.38 and Phe 7.39 (Supplementary table S1).

5.3.6) Antagonists:
We observed that of the sixteen antagonists used for docking with α2c receptor, WB-4101 had the highest glide gscore (-9.87093) and showed π-π with Phe 6.51. The other antagonists were bound to the α2c receptor with Glide scores in the range-8.176154 to-0.558762 (Supplementary table S1). The antagonists BRL-44408 and Atipamezole showed hydrogen bonding interaction with Asp 3.32 while Prazosin showed hydrogen bonding with Tyr 3.28 (Supplementary table S1). Most of the antagonists were involved in π-π interactions with Phe 6.51, few with Phe 7.39 and Phe 5.48 (Supplementary table  S1).

The orientation of common residues in the models and the prediction of their importantance in residue ligand interactions
Most of the residue, that were predicted facing the ligand binding cavity by earlier investigators, based on bovine Rhodopsin [22,25] and β2-adrenergic receptor [44][45][46][47] models, were occupying the same position in the models based on human D3 Dopamine receptor also. The amino acid residues predicted to be the key determinants of agonist binding, D3.32, W6.48, S5.43 and S5.46, were well positioned to interact with the protonated nitrogen, aromatic ring, and catecholic hydroxyls of catecholamine agonists (Figure 3). The position of the binding site was comparable to what was observed by previous studies, which is a validation of experimental and theoretical findings by different groups [18,21,22,25,[44][45][46][47]54]. Furthermore, as described before, the binding pockets were mainly formed by conserved hydrophobic amino acids of TM5 (F5.47) and TM6 (W6.48, F6.51, F6.52 and Y6.55). The aromatic rings of the hydrophobic side chains of these residues may interact with the aromatic rings (or other hydrophobic structures) of ligands through π-π stacking interaction [21,22,25]. The binding modes of the ligands in our study (Figures 3 and 4) were generally consistent with previous results [18,[20][21][22][56][57][58][59] with some differences.However, the model based on the Human Dopamine D3 receptor (3PBL) was useful in identifying residues which may be important in conferring subtype specificity.
From our analysis of binding site residues, we observed that V (2.53) and V (2.57) were not facing ligand binding cavity and hence could not influence ligand binding as predicted [22]. F (3.35) [18,21,22] was projecting away from binding site and hence may not have any role in binding. This observation is in agreement with Frang et al. [60]. This group proved, by experimental studies, that F (3.35) in α2a is not exposed in the binding pocket and thus not accessible for Phenoxybenzene and other receptor ligands. However, the orientation of C (3.36), important in ligand binding and earlier predicted by Frang et al. [60], was towards binding site proving its involvement in ligand binding. This observation supports experimental finding [60] and theoretical finding [22] hence proving that our models are suitable for further studies.
The other residues, T (3.37) [18,[20][21][22] and I (3.40) [22], also may not be able to make any interaction as they are pointed downwards below the binding cavity. Nonetheless, we suggest that hydrogen bonding interaction between S (5.46) and T (3.37) along with van der Waals interaction with I (340), may be responsible for maintaining the receptor in inactive state as reported earlier [61] (Figure S3). We observed this interaction, in the receptor models of α2a, α2b and α2c adrenergic receptors ( Figure S3). While in docked protein models, the distance between S (5.46) and T (3.37) and the orientation of the bonding residues were different, fitting well with the accepted mechanism of receptor activation [62]. Furthermore, we observed that the distance between T (3.37) and S (5.46) in α2a is around 4.8 Å and hence the hydrogen bonding interaction was impossible between these two residues, therefore less stability to an inactive state and thus can be easily activated. However, the functional relevance of this observation along with the cavity size, in the regulation of receptor needs further investigation. With a similar reasoning we propose that in α2c, where the distance between S (5.56) and T (3.37) is 2.38 Å, the inactive state is more stabilized than α2a. However in α2b-AR the distance between T (3.37) and S (5.46) is 1.9 Å and hence the binding affinity is more, which results in the highly stabilized inactive state. Our results provide evidence of receptor activation, where it was suggested that T (3.37) interacts with TM5, stabilizing the inactive state of the receptor [61]. We would like to mention here that we observed inverse relationship between binding site cavity volume and stability of the inactive state of α2-AR as enforced by hydrogen bond formation, while I (3.40), highly conserved in the whole class A GPCR family, facilitates the reorientation of TM5. It was proposed that the structural change of TM5 during the process of GPCR activation involves a local P (5.50) induced unwinding of the helix, acting as a hinge, and the highly conserved hydrophobic I (3.40) side chain, acting as a Pivot [61]. We also agree that I (3.40) ( Figure S3) is not involved in the initial binding step but participates in the subsequent signal propagation as was observed in mutational analysis of I (3.40) in Histamine H1 receptor upon histamine binding [61]. The difference in the orientation of the residues, namely T (3.37) and S (5.46), before and after activation supports the findings of Sansuk et al. [61] that there is a minor, but significant, clockwise rotation of TM3 during the process of receptor activation.
The cationic nitrogen present in agonists formed a salt bridge with negatively charged side chain of carboxylate of D (3.32) in α2-ARs ( Figure 3). Agonists were docked close to S (5.42) and S (5.46), the residues predicted to be important for ligand binding [18,[20][21][22]. However, not all of the antagonist docking modes formed the required interaction with D (3.32), and many were shifted away or remained at the top of the binding cavity as reported earlier [22,46]. We relate this docking mode to the residues T (3.37) and S (5.46) which are responsible for the rotation of TM3 in conjunction with I (3.40) with P (5.50) at pivot [61]. This rotation of TM3 may create new space near TM2/TM7 such that antagonists would be shifted away from TM5. Moreover, the absence of interaction of D (3.32) with antagonist may be due to alternate binding mode, where in the antagonists interact with D (3.32) via carboxylate-aromatic interactions as reported earlier [22]. The position of S (5.46), which was predicted to be important for binding, appears to be away from ligand binding cavity and would not be able to form any interaction with ligands but may play important role in receptor activation [61]. F (5.47) was also unable to form any interaction with ligand as predicted before [18,20,22]. We believe that stronger interactions with D (3.32) and W (6.48) occur in α2-ARs and the hydrogen bonding network is different in α2-ARs. We propose the interaction D (3.32) or/and W (6.48) may be the original driving force during the whole activation process similar to the observation made by Gong and Wang [63]. F (6.51) [18,[20][21][22] appeared to interact while F (6.52) [18,20,22] and F (6.53) ( Figure S4) [21] were not interacting at all in α2a/b/c-AR. W (6.48) was oriented in similar position but was pointed towards ligands and may play key role in activation even though no hydrogen bonding was visible in the models. On binding, the aromatic catechol ring of catecholamines presumably has a direct interaction with the aromatic residues of the rotamer toggle switch, W (6.48) and F (6.52). Previous studies using Monte Carlo simulations suggested that rotamer configurations of C (6.47), W (6.48) and F (6.52), the residues that comprise the rotamer toggle switch, were involved in the movement of the cytoplasmic end of the TM6 by coupling and modulating around highly conserved Proline kink (6.50) [64]. The next residue predicted earlier was F (7.38) [18,20] but it was pointing away from ligand binding cavity in our models and hence was unable to interact with ligands in α2a /b/c-ARs. This observation is in agreement with the findings of Balogh et al. [25].

Residues which may play important role in subtype specificity
Recently, Ostopovici-Halip et al. [45] have analyzed the role of different amino acids occupying the same position around the endogenous ligand within 6 Å distance and gave a detailed explanation of the variations and their role in the binding of subtype specific ligands. Our findings were in agreement with their analysis to some extent and differed with them on the role of some residue variation in giving subtype specificity. Apart from this we have been able to identify novel residues which may be important for subtype specificity.
Our models also depict that the variation at positions 2.57 (Val86/ Ile65/Val104) and 5.39 (Val197/ Ile173/ Ile211) in α2a, α2b and α2c, would not make a big difference in the properties of the binding site as explained by earlier investigators [45,46]. We agree with the earlier studies in terms of the similarity and differences shared by valine and isoleucine with regard to topology, physicochemical properties, size and volume. Though the presence of valine residue would allow for ligands with larger substituent, ligands that can distinguish between the two residues are thought to be rare [45,46]. However, the orientation of residue V/I (5.39) in α2a-AR and 2b-AR was towards center of binding cavity whereas the I (5.39) in α2c was facing outside the binding cavity in our models, which can play key role in subtype specificity. This observation disagrees with Laurila et al. [46] and Ostopovici-Halip et al. [45], who have predicted that this change may not make a difference in binding properties.
From analysis of our models, we observed that Y (5.38) [43] was positioned between TM5 and TM4 in α2a and α2b whereas it is positioned within the binding cavity in α2c. Even though its role in ligand binding was predicted earlier, its role in subtype specificity was not suggested before. We propose that the orientation of this residue should be considered in designing subtype specific drugs ( Figure S5).
The next residue which may contribute to subtype specificity is C (5.43)/S (5.43)/C (5.43) [45]. Our models show that the cysteine (5.43) in α2a and 2c would not interact with the ligands but S (5.43) in α2b is close to ligand binding site. The hydroxyl group in the serine sidechain could be involved in the formation of hydrogen bonds with the protein backbone or with diverse polar functional groups from ligands ( Figure S6), confirming earlier observation of Ostopovici-Halip et al. [45]. The orientation of S (5.43) and its position, close to agonist, in our models also confirms the findings of Balogh et al. [25]. As described by Ostopovici-Halip et al. [45], due to its position C (5.43) in α2a-AR and α2c-AR, cannot engage in disulfide bond formation with the other cysteines. Thus, this residue confers subtype selectivity which can be exploited during ligand design through incorporation of a substituent that can interact with the side chain of S (5.43) for α2b-AR and C (5.43) for α2a-AR and α2c-AR. Balogh et al. [25] also reported the role of this residue in conferring subtype selectivity during agonist binding.
The binding site residues in xl2 loop at the top of binding cavity, which we predict could give rise to subtype specificity, are E (xl2.51)/K (xl2.51)/G (xl2.51). Thompson et al. [47] and Laurila et al. [59] have reported that this residue along with xl 2.50 and xl 2.52 acts as a lid covering the binding cavity and may interact with certain ligands to influence the binding mode. Our models show that the position of lysine in α2b-AR (K (xl2.51)) was at the top of binding cavity, going across the cavity. N (xl2.53) appears to be the major player along with E (xl2.51)/K (xl2.51)/G (xl2.51) influencing ligand entry ( Figure 5). The variations in xl2, described by previous investigators, which could give rise to subtype specificity, were K174/D153/R192 (xl2), I190/L166/ L204 (xl2.52), and D192/Q168/ D206 (xl2.55) in α2a-AR, α2b-AR and α2c-AR respectively [45]. Our models show that the orientation of K174/D153/R192 (xl2) and D1192/Q168/D206 (xl2.55) (Figures 6b,  6c and 6d) were pointed upwards and hence may not interact with ligands directly but may play role in the entry of ligand as they are close to binding site cavity, and interaction with legends before they enter binding cavity would play crucial role in ligand recognition specific to subtype α2a and signalling ability of a ligand [65]. However, the other residue at position xl2.52, (I190/L166/L204) pointed downwards, may give rise to subtype specific binding as suggested by Ostopovici-Halip et al. [45] (Figures 6b, 6c, 6d). As stated by Ostopovici-Halip et al. [45], the negatively charged carboxyl side chain of Asp192 (α2a-AR) and Asp 206 (α2c-AR) and the polar side chain of Gln 168 (α2b-AR) could be used in designing ligands with substituent of opposite charge that can interact with these residue side chains. We suggest that the prolines may have important role in giving structural stability to secondary structural elements. Thus, we propose that mutation of proline at (xl2.48) may alter loop architecture of xl2 and may alter ligand binding affinity of subtype b. Recently, single proline was reported to be associated with some architectural pattern rather than longer prolines [66].
The cysteine at position (xl2.50) has similar roles in the β 2adrenoceptor and in rhodopsin, forming a disulphide bond with the cysteine at position 3.25 in TM3 and constraining xl2, to fold on top of the binding cavity. The other disulphide bridge was observed between cysteines in xl3 region only in α2b (397-401) and α2c (408-412) adrenergic receptors. In α2a one of the cysteine residues in xl3 was replaced by glycine and hence disulphide bridge could not be established. We relate substitution of cysteine by glycine in Human α2a receptor a process of creating specificity to receptor, thus giving more space for the entry of ligands. The functional correlation of this structural feature needs further evaluation.
From our modelling and docking studies we show that residue variation at R405/H404/G416 (7.32), present in the beginning of seventh helix, would play key role in subtype specificity. R (7.32) in α2a-AR was observed to interact with antagonist docked in the top of binding cavity through hydrogen bond formation as it is positioned close to TM7 ( Figure S7). Thus, we believe that R (7.32) could be related to subtype specificity of α2a-AR.. Earlier, Laurila et al. [59] have also suggested that R7.32 could indirectly influence ligand binding through interaction with xl2. Bokoch et al. [67] reported that a salt bridge is formed between Lys305 (7.32) and Asp 192 (xl2) in β2-adrenoceptor and this bridge is displaced by TM6 during activation of the receptor.
Ostopovici-Halip et al. [45], in their molecular docking experiments with JP-1302, a selective antagonist of the α2c-AR, have mentioned that presence of glycine (7.32) in α2c-AR subtype allows for the accommodation of the ligand's acridine ring into a hydrophobic pocket located in the extracellular part of the receptor, between the upper parts of helices 6 and 7. This position is occupied by larger residues-arginine, in the case of α2a-AR and histidine in the case of α2b-AR, which obstruct the acridine ring, and implicitly the entire ligand, to adopt a similar orientation as in the α2c-AR binding site.
The next residue which was identified from molecular models based on Human dopamine D3 with D2/D3 antagonist was F (7.39). Even though its role in ligand binding was predicted earlier [18,[20][21][22], our models show that its position within the binding cavity may determine the binding cavity volume of the associated subtype. The orientation of this residue is same in α2a and α2c while in α2b it is pointed towards the membrane ( Figure S8). G (7.42), predicted by Gentili et al. [21], showed no interaction with ligands in all the α2-AR. However, the orientation of G (7.42) was towards ligand binding cavity in α2a-AR and outside in α2b and 2c ( Figure S4). Following residue was Y (7.43), which was penetrated deep into the ligand binding cavity and was in close contact with the ligands even though no hydrogen bonds were seen. Tyrosine at this position may interact with aromatic rings of ligands through π-π stacking interactions as reported by Gentili et al. [21] and Xhaard et al. [22].

Docking studies with subtype specific ligands
We docked subtype specific drug BRL-44408, with the models of α2-ARs. BRL-44408 docking was reported earlier in the model based on β2 adrenergic receptor [45]. BRL-44408 is stabilized by cation-π interactions with K (175) (xl2), whereas, K175 is stabilized by hydrogen bond with E (189) (xl2) (Figure 6a). The residues identified for binding of BRL-44408 by earlier study [45], K174/ D153/ R192 (xl2), are present on the top of binding cavity of the models and may not play any role in ligand binding directly, but interactions between K174/D192/I109 (xl2) in α2a, D153/Q168/L166 (xl2) in α2b and R192/D206/L204 (xl2) in α2c can be attributed to creating subtype specificity as their interactions may enclose a characteristic space. The specificity of BRL-44408 was further confirmed by higher glide score for α2a-AR (Table S1). A more negative value indicates a more favorable binding. Glide score is an empirical scoring function that approximates the ligand binding free energy. It has many terms, including force field (electrostatic, van der Waals) contributions and terms rewarding or penalizing interactions known to influence ligand binding. It has been optimized for docking accuracy, database enrichment, and binding affinity prediction [54,68,69]. Docking of α2c-receptor with JP-1302 resulted in higher glide score for α2c-AR compared to α2a-AR and α2b-AR ( Figure 7) in agreement with Ostopovici et al. 44] but we suggest that the residue F (7.35) involved in π-π stacking interactions may stabilize the ligands [22] and these interactions may be influenced by the side chains of the residue R/H/G at 7.32 giving subtype specificity to JP-1302. We have identified F (7.35) to be important in ligand binding in our previous paper [69]. However, the amino acid residue variation described by Ostopovici-Halip et al. [45] indicated for subtype specificity, was not consistent with our observation from our models. We suggest that the orientation of T/G/Y (6.58) and K/Q/K (7.36) may not play important role in giving rise to subtype specificity to antagonist as they are pointed away from the docked position of the ligand.
The method of docking (Glide) was tested by docking eticlopride (antagonist) to crystal structure of Human dopamine D3 Receptor. The results obtained reproduced the interactions reported in the original study by Chien et al. [31].

Discussion
GPCR crystallization is a challenging task by itself as they are unstable outside the membrane and adopt many conformational states and besides all this, loops add to the structural diversity. The fact that they are involved in many physiological processes, make them important targets in GPCR superfamily of proteins [1,70]. However, the absence of high resolution crystal structure of α2-AR subtypes is a limitation in understanding atomic details. We have modelled α2-adrenergic receptor subtypes using Human Dopamine D3 receptor (3PBL) as template. Our models reproduced most of the key interactions reported by experimental and theoretical investigators [18,[20][21][22]58,59] with some novel findings.
The first important finding was with respect to binding site cavity volume. While the binding site volume in our models was largest in α2b-AR, as predicted by previous studies based on the crystal structure of Rhodopsin and β2 adrenergic receptor, the smallest cavity was observed in α2a-AR, instead of α2c-AR, as reported by previous results. Another important observation was comparable binding site volume of α2b-and 2c-AR. This observation can be correlated to the shared specificity of ligands between α2b and α2c subtype. We propose that small binding site cavity volume could be the reason for the observed 10-100 fold lower affinity binding of bulky antagonist with an extended side chain to α2a in comparison to α2b and α2c. This experimental observation could not be explained by studies based on previous models which concentrated on differences in TM5 and xl2 and involvement of TM1 [45]. From our analysis of binding cavity residues it was observed that, the models agree that residues D (3. The next important result from our studies was the finding that the residues T (3.37) [18,[20][21][22] and I (3.40) [22] may not be able to participate in ligand binding as predicted by previous investigators. We believe that interaction between T (3.37) and S (5.46) may play an important role in inactive state stabilization along with I (3.40) in signal propagation [59]. We extended our study to analyze the residues predicted to play a role in subtype specificity. Here, our models support the findings of Xhaard et al. [22] and Ostopovici-Halip et al. [45] that variation in position 2.57 (Vl86/I65/V104), 5.39 (V197/ I173/ I211), xl2.52 (I190/L166/L204) would not affect the properties of binding sites. Our models also support earlier findings that the residue at position 5.43, C(5.43)/S(5.43)/C(5.43), could be considered for subtype specificity as described by Ostopovici-Halip et al [45].
From our modelling studies we have observed that the residues E (xl2.51)/K (xl2.51)/G (xl2.51) and N (xl2.53) together play key role in subtype specificity, as their interactions would influence the available space for ligand binding. Mutational analysis of N (xl2.53) will elucidate the role of these residues in subtype specificity. The length and the charge of Lysine may influence the ligand entry and hence create specificity for α2b-subtype. Hydrogen bonding interactions in xl2 may play important role in stabilizing xl2 loop and influence the affinity of ligand binding by acting as important constraint, thereby limiting the availability of binding cavity for ligand binding.
In our models also such interaction between K(175)-E189 (xl2.51) in α2a-AR and between Q (193)-Q201(xl2.54) in α2c and between Q(154)-K165 (xl2.51) in α2b-AR was observed ( Figure S9). It was observed that the hydrogen bond formed between Q154-K165 (xl2.51) in α2b was missing after docking. The functional relevance of this change needs to be correlated further but this interaction appears to play important role between active and inactive state. Studies on GPCRs have indicated the functional roles of loops in receptor activation and ligand binding [71]. The second extracellular loop was reported to be important in ligand selectivity in aminergic and other small moleculebinding GPCRs as described in a review [72]. We propose that mutation of proline at (xl2.48) may alter loop architecture of xl2 and may alter ligand binding affinity of subtype b. Constraining the loop seems to be essential for receptor activation among all class-A GPCRs because disturbance of the conserved disulfide bridge between xl2 and TM3 largely diminishes receptor function [73]. As stated earlier, xl2 is an important determinant in the subtype selectivity of ligands. The conformation that xl2 can adopt to accommodate these large side chains might be receptor-specific and could be used in the design of subtype-selective ligands [72].
R405/H404/G416 (7.32) is present in the beginning of seventh helix. We have observed that variation at 7.32 may interfere with interaction between ligands and F (3.35) and this interaction was identified by our modelling studies. F (3.35) was identified to stabilize ligands with π-π stacking interactions in docking of JP-1302, an α2c specific antagonist. Our models have also shown that R (7.32) could make specific interactions with antagonist at the top of binding cavity in α2a-AR as reported earlier [46].
Our models have identified the key interactions of BRL-44408, an α2a specific antagonist with K 175 (xl2). Furthermore, we believe that these interactions are stabilized with the formation of hydrogen bond with E189. Further analysis of K175 (xl2) and E189 (xl2) is required for understanding its role in subtype specificity.
The third extracellular loop has been proposed to be important in GPCR signaling [74]. In Human α2a-AR one of the cysteine residues in xl3 was replaced by glycine and hence disulfide bridge could not be established. We relate substitution of cysteine by glycine in α2a-AR receptor a part of creating specificity to receptor, thus giving more space for the entry of ligands. The functional correlation of this structural feature needs further evaluation. However, there are evidences based on mutagenesis studies, in which extracellular non conserved cysteine molecules were predicted to be important for other aspects like kinetics of ligand binding [75].
We believe that comparable size of the binding cavity of α2b-AR (503 Å 3 ) and α2c-AR (471 Å 3 ) is one of the factors influencing the similar binding affinity of ligands studied. The funnel shaped geometry of binding cavity was observed as observed by Balogh et al. [25].

Conclusion
Thus, in our study, we obtained Human alpha2-adrenergic receptors (α2a, α2b and α2c) homology models, based on crystal structure of Human Dopamine D3 receptor as it showed highest transmembrane identity in comparison to available crystal structures. We suggest that these models would prove useful in structure based drug design studies, as they agreed with experimental findings regarding the role of residues important for binding and showed correct orientation of the conserved residues involved in binding as reported before. Based on this, we suggest that the predictions for the residues critical for subtype specificity may be important.