A Comparison in the Use of the Crystallographic Structure of the Human A 1 or the A 2 A Adenosine Receptors as a Template for the Construction of a Homology Model of the A 3 Subtype

In the last decades, the field of therapeutic application in targeting the human A3 adenosine receptor has represented a rapidly growing area of research in adenosine field. Both agonists and antagonists have been described to have a potential application in the treatment of several diseases, including, for example, glaucoma, cancer, and autoimmune inflammations. To date, the most severe factor limiting the accuracy of the structure-based molecular modeling approaches is the fact that the three-dimensional human A3 structure has not yet been solved. However, the crystallographic structures of either human A1 or A2A subtypes are available as potential templates for the construction of its homology model. In this study, we have compared the propensity of both models to accommodate a series of known potent and selective human A3 agonists and antagonists. As described, on the basis of the results obtained from this preliminary study, it is possible to affirm that the human A3 receptor model based on the crystallographic structure of the A1 subtype can represent a valid alternative to the one conventionally used today, based on the available A2A structures.


Introduction
Adenosine is a key extracellular signaling molecule that regulates the cellular responses to tissue damage, hypoxia, and energy depletion, through activation of G protein-coupled receptors [1].To date, four adenosine receptor (AR) subtypes have been identified and pharmacologically characterized: A 1 A 2A , A 2B , and A 3 .These receptors are ubiquitously expressed on almost all cell types [1].Retrospectively, the human A 3 AR subtype (hA 3 AR) was the last member of the adenosine family to have been cloned.It was originally described in 1991 by Meyerhof and collaborators as an orphan receptor from rat testis and coded as tgpcr1, sharing 40% of sequence similarity with canine A 1 and A 2A ARs [2].One year later, Zhou and collaborators described the cDNA sequence, initially named R226 and extracted from the rat striatum, that encoded for a G protein-coupled receptor with an identical sequence of tgpcr1 and able to bind adenosine [3].This experimental evidence led to the conclusion that it was a new AR subtype, namely A 3 .In the last decades, the field of the therapeutic application in targeting the hA 3 AR represents a rapidly growing area of research in adenosine field, both agonists and antagonists have been described to have a potential application in the treatment of several diseases, e.g.glaucoma, cancer, and autoimmune inflammations [4].In particular, hA 3 AR agonists have been recently described as promising antinociceptive agents in different preclinical models of chronic pain [5] and in clinical trials as drug candidates for treating psoriasis, rheumatoid arthritis, and hepatocellular carcinoma [6].In addition, hA 3 AR antagonists were described as potential drug candidates for the treatment of respiratory tract inflammations, such as asthma [7] and glaucoma [8].
As recently reviewed by several authors, in the recent past structure-based molecular modeling approaches, above all molecular docking, have been increasingly applied in rationalizing the structure-activity relationships (SARs) of both agonists and antagonists and in supporting the design of novel potent and selective hA 3 AR ligands [9].More recently, molecular dynamics (MD) simulations have been performed to elucidate the ligand recognition at the hA 3 AR [10].
To date, the most severe fact that limits the accuracy of the structure-based molecular modeling approaches is due to the fact that the three-dimensional hA 3 AR structure has not been solved yet.The first fundamental breakthrough in this field is the disclosure of the first structure of the hA 2A AR, a close homolog, which allowed the construction of three-dimensional models for homology by an accuracy level much higher than before, in terms of sequence identity and resolution of the available templates [11].Following the publication of the first structure, the availability of novel crystallographic information of hA 2A AR in complexes with both agonists and antagonists are contributing to increase the knowledge at the molecular level of ligand recognition at ARs, leading, finally, to an increased applicability domain and accuracy of structure-based ligand design (SBLD) approaches (not only in the adenosine field).Nowadays, just 46 crystal structures are already available, bound to a wide variety of agonists and antagonists that are either natural or synthetic, and representing the majority of structural determinants on adenosine receptors ligands' activity.However, dealing with the hA 3 AR at a molecular level is still challenging, because the accuracy of the hA 3 AR model strongly depends on the sequence identity/similarity and resolution of the available crystallographic templates used for its construction.These bioinformatics and structural details are particularly important especially if we consider their impact in guaranteeing an accurate structural description of extra-cellular domains, for which it has already been described as having a key role in defining both the ligand orthosteric binding site and the ligand meta-binding sites [12][13][14].
As expected, several remarkable differences to the previously solved hA 2A AR structures have been underlined, in particular, the hA 1 AR revealed a peculiar conformation of the second extracellular loop and a wider ligand binding site (Figure 1).Intriguingly, from a phylogenetic and functional point of view, the hA 3 AR is much closer to the hA 1 AR than to the hA 2A AR.Indeed, hA 3 AR shows a higher sequence identity with the hA 1 AR (54%) than the hA 2A AR (49%), and both A 3 AR and A 1 AR, coupled to G i/o proteins share, from a functional point of view, the crucial biochemical pathway based on the inactivation of the adenyl cyclase activity.Structurally, another important element of similarity between the hA 3 AR and the hA 1 AR is related to both receptors having only a single cysteine residue on the second extracellular loop (EL2) which establishes a disulfide bridge, largely conserved among the Rhodopsin-family of GPCRs, along with a second cysteine located on the transmembrane helix 3 (TM3) [15,16].Moreover, the lack of the additional disulfide bridges at the EL2 level, as present in the A 2A AR subtype, involves a significant reorganization of the transmembrane helices 1, 2, 3 and 7, leading to a different shape of the orthosteric binding cavity.An additional disulfide, C260-C263, present in both A 1 AR and A 2A AR, staples the conformation of the EL3.Taking all the above arguments into consideration, from a computational point of view the question arises spontaneously: "could the crystallographic structure of the human A1 AR be a better template than the human A2A AR for the construction of a homology model of the A3 subtype?", the question which provided the scientific motivation in taking up the present research work.To appreciate the qualitative and the quantitative differences obtained by using the two alternative templates, we decided to retrospectively compare the posing and scoring performance of molecular docking for a selection of well-known potent selective agonists and antagonists, summarized in Figures 2 and 3.  Taking all the above arguments into consideration, from a computational point of view the question arises spontaneously: "could the crystallographic structure of the human A 1 AR be a better template than the human A 2A AR for the construction of a homology model of the A 3 subtype?", the question which provided the scientific motivation in taking up the present research work.To appreciate the qualitative and the quantitative differences obtained by using the two alternative templates, we decided to retrospectively compare the posing and scoring performance of molecular docking for a selection of well-known potent selective agonists and antagonists, summarized in Figures 2 and 3. Taking all the above arguments into consideration, from a computational point of view the question arises spontaneously: "could the crystallographic structure of the human A1 AR be a better template than the human A2A AR for the construction of a homology model of the A3 subtype?", the question which provided the scientific motivation in taking up the present research work.To appreciate the qualitative and the quantitative differences obtained by using the two alternative templates, we decided to retrospectively compare the posing and scoring performance of molecular docking for a selection of well-known potent selective agonists and antagonists, summarized in Figures 2 and 3.

Computational Facilities
The MOE suite (Molecular Operating Environment, version 2016.0801)[17] was exploited to perform general molecular modeling operations.All computations were carried out on a 12 CPU (Intel® Xeon® CPU E5-1650 3.80GHz) Linux workstation (distribution 14.04, 64 bit).Docking simulations were performed exploiting the docking program GOLD 5.4.1 [18] and the Scoring Function GoldScore.Plots were generated using Gnuplot 4.6 [19].2D chemical structures were drawn with Marvin Sketch [20].Molecular graphics were performed with the UCSF Chimera package [21].

Protein Preparation
The hA2A and hA1 AR crystallographic structures were retrieved from the RCSB Protein Data Bank (PDB) database [22].In the case of inactive models, the A3-3PWH model was retrieved from the Adenosiland web-platform [23].Fusion proteins and antibody portions were removed, along with ions and any other co-crystallized molecule.Protein structures were ionized using MOE Protonate 3D tool [24], at physiological pH 7.4, temperature 310 K and salt concentration 0.1 M. The Generalized Born (GB) solvation model [25] was used, with an inner dielectric constant value of 2 and the 800 R3 van der Waals (vdW) energy function.MOE energy minimization tool was exploited to minimize contacts between hydrogen atoms, using the Amber12:EHT force field [26], until the RMS of the conjugate gradient was < 0.05 kcal•mol -1 x Å -1 .Sodium ion and its first hydration shell were included during docking simulations on antagonists, basing on one of the highest resolution hA2A AR crystal structures: PDB code 5NM2 [27].The 3D-structures of ligands were built by the MOE-builder tool.Tautomeric states and atoms' hybridization were checked.Geometry minimization was performed by the MMFF94x [28], setting the root mean square gradient < 0.05 kcal•mol -1 •Å -1 .Ligands partial

Computational Facilities
The MOE suite (Molecular Operating Environment, version 2016.0801)[17] was exploited to perform general molecular modeling operations.All computations were carried out on a 12 CPU (Intel®Xeon®CPU E5-1650 3.80GHz) Linux workstation (distribution 14.04, 64 bit).Docking simulations were performed exploiting the docking program GOLD 5.4.1 [18] and the Scoring Function GoldScore.Plots were generated using Gnuplot 4.6 [19].2D chemical structures were drawn with Marvin Sketch [20].Molecular graphics were performed with the UCSF Chimera package [21].

Protein Preparation
The hA 2A and hA 1 AR crystallographic structures were retrieved from the RCSB Protein Data Bank (PDB) database [22].In the case of inactive models, the A 3 -3PWH model was retrieved from the Adenosiland web-platform [23].Fusion proteins and antibody portions were removed, along with ions and any other co-crystallized molecule.Protein structures were ionized using MOE Protonate 3D tool [24], at physiological pH 7.4, temperature 310 K and salt concentration 0.1 M. The Generalized Born (GB) solvation model [25] was used, with an inner dielectric constant value of 2 and the 800 R3 van der Waals (vdW) energy function.MOE energy minimization tool was exploited to minimize contacts between hydrogen atoms, using the Amber12:EHT force field [26], until the RMS of the conjugate gradient was < 0.05 kcal•mol −1 •Å −1 .Sodium ion and its first hydration shell were included during docking simulations on antagonists, basing on one of the highest resolution hA 2A AR crystal structures: PDB code 5NM2 [27].The 3D-structures of ligands were built by the MOE-builder tool.Tautomeric states and atoms' hybridization were checked.Geometry minimization was performed by the MMFF94x [28], setting the root mean square gradient < 0.05 kcal•mol −1 •Å −1 .Ligands partial charges were calculated by means of the PM3/ESP semiempirical Hamiltonian [29].Protein atomic partial charges were computed according to the Amber12:EHT force field.

Sequence Alignment
The hA 1 and hA 2A AR crystallographic structures (5UEN and 3PWH, respectively) were aligned and superposed by MOE Sequence Editor.Structural domains were annotated basing on visual inspection of each TM region.The canonical fasta sequence of the hA 3 adenosine receptor (ID P0DMS8-1) was retrieved from the UniProt database [30] and aligned on the sequences of hA 1 and hA 2A AR (see Table 1).Sequence identity and similarity percentages were calculated for all domains using BLOSUM62 matrix (see Results section).To perform a more robust comparative analysis, we calculated the deviation values between either sequence identity (SI) or sequence similarity (SS) of each domain of the hA 3 AR calculated on the hA 1 (A3/A1) from the one calculated on the hA 2A AR (A3/A2a), as reported in Equations ( 1) and ( 2): Negative values are in favor of the hA 2A subtype, while positive values are of hA 1 .

Homology Modeling
Homology models of both the active and the inactive state were built using the MOE Homology Model Tool.Only aligned residues were overridden.C/N-terminus outgap modeling was disabled, while automatic detection of disulfide bonds was enabled.The maximum number of main chain models was set to 10 and side chain sampling was set to 1 at 300 K. Differently from the intermediates, the final state was ionized at physiological pH (7.4 units) and subjected to medium minimization in order to moderately relieve the steric strain, according to the AMBER12:EHT force field.The final output model was checked, inspecting the wellness of its assigned biophysical parameters and scores, along with good superposition with the reference structure, transmembrane domain side chain rotamers and exposure, dissociation state of charged residues and tautomeric state of the histidine residues.For inactive state model of the hA 3 AR, the following template PDB structures were adopted and then subjected to the comparison (see Section 3): 3PWH, 5NM2 (A 2a ) and 5UEN (A 1 ).For the active state model 2YDO (A 2A -based) and 6D9H (A 1 -based) PDB templates were used.All generated homology models passed with all major biophysical checkpoints, basing on the MOE tool.For the adenosine bound hA 1 AR crystal structure (PDB code: 6D9H), two homology models were carried out by disabling or enabling the Induced Fit option, respectively.This option permits the side chains of the amino acids defining the ligand binding pocket to adapt their conformational position around the ligand during the construction of the homology model.

Molecular Docking
GOLD (Genetic Optimization for Ligand Docking, version 5.4.1)docking algorithm was selected as a conformational search program and Goldscore as a scoring function, referring to previous docking benchmark studies performed by our research group [49][50][51].For any ligand, 20 docking simulation runs were performed on each receptor subtype, searching on a sphere of 20Å radius, centered on the sidechain nitrogen of the conserved N 6.55 .The RMSD threshold for the conformational search was set to 1.0 Å, and 20 poses were collected.

Docking Analysis
Ligand and protein partial charges were calculated using the PM3 method and AMBER12EHT force field respectively.Then, overall electrostatic (Ele) and van der Waals (vdW) energy contributions to the binding energy were calculated by MOE, together with per residue electrostatic and hydrophobic interactions, giving the so-called "Interaction Energy Fingerprints" (IEF) [49].The interactions of most relevant residues were reported in histograms, whose height is proportional to the strength of the interaction.All the plots were produced using GNUPLOT 4.6 [19].

Docking-based Homology Models' Comparison
To appreciate the qualitative and the quantitative differences derived by using the two alternative templates, we decided to retrospectively compare the performance of molecular docking in posing and scoring for a selection of known potent and selective either agonists or antagonists.For each ligand and each homology model, 20 docking poses were retained and the representative one has been selected using the combination of two conditions: (1) the best match in reproducing the two crucial and conserved ligand-receptor interactions (hydrogen-bond with N 6.55 and π-π stacking with F EL2 ); (2) for all poses in agreement with the previous condition, the best combination of the force field-driven energetic Ele and vdW contributions.
The EM correlation coefficient (whose value ranges between -2 and 2) was recently implemented as a statistical tool to evaluate the correlation between two variables dependent on one or more of interest, and to estimate the commutative influence of the latter without statistical bias [52].It quantifies simultaneously (Equation (3), Figure 4): if two variables (numerical results), x and y, have the same absolute value and the same sign (positive or negative); the sign of each variable (sgnx, sgny); the ratio between x and y (|x/y|).Here, we demonstrate that the EM coefficient can be solved as a function of the ratio x/y: In the latter case of white spots, the exponential form of the EM coefficient is expressly related to the ratio of x and y, with a view to individuate the most favored pose.In this way, we determined the models' influence on the binding mode of ligands and, finally, their wellness.one pose is better than the other (x model pose is better than y model pose, blue spots, −+; y model pose is better than x model pose, red spots, +−); 2.
In the latter case of white spots, the exponential form of the EM coefficient is expressly related to the ratio of x and y, with a view to individuate the most favored pose.In this way, we determined the models' influence on the binding mode of ligands and, finally, their wellness.

MMSDocking Video Maker
The in-house MMsDocking video maker tool was exploited to produce videos showing the docking poses, per residue IEFhyd and IEFele data for selected residues.Representations of docking poses were produced using CHIMERA52 [21], two-dimensional depictions were constructed through the cheminformatics toolkit RDKit [52], the heat maps were obtained through GNUPLOT 4.6 [19]; in the end, videos were mounted using MEncoder [53].

Sequence Comparison and Homology Modeling
The bioinformatics analysis of the protein sequences of the four adenosine receptor subtypes has already been presented and discussed in the literature [54][55][56].The overall sequence identity/similarity between human adenosine receptors is relatively high.Interestingly, the adenosine A 3 receptor reveals higher sequence identity/similarity to the adenosine A 1 receptor (48/69%) than to the adenosine A 2A (37/48%) receptor, as summarized in Tables 1 and 2. In homology modeling, the simplest template selection rule is to select the structure with the highest sequence similarity to the modeled sequence [57].Moreover, when it is possible, a template bound to the same or similar ligands as the modeled sequence should generally be used.By the analysis of data shown in Table 2 and Figure 5, i.e., the SI% and SS% for all the receptor domains involved in the recognition of agonists and antagonists, the hA 1 AR can be considered a better possible template than the hA 2A AR.As anticipated in the Introduction, the A3 AR has not yet been crystallized.To appreciate the qualitative and the quantitative differences derived by using hA1 AR as alternative templates of the conventional hA2A AR, we decided to retrospectively compare the performance of molecular docking in posing and scoring a selection of known potent and selective agonists and antagonists, summarized in Figures 2 and 3.Moreover, we have recently reported that the presence of the sodium ion and its first sphere of hydration improves the accuracy and precision of positioning of the antagonist, with respect to the crystallographic poses.This occurs through most of the docking programs and regardless of the chemical nature of the ligand [51].These data are strongly consistent with the crystallographic evidence which would seem to indicate sodium and his first hydration sphere being an integral part of the orthosteric site that recognizes the structure of antagonist, even for those crystallographic structures in which these have not been solved.By logical extension, sodium and its first hydration sphere should be incorporated into all docking studies, taking into account the interaction between an antagonist and its orthosteric receptor cavity.On the contrary, the presence of sodium and its first hydration sphere, as observed crystallographically, are not an integral part of the orthosteric site that recognizes the structure of the agonist.Taking into account the evidence described above, in this work we have constructed two different hA3 AR homology models: one used for the recognition of agonists (in which the sodium cation and its first hydration sphere are absent), and another used for the recognition of antagonists (in which the sodium cation and its first hydration sphere are an integral part of the orthostatic site).

"Agonist-driven" hA3 AR Models
The first comparative validation was performed using the structure of adenosine, its endogenous agonist.We have compared the adenosine binding motif in both hA1 AR and hA2A AR crystallographic structures with those obtained through our docking procedure on the hA3 AR models derived from the use of the two different crystallographic templates, hA1 AR and hA2A AR respectively.As shown in Figure 6, the best docking poses of adenosine in both the hA3 ARs are geometrically similar to the crystallographic poses in hA1 AR and hA2A AR.In particular, in each model, adenosine engages its canonical stabilizing bidentate H-bond of N7 and N 6 H with Asn250 (6.55) and the π-π stacking interaction with Phe168 (EL2).As anticipated in the Introduction, the A 3 AR has not yet been crystallized.To appreciate the qualitative and the quantitative differences derived by using hA 1 AR as alternative templates of the conventional hA 2A AR, we decided to retrospectively compare the performance of molecular docking in posing and scoring a selection of known potent and selective agonists and antagonists, summarized in Figures 2 and 3.Moreover, we have recently reported that the presence of the sodium ion and its first sphere of hydration improves the accuracy and precision of positioning of the antagonist, with respect to the crystallographic poses.This occurs through most of the docking programs and regardless of the chemical nature of the ligand [51].These data are strongly consistent with the crystallographic evidence which would seem to indicate sodium and his first hydration sphere being an integral part of the orthosteric site that recognizes the structure of antagonist, even for those crystallographic structures in which these have not been solved.By logical extension, sodium and its first hydration sphere should be incorporated into all docking studies, taking into account the interaction between an antagonist and its orthosteric receptor cavity.On the contrary, the presence of sodium and its first hydration sphere, as observed crystallographically, are not an integral part of the orthosteric site that recognizes the structure of the agonist.Taking into account the evidence described above, in this work we have constructed two different hA 3 AR homology models: one used for the recognition of agonists (in which the sodium cation and its first hydration sphere are absent), and another used for the recognition of antagonists (in which the sodium cation and its first hydration sphere are an integral part of the orthostatic site).

"Agonist-driven" hA 3 AR Models
The first comparative validation was performed using the structure of adenosine, its endogenous agonist.We have compared the adenosine binding motif in both hA 1 AR and hA 2A AR crystallographic structures with those obtained through our docking procedure on the hA 3 AR models derived from the use of the two different crystallographic templates, hA 1 AR and hA 2A AR respectively.As shown in Figure 6, the best docking poses of adenosine in both the hA 3 ARs are geometrically similar to the crystallographic poses in hA 1 AR and hA 2A AR.In particular, in each model, adenosine engages its canonical stabilizing bidentate H-bond of N7 and N 6 H with Asn250 (6.55) and the π-π stacking interaction with Phe168 (EL2).However, in the case of the hA2A AR model, the ribose ring is not anchored to Thr94 (3.36) by hydrogen bond, conversely, in the hA1 AR model, it is, by a pretty network.This is also confirmed by corresponding energy histograms, in which the hA2A AR model is shown to be highly disfavored not only because of Thr94 (3.36) but also of critical residues such as Val169 (EL2).Moreover, the π-π stacking interaction with Phe168 (EL2) is not verified, as no hydrophobic interaction is detected.It is evident that an increased stabilization of adenosine occurred in the hA1 AR-based model due to some crucial amino acids, i.e.Thr94 (3.36), Phe168 (EL2) and the not conserved residues Val169 (EL2) and Leu90 (3.32).
As for adenosine, all other agonists have been docked to any of the hA3 ARs, and the best docking poses have been summarized in Video S1 (available as Supplementary Information).
Summarizing the docking results, all analyzed agonists were preferably fitted in the orthosteric binding site of the hA1-based models of hA3, as graphically shown in the correlation matrix of Figure 7, where the energetic contributions (electrostatic and van der Waals) assigned to the hA1-based model are compared with the hA2A-based one.The chromatic scale must be interpreted as follows: a) the white color indicates favorable energetic contributions of the specific ligand (reported in the yaxis) for both hA1-based model and hA2A-based models (as for example calculated for ADO and MRS5127); b) the orange color indicates unfavorable energetic contributions of the specific ligand for both hA1-based model and hA2A-based models; c) the blue color indicates favorable energetic contributions of the specific ligand (reported in the y-axis) for the hA1-based model and unfavorable for the hA2A-based model (as for example calculated for MRS7110 and MRS5644); and, finally, d) the red color indicates unfavorable energetic contributions of the specific ligand for the hA1-based model and favorable for the hA2A-based model.Moreover, the numerical value inside each cell describes, in logarithmic form, the ratio of the specific energy contribution calculated between the hA1-based model and the hA2A-based reference model (see Equation ( 2), x and y respectively).From an interpretative point of view, if this value is more negative it means that this energetic contribution is much more favorable for the binding of the ligand versus the hA1-based model compared to the hA2Abased model, whereas when this value is more positive, the energetic contribution is much more However, in the case of the hA 2A AR model, the ribose ring is not anchored to Thr94 (3.36) by hydrogen bond, conversely, in the hA 1 AR model, it is, by a pretty network.This is also confirmed by corresponding energy histograms, in which the hA 2A AR model is shown to be highly disfavored not only because of Thr94 (3.36) but also of critical residues such as Val169 (EL2).Moreover, the π-π stacking interaction with Phe168 (EL2) is not verified, as no hydrophobic interaction is detected.It is evident that an increased stabilization of adenosine occurred in the hA 1 AR-based model due to some crucial amino acids, i.e.Thr94 (3.36), Phe168 (EL2) and the not conserved residues Val169 (EL2) and Leu90 (3.32).
As for adenosine, all other agonists have been docked to any of the hA 3 ARs, and the best docking poses have been summarized in Video S1 (available as Supplementary Information).
Summarizing the docking results, all analyzed agonists were preferably fitted in the orthosteric binding site of the hA 1 -based models of hA 3 , as graphically shown in the correlation matrix of Figure 7, where the energetic contributions (electrostatic and van der Waals) assigned to the hA 1 -based model are compared with the hA 2A -based one.The chromatic scale must be interpreted as follows: (a) the white color indicates favorable energetic contributions of the specific ligand (reported in the y-axis) for both  2), x and y respectively).From an interpretative point of view, if this value is more negative it means that this energetic contribution is much more favorable for the binding of the ligand versus the hA 1 -based model compared to the hA 2A -based model, whereas when this value is more positive, the energetic contribution is much more favorable for the binding of the ligand versus the hA 2A -based model and compared to the hA 1 -based model.From the analysis of the correlation matrix is deducible that all homology models, independently from the original template, are able to properly accommodate all selected agonists inside the orthosteric binding site of the hA3 AR.However, in general, the energetic contributions for agonist-receptor interaction (electrostatic and/or van der Waals) indicates a higher stabilization in the hA1-based model compared to a model based on A2A structure, in particular for the agonists MRS7110 and MRS5644.Interestingly, when docked to the hA1-based model, MRS-7110 achieves a very suitable accommodation for the 2-(4-aryl-1,2,3-triazolyl) substituent, in the proximity to both EL2 and TM2 (Figure 8).The 2-aryl-triazolyl derivative shows the conserved binding pattern of agonists and a ribosyl C2′-endo conformation.Conversely, MRS-5644 is partially displaced from the binding site.MRS-5644 presents a huge planar 2-arylethynyl substituent, which is thought to be responsible for selectivity towards the hA3 AR.A favorable as possible accommodation of the substituent (without prohibitive sterical hindrance) must be determinant for hA3 selectivity.None of the built models is able to predict the most suitable pose for this ligand, but the hA1-based model can accommodate its analog MRS-7710, which shares a geometrically similar substituent in C2.This seems to be consistent with the evidence that the different conformation of the ELs and the TMs' shift observed for the hA1 AR, as compared with the hA2A subtype, leads to a different shape of the binding cavity.It is reasonable to assume that peculiar residues of the hA3 AR sequence could further enhance this aspect, making the binding cavity more prone to locate even rigid bulky substituents, especially towards the TM2 and EL2.In order to explore this hypothesis, we built up another structure, based on the already obtained 6D9H-A3 model, exploiting the MOE Induced Fit functionality in presence of the MRS-7110 ligand pose, whose coordinates had been previously superposed to the MRS-5644 pose (indeed, these ligands differ only for the respective functionalities in C2).
Docking results on the 6D9H-in A3 AR model are visibly better: 5 out 5 poses are in line with the canonical binding interaction geometry of AR ligands, as shown in Video S2.Furthermore, interaction with Thr94 (3.36) and Asn250 (6.55) is significantly more stable for each ligand.From the analysis of the correlation matrix is deducible that all homology models, independently from the original template, are able to properly accommodate all selected agonists inside the orthosteric binding site of the hA 3 AR.However, in general, the energetic contributions for agonist-receptor interaction (electrostatic and/or van der Waals) indicates a higher stabilization in the hA 1 -based model compared to a model based on A 2A structure, in particular for the agonists MRS7110 and MRS5644.Interestingly, when docked to the hA 1 -based model, MRS-7110 achieves a very suitable accommodation for the 2-(4-aryl-1,2,3-triazolyl) substituent, in the proximity to both EL2 and TM2 (Figure 8).The 2-aryl-triazolyl derivative shows the conserved binding pattern of agonists and a ribosyl C2 -endo conformation.Conversely, MRS-5644 is partially displaced from the binding site.MRS-5644 presents a huge planar 2-arylethynyl substituent, which is thought to be responsible for selectivity towards the hA 3 AR.A favorable as possible accommodation of the substituent (without prohibitive sterical hindrance) must be determinant for hA 3 selectivity.None of the built models is able to predict the most suitable pose for this ligand, but the hA 1 -based model can accommodate its analog MRS-7710, which shares a geometrically similar substituent in C2.This seems to be consistent with the evidence that the different conformation of the ELs and the TMs' shift observed for the hA 1 AR, as compared with the hA 2A subtype, leads to a different shape of the binding cavity.It is reasonable to assume that peculiar residues of the hA 3 AR sequence could further enhance this aspect, making the binding cavity more prone to locate even rigid bulky substituents, especially towards the TM2 and EL2.In order to explore this hypothesis, we built up another structure, based on the already obtained 6D9H-A 3 model, exploiting the MOE Induced Fit functionality in presence of the MRS-7110 ligand pose, whose coordinates had been previously superposed to the MRS-5644 pose (indeed, these ligands differ only for the respective functionalities in C2).Docking results on the 6D9H-in A 3 AR model are visibly better: 5 out 5 poses are in line with the canonical binding interaction geometry of AR ligands, as shown in Video S2.Furthermore, interaction with Thr94 (3.36) and Asn250 (6.55) is significantly more stable for each ligand.Importantly, as hypothesized, MRS-5644 reaches a reliable orientation (Figure 8), with the 2-aryl-ethynyl substituent directed towards the TM2, on the side of EL2.Indeed, the model, as compared to the template 6D9H, shows a slight outward shift of the TMs 2, 3 and 7 upper regions, while an inward shift of the EL2 is observed (the measured Cα-RMSD is about 0.28 Å).
These results indicate that, ideally, the orthosteric site of the hA 3 AR active state may accommodate the peculiar substituents of selective agonists in the proximity of EL2 and TMs 2-3 through some expansion of the cavity, favored by the conserved and not conserved residues located in these regions, as Leu90 (3.32) [58], Gln167 (EL2) [59] and Thr94 (3.36).Importantly, as hypothesized, MRS-5644 reaches a reliable orientation (Figure 8), with the 2-arylethynyl substituent directed towards the TM2, on the side of EL2.Indeed, the model, as compared to the template 6D9H, shows a slight outward shift of the TMs 2, 3 and 7 upper regions, while an inward shift of the EL2 is observed (the measured Cα-RMSD is about 0.28 Å).
These results indicate that, ideally, the orthosteric site of the hA3 AR active state may accommodate the peculiar substituents of selective agonists in the proximity of EL2 and TMs 2-3 through some expansion of the cavity, favored by the conserved and not conserved residues located in these regions, as Leu90 (3.32) [58], Gln167 (EL2) [59] and Thr94 (3.36).

"Antagonist-driven" hA3 AR Models
In contrast to the construction of the "agonist-driven" hA3 AR models, in the case of "antagonistdriven" hA3 AR, it has been decided to build up three different homology models, with two of these models using as templates the two hA2A AR structures in complex with the antagonist ZM241385 (PDB codes: 3PWH and 5NM2, respectively), and the third one using as a template the hA1 AR structure in complex with the antagonist DU1 (PDB code: 5UEN).As already anticipated, all "antagonist-driven" hA3 AR models have been built up considering the sodium cation and its first hydration sphere as an integral part of the orthostatic site.
All other antagonists, as listed in Figure 3, have been docked in any of the three hA3 ARs and the best docking poses have been summarized in Video S2 (available as Supplementary Information).
As for the comparison of agonists-driven docking results, even for antagonists the correlation matrix, shown in Figure 9, summarizes the energetic contributions (electrostatic and van der Waals) assigned to the hA1-based and the hA2A 5NM2 model compared to the hA2A 3PWH model.

"Antagonist-driven" hA 3 AR Models
In contrast to the construction of the "agonist-driven" hA 3 AR models, in the case of "antagonist-driven" hA 3 AR, it has been decided to build up three different homology models, with two of these models using as templates the two hA 2A AR structures in complex with the antagonist ZM241385 (PDB codes: 3PWH and 5NM2, respectively), and the third one using as a template the hA 1 AR structure in complex with the antagonist DU1 (PDB code: 5UEN).As already anticipated, all "antagonist-driven" hA 3 AR models have been built up considering the sodium cation and its first hydration sphere as an integral part of the orthostatic site.
All other antagonists, as listed in Figure 3, have been docked in any of the three hA 3 ARs and the best docking poses have been summarized in Video S2 (available as Supplementary Information).
As for the comparison of agonists-driven docking results, even for antagonists the correlation matrix, shown in Figure 9, summarizes the energetic contributions (electrostatic and van der Waals) assigned to the hA 1 -based and the hA 2A 5NM2 model compared to the hA 2A 3PWH model.From the analysis of the correlation matrix, in this case, it is deducible that not all homology models equally stabilize the accommodation of the different antagonists in the orthosteric binding site of hA3 AR.In particular, the energetic contributions for antagonist-receptor interaction (electrostatic and/or van der Waals) would seem to show a much higher stabilization in the hA1-based model, compared to models based on both the A2A structures, for several antagonists such as the highly selective 1,4-dihydropyridine derivative MRS-1334 (especially the S-enantiomers) or the isoquinoline VUF-5455.
As an example, the highly selective 1,4-dihydropyridine derivative MRS-1334 exists in two possible enantiomers, (4R) and (4S).It has been already demonstrated that the (S)-enantiomer is the main responsible for selectivity for this class of compounds [47].The (4S)-enantiomer orients the pnitrobenzyl substituent in the same direction of the aryl-ethynyl one, by intramolecular aromatic stacking interaction, conversely, the (4R)-enantiomer does to the top between TMs 1 and 7 (Figure 10).Enantiomers' poses differ mainly in the orientation of the ethyl ester substituent, up (4S) and down (4R), leading to different interactions as well.The (4R)-enantiomer shows very unfavorable interactions of the 2-phenyl moiety with the EL2 residue Val169.Conversely, the (4S)-enantiomer ethyl ester substituent contracts favorable interactions with Val169 and strong hydrophobic interactions with Phe168 (EL2).Hydrophobic interactions are established also with Val65 and Leu90.In general, the latter gives the main contribution to both enantiomers of MRS-1344 poses.These data are largely in favor of the (4S)-enantiomer, consistently with enantioselectivity observed.Asn250 (6.55) contribution is not determinant in this case, suggesting that excluding inevitable limitations of homology modeling, the atypical dihydropyridine scaffold in its S-configuration could be unable to engage hydrogen bonds with Asn250 (6.55), except if we consider a rotation of the amidic side chain, i.e. the driving component of the final state may be mainly hydrophobic.From the analysis of the correlation matrix, in this case, it is deducible that not all homology models equally stabilize the accommodation of the different antagonists in the orthosteric binding site of hA 3 AR.In particular, the energetic contributions for antagonist-receptor interaction (electrostatic and/or van der Waals) would seem to show a much higher stabilization in the hA 1 -based model, compared to models based on both the A 2A structures, for several antagonists such as the highly selective 1,4-dihydropyridine derivative MRS-1334 (especially the S-enantiomers) or the isoquinoline VUF-5455.
As an example, the highly selective 1,4-dihydropyridine derivative MRS-1334 exists in two possible enantiomers, (4R) and (4S).It has been already demonstrated that the (S)-enantiomer is the main responsible for selectivity for this class of compounds [47].The (4S)-enantiomer orients the p-nitrobenzyl substituent in the same direction of the aryl-ethynyl one, by intramolecular aromatic stacking interaction, conversely, the (4R)-enantiomer does to the top between TMs 1 and 7 (Figure 10).Enantiomers' poses differ mainly in the orientation of the ethyl ester substituent, up (4S) and down (4R), leading to different interactions as well.The (4R)-enantiomer shows very unfavorable interactions of the 2-phenyl moiety with the EL2 residue Val169.Conversely, the (4S)-enantiomer ethyl ester substituent contracts favorable interactions with Val169 and strong hydrophobic interactions with Phe168 (EL2).Hydrophobic interactions are established also with Val65 and Leu90.In general, the latter gives the main contribution to both enantiomers of MRS-1344 poses.These data are largely in favor of the (4S)-enantiomer, consistently with enantioselectivity observed.Asn250 (6.55) contribution is not determinant in this case, suggesting that excluding inevitable limitations of homology modeling, the atypical dihydropyridine scaffold in its S-configuration could be unable to engage hydrogen bonds with Asn250 (6.55), except if we consider a rotation of the amidic side chain, i.e. the driving component of the final state may be mainly hydrophobic.The potent antagonist VUF-5455, characterized by a not xanthinic scaffold as well, reveals a peculiar accommodation when docked to the hA1 AR model.As shown in Figure 11 and Video S2, in this case, this antagonist establishes a highly favorable hydrogen bond with the canonical Asn250 (6.55), which is unfavorable instead in other models.Interaction with Leu90 (3.32), Phe168 (EL2) and Thr94 (3.36) is crucial, suggesting that also in this case these residues can play a major role in ligand recognition, creating highly selective hydrophobic locations for selective substituents.The potent antagonist VUF-5455, characterized by a not xanthinic scaffold as well, reveals a peculiar accommodation when docked to the hA 1 AR model.As shown in Figure 11 and Video S2, in this case, this antagonist establishes a highly favorable hydrogen bond with the canonical Asn250 (6.55), which is unfavorable instead in other models.Interaction with Leu90 (3.32), Phe168 (EL2) and Thr94 (3.36) is crucial, suggesting that also in this case these residues can play a major role in ligand recognition, creating highly selective hydrophobic locations for selective substituents.

Conclusions
On the basis of the results obtained from this preliminary study, it is possible to affirm that the human A3 receptor model based on the crystallographic structure of the A1 subtype can represent a valid alternative to the models conventionally used today based on the available A2A structures, as supported by the comparative analysis of agonists and antagonists docking simulations on both models.A new series of molecular dynamics simulations (unbiased and supervised) are ongoing in our laboratory to accurately explore the time-dependent behavior of both models in terms of threedimensional structure stability and in their capability to recognize both agonists and antagonists starting from an unbound state.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1,Video S1: Agonist docking poses and energy fingerprint, Video S2: Antagonist docking poses and energy fingerprint.
Author Contributions: The overall work was driven within a strong joint framework between E.M. and S.M.

Conclusions
On the basis of the results obtained from this preliminary study, it is possible to affirm that the human A 3 receptor model based on the crystallographic structure of the A 1 subtype can represent a valid alternative to the models conventionally used today based on the available A 2A structures, as supported by the comparative analysis of agonists and antagonists docking simulations on both models.A new series of molecular dynamics simulations (unbiased and supervised) are ongoing in our laboratory to accurately explore the time-dependent behavior of both models in terms of three-dimensional structure stability and in their capability to recognize both agonists and antagonists starting from an unbound state.
Author Contributions: The overall work was driven within a strong joint framework between E.M. and S.M.
Funding: This research received no external funding.Acknowledgments: S.M. is very grateful to Chemical Computing Group and Acellera for the scientific and technical partnerships.MMS lab gratefully acknowledges the support of NVIDIA Corporation with the donation of the Titan V GPU used for this research.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Structural superposition of hA1 (dark blue) and hA2a (gold) ARs crystal structures in their active state (PDB code: 6D9H and 2YDO, respectively).Panel A and B representations differ only for camera orientation.Additional A2A disulfide bridges are highlighted in red.

Figure 2 .
Figure 2. Chemical structures of agonists selected for docking studies.Experimental Ki values (in nM or activity assay percentages) for each receptor subtype are also reported for comparison.

Figure 1 .
Figure 1.Structural superposition of hA 1 (dark blue) and hA 2a (gold) ARs crystal structures in their active state (PDB code: 6D9H and 2YDO, respectively).Panel (A) and (B) representations differ only for camera orientation.Additional A 2A disulfide bridges are highlighted in red.

Figure 1 .
Figure 1.Structural superposition of hA1 (dark blue) and hA2a (gold) ARs crystal structures in their active state (PDB code: 6D9H and 2YDO, respectively).Panel A and B representations differ only for camera orientation.Additional A2A disulfide bridges are highlighted in red.

Figure 2 .
Figure 2. Chemical structures of agonists selected for docking studies.Experimental Ki values (in nM or activity assay percentages) for each receptor subtype are also reported for comparison.

Figure 2 .
Figure 2. Chemical structures of agonists selected for docking studies.Experimental Ki values (in nM or activity assay percentages) for each receptor subtype are also reported for comparison.

Figure 3 .
Figure 3.Chemical structures of antagonists selected for docking studies.Experimental Ki values (in nM or activity assay percentages) for each receptor subtype are also reported for comparison.

Figure 3 .
Figure 3.Chemical structures of antagonists selected for docking studies.Experimental Ki values (in nM or activity assay percentages) for each receptor subtype are also reported for comparison.

Figure 4 .
Figure 4. Graphical representation of the EM correlation coefficient (Equation (3)), a Gaussian function of the logarithmic ratio between x and y absolute values (N).

Figure 4 .
Figure 4. Graphical representation of the EM correlation coefficient (Equation (3)), a Gaussian function of the logarithmic ratio between x and y absolute values (N).In this work, y, specifically, represents the energy contribution of the ligand (Ele or vdW) calculated on the reference model complex (A2a-based model), while x represents the energy contribution of the same ligand calculated on the model subjected to the comparison (A1-based model).Positive energy values indicate disfavored binding poses, negative values in contrast indicate favorable binding poses.Four combinations of x and y energy values are possible, giving four different possible behaviors: 1.one pose is better than the other (x model pose is better than y model pose, blue spots, −+; y model pose is better than x model pose, red spots, +−); 2.both the poses are disfavored (++, orange spots); 3.both the poses are favored (−−, white spots).

Figure 5 .
Figure 5. Sequence comparison calculations on each topological position and overall alignment.The percentage reported on the y-axis indicates how much the difference in identity (ΔSI, cyan), or in similarity (ΔSS, orange), favors alignment towards the hA1 or hA2a AR.

Figure 5 .
Figure 5. Sequence comparison calculations on each topological position and overall alignment.The percentage reported on the y-axis indicates how much the difference in identity (∆SI, cyan), or in similarity (∆SS, orange), favors alignment towards the hA1 or hA2a AR.

19 Figure 6 .
Figure 6.Per residue energy histograms and selected docking pose of adenosine on either the hA2A based model (a) or the induced fit hA1-based model (b).

Figure 6 .
Figure 6.Per residue energy histograms and selected docking pose of adenosine on either the hA 2A based model (a) or the induced fit hA 1 -based model (b).
hA 1 -based model and hA 2A -based models (as for example calculated for ADO and MRS5127); (b) the orange color indicates unfavorable energetic contributions of the specific ligand for both hA 1 -based model and hA 2A -based models; (c) the blue color indicates favorable energetic contributions of the specific ligand (reported in the y-axis) for the hA 1 -based model and unfavorable for the hA 2A -based model (as for example calculated for MRS7110 and MRS5644); and, finally, (d) the red color indicates unfavorable energetic contributions of the specific ligand for the hA 1 -based model and favorable for the hA 2A -based model.Moreover, the numerical value inside each cell describes, in logarithmic form, the ratio of the specific energy contribution calculated between the hA 1 -based model and the hA 2A -based reference model (see Equation ( Appl.Sci.2018, 8, x FOR PEER REVIEW 11 of 19 favorable for the binding of the ligand versus the hA2A-based model and compared to the hA1-based model.

Figure 7 .
Figure 7. EM correlation matrix of energy contributions (electrostatic and van der Waals) for agonist poses (y-axis) obtained from the hA1-based model with induce fit (6D9H-in) and without induce fit (6D9H), as compared to the reference model (2YDO).

Figure 7 .
Figure 7. EM correlation matrix of energy contributions (electrostatic and van der Waals) for agonist poses (y-axis) obtained from the hA 1 -based model with induce fit (6D9H-in) and without induce fit (6D9H), as compared to the reference model (2YDO).

Figure 8 .
Figure 8. Per residue energy histograms and selected docking poses of MRS7110 (panel A) and MRS5644 (panel B) on either the hA2A based model (a) or the hA1 based models (b, c).

Figure 8 .
Figure 8. Per residue energy histograms and selected docking poses of MRS7110 (panel A) and MRS5644 (panel B) on either the hA 2A based model (a) or the hA 1 based models (b,c).

19 Figure 9 .
Figure 9. EM correlation matrix of energy contributions (electrostatic and van der Waals) for antagonists poses (y-axis), obtained from the hA2A-based model (5NM2) and the hA1based model (5UEN), in comparison with the reference model (3PWH).

Figure 9 .
Figure 9. EM correlation matrix of energy contributions (electrostatic and van der Waals) for antagonists poses (y-axis), obtained from the hA 2A -based model (5NM2) and the hA 1 -based model (5UEN), in comparison with the reference model (3PWH).

Figure 11 .
Figure 11.Per-residue energy histograms and selected docking poses of VUF-5455 on either the hA2A-based models (A, B) or the hA1-based model (C).

Figure 11 .
Figure 11.Per-residue energy histograms and selected docking poses of VUF-5455 on either the hA2A-based models (A,B) or the hA 1 -based model (C).

Table 1 .
Sequence alignment of the hA3 AR canonical isoform primary sequence (second row) on the PDB crystallographic structures 5UEN (hA 1 AR) and 3PWH (hA 2A AR), basing on respective structural domains.

Table 2 .
Sequence identity (SI) and sequence similarity (SS) percentages calculated for each alignment performed on structural domains of crystal structures 5UEN (hA 1 AR) and 3PWH (hA 2A ).Overall values are also reported along with domains showing the highest SI and SS percentages.