Pharmacological Characterisation of Novel Adenosine Receptor A3R Antagonists

Background and Purpose The adenosine A3 receptor (A3R) belongs to a family of four adenosine receptor (AR) subtypes which all play distinct roles throughout the body. A3R antagonists have been described as potential treatments for numerous diseases including asthma. Given the similarity between ARs orthosteric binding sites, obtaining highly selective receptor antagonists is a challenging but critical task. Experimental approach 39 potential A3R, antagonists were screened using agonist-induced inhibition of cAMP. Positive hits were assessed for AR subtype selectivity through cAMP accumulation assays. The antagonist affinity was determined using Schild analysis (pA2 values) and fluorescent ligand binding. Further, a likely binding pose of the most potent antagonist (K18) was determined through molecular dynamics (MD) simulations and consistent calculated binding free energy differences between K18 and congeners, using a homology model of A3R, combined with mutagenesis studies. Key Results We demonstrate that K18, which contains a 3-(dichlorophenyl)-isoxazole group connected through carbonyloxycarboximidamide fragment with a 1,3-thiazole ring, is a specific A3R (<1 µM) competitive antagonist. Structure-activity relationship investigations revealed that loss of the 3-(dichlorophenyl)-isoxazole group significantly attenuated K18 antagonistic potency. Mutagenic studies supported by MD simulations identified the residues important for binding in the A3R orthosteric site. Finally, we introduce a model that enables estimates of the equilibrium binding affinity for rapidly disassociating compounds from real-time fluorescent ligand-binding studies. Conclusions and Implications These results demonstrate the pharmacological characterisation of a selective competitive A3R antagonist and the description of its orthosteric binding mode.


INTRODUCTION
The adenosine A3 receptor (A3R), belongs to a family of four adenosine receptor (AR) subtypes (A1R, A2AR, A2BR and A3R), and is involved in a range of pathologies including cardiovascular, neurological and tumour-related diseases. In particular, mast cell regulation and myocardial preconditioning are key physiological processes regulated by the A3R (Fredholm et al., 2011). Unsurprisingly therefore, A3R is a pharmaceutical target. Interestingly, the A3R has been described as enigmatic, whereby many of the effects attributed to A3Rs are contradictory (Gessi et al., 2008). Despite this, A3R antagonists having been described as potential treatments of asthma, chronic obstructive pulmonary disease (COPD) and glaucoma (Miwatashi et al., 2008, Okamura et al., 2004, Haeusler et al., 2015, to name a few, and continuous research into both agonists and antagonists at the A3R are warranted. One of the challenges associated with the druggability of the AR family has been the targeting of individual subtypes with sufficient specificity to limit off-target side effects (Chen et al., 2013). In silico screening of vast compound libraries against receptor structures, known as structuralbased drug design, offers huge potential in the development of highly selective ligands.
Although all AR members are activated by the endogenous agonist adenosine, the A2AR and A2BR are predominantly Gs-coupled whereas A1R and A3R generally couple to Gi/o. This classical pathway following A3R activation and Gi/o coupling is the inhibition of adenylate cyclase (AC) resulting in a decrease in cAMP levels. Extracellular signal-regulated kinase 1/2 (ERK1/2) activation has also been described as downstream of A3R and is reported to be dependent on βγ-subunits released from pertussis toxin (PTX)-sensitive Gi/o proteins, phosphatidylinostitol-3-kinase (PI3K), the small GTP binding protein Ras, and MAP/ERK kinase (MEK) (Schulte and Fredholm, 2002). In addition to Gi/o, A3R has also been reported to couple to Gq, leading to phospholipase C (PLC) activation and ultimately elevation of intracellular inositol 1,4,5-trisphosphate (IP3) and calcium (Ca 2+ ) levels (Gessi et al., 2008).
The A2AR is one of the best structurally characterised G protein-coupled receptors (GPCRs), with multiple crystal structures (both active and inactive) available including that bound to an engineered G protein (Carpenter et al., 2016) and the A2AR bound to the agonists (5′-(Nethylcarboxamido)adenosine) (NECA) and adenosine (Lebon et al., 2011), CGS 21689 (Lebon et al., 2015, UK-432097 (Xu et al., 2011) and the antagonists ZM241385 (Liu et al., 2012, Doré et al., 2013, Jaakola et al., 2008, PSB36, caffeine and theophylline (Cheng et al., 2017). Although the remaining AR subtypes have proven more difficult to crystallise with the A3R structure yet to be resolved, there are a number of A1R structures published including the adenosine-bound A1R-Gi complex (Draper-Joyce et al., 2018) and antagonist bound structures; DU172 (Glukhova et al., 2017) and PSB36 (Cheng et al., 2017). Thus, structuralbased drug design offers huge potential in the development of highly selective ligands (Carlsson et al., 2010, Katritch et al., 2010, Lenselink et al., 2016, Lagarias et al., 2018. The limited availability of diverse high-resolution structures of the remaining AR subtypes bound to pharmacologically distinct ligands has meant there is a discrepancy between the capability to predict compound binding versus pharmacological behaviour; partial agonism, inverse agonism, biased agonist, antagonism, allosteric modulation etc (Sexton and Christopoulos, 2018). With this in mind, the potential antagonists (K1-K25, K28 and K35) identified in our previously published virtual screening investigation and binding experiments (Lagarias et al., 2018) and some newly identified potential antagonists (K26, K27, K29-K34 and K36-K39) were pharmacologically characterised using cAMP accumulation and ERK1/2 phorphorylation assays. We were able to identify a potent and selective A3R antagonist, K18 (O4-{[3-(2,6dichlorophenyl)-5-methylisoxazol-4-yl]carbonyl}-2-methyl-1,3-thiazole-4carbohydroximamide) and, using molecular dynamic (MD) simulations combined with sitedirected mutagenesis, elude to the potential binding site. Binding free energy calculations of similar in structure analogs of K18 were consistent with the proposed A3R orthosteric binding area. Kinetic binding experiments of K5, K17 and K18 using a bioluminescence resonance energy transfer (BRET) method combined with functional assays led to the identification of important structural features of K18 for binding and activity. Further evaluation of this compound (and structurally related analogues) may afford a novel therapeutic benefit in pathologies such as inflammation and asthma.

Identification of A3R selective antagonists
We set out to conduct a functional screen of 39 compounds for the identification of A3R antagonists, some of which have previously been identified to bind one of the three AR subtypes; A1R, A3R or A2AR (Lagarias et al., 2018). Our screen was conducted using A3R expressing Flp-In™-Chinese hamster ovary (CHO) cells where cAMP accumulation was detected following a combined stimulation of 1 µM forskolin (to allow A3R mediated Gi/o response to be observed), 1 µM tested compound and the predetermined IC80 concentration of NECA (3.16 nM). This initial screen was blinded with each compound numbered without the corresponding name or chemical structure (K1-39). Compound K1-39 were identified by unblinding (Table 1 and Supplementary Table 1) but are hereinafter referred to as their denoted 'K' number. For the purpose of structure-activity relationships studies, the new compounds (K26, K27, K29-K34 and K36-K39), were tested through this functional assay and radioligand binding (Supplementary Table 1). As expected, co-stimulation with 10 µM of both forskolin and NECA reduced the cAMP accumulation when compared to 10 µM forskolin alone and this was reversed with the known A3R antagonist MRS 1220 (Table 1 and Supplementary  Fig 1). Compounds K1, K10, K11, K17, K18, K20, K23, K25 and K32 were identified as potential antagonists at the A3R through their ability to elevate cAMP accumulation when compared to forskolin and NECA co-stimulation. Of the nine potential A3R antagonists, eight (K1, K10, K17, K18, K20, K23, K25 and K32) were confirmed as antagonists at the tested concentration of 1 µM (Supplementary Fig. 2 and Supplementary Table 2). K8, despite showing no binding at any AR subtype (Lagarias et al., 2018), showed a reduced cAMP accumulation. We tested K8 for agonist activity at the A3R but was found to be no different to DMSO ( Supplementary Fig. 3).
A number of compounds previously documented (K5, K9, K21, K22 and K24;Lagarias et al., 2018) or determined in this study (K26, K27 and K34) to have micromolar binding affinity for A3R showed no activity in our functional screen (Table 1, Supplementary Table 1). These compounds, with a Ki in the low micromolar range, were further tested to ensure our functional screen was robust. In addition, compound K11 with a previously determined low micromolar Ki and a similar structure to the active K10 and K32 was also tested at the higher concentration. Full inhibition curves of NECA in the presence or absence of tested compounds (1 µM or 10 µM) were determined in A3R Flp-In CHO cells ( Supplementary Fig. 4, Supplementary Table. 3). All nine compounds (K5, K9, K11, K21, K22, K24, K26, K27 and  K34) reduced the NECA potency at the highest tested concentration (10 µM) but showed no effect at 1 µM and thus appear to be low potency antagonists at the A3R.

AR subtype selectivity and specificity
Stimulation of A3R Flp-In CHO or CHO-K1 cells expressing one of the remaining AR subtypes (A1R, A2AR or A2BR) with a single high concentration of antagonist (10 μM) and increasing concentrations of NECA identified K10, K17, K18 and K25 as A3R selective antagonists, with no apparent antagonism at the remaining AR subtypes (Fig. 1). K20 and K23 were antagonists at both the A1R and A3R ( Fig. 1 and Table 2). K1, K20 and K23 showed weak antagonism of the A2AR and K32 was the only tested antagonist which showed any A2BR activity. These selectivity findings agree with our previously published radioligand binding data (Lagarias et al., 2018) and are summarised in Table 2.
Characterisation of competitive antagonists at the A3R All eight A3R antagonists were confirmed to antagonise IB-MECA agonism ( Fig. 2 and Table  3) and NECA agonism (Supplementary Fig. 5 and Supplementrary Table 4) in a concentrationdependent manner. Schild analysis of the antagonism of both NECA or IB-MECA stimulated cAMP inhibition characterised K10, K17, K18, K20, K23 and K32 as competitive antagonists at the A3R with a slope not significantly different from unity (Supplementary Fig. 5 and Fig. 2). Interestingly, the slope deviated from unity for K1 (in experiments looking at competition with NECA but not IB-MECA) and K25 suggesting a more complicated mechanism of antagonism at the A3R is in play. K20 and K23 were also characterised as competitive antagonists at the A1R with a Schild slope not significantly different from unity (Supplementary Fig. 6 and Supplementary Table 5).
When comparing the activity of A3R selective antagonists (K10, K17, K18 and K25), K18 was the most potent, showed A3R specificity and greater A3R binding affinity (Table 2). It should be noted however, that the original competition binding experiments that identified the panel of antagonist was performed using [ 3 H]HEMADO (Lagarias et al., 2018). To ensure that the different ligand used in our studies was not influencing our characterisation of the compounds we assessed the ability of K18 to antagonise cAMP inhibition by HEMADO at the A3R and compared its potency to K17(Supplementary Fig. 7 and Table 6). As we observed for both NECA and IB-MECA, K18 remained the most potent antagonist at the A3R and we propose it as our lead compound.
We wanted to determine if our lead A3R antagonist, K18, could also reduce the potency of IB-MECA when an alternative downstream signalling component was measured; ERK1/2 phosphorylation (Fig. 3). In line with previously reported findings (Schulte and Fredholm, 2002), agonists at the A3R caused an increase in ERK1/2 phosphorylation after 5 minutes, with IB-MECA 10-fold more potent than NECA ( Supplementary Fig. 8). As previously reported (Graham et al., 2001, Schulte andFredholm, 2002), this was entirely Gi/o-mediated, as demonstrated by the abolished pERK1/2 level in PTX treated A3R Flp-In™-CHO stimulated with NECA/IB-MECA ( Supplementary Fig. 8). The pERK1/2 level following Phorbol 12myristate 13-acetate (PMA) stimulation was entirely unaffected by PTX treatment (Supplementary Fig. 8). Perhaps unsurprisingly, K18 was able to antagonise A3R-mediated evaluate the energetic contributions for their binding (Table 4). The calculated ranking in the binding free energies were in agreement with experimental differences in potencies. Binding free energies (ΔGeff) is calculated as the difference in energetic components between the complex, the apoprotein and the ligand. These components include the difference in electrostatic energy of binding interactions (Eelec), the difference in the van der Waals energy of binding interactions (EvdW) and the difference in the solvation energy (ΔGsolv) ( Table 4). The calculations suggested that the major difference between the energetic components of ΔGeff values for K5, K17 and K18 is on the solvation energies. The two chlorine atoms make K18 more lipophilic and reduce the energy required to transfer the compound from solution to the binding area, increasing the free energy of binding and activity compared to K17 and K5. Interestingly, following MD simulations of the unpublised compounds (K26, K27, K29-K34 and K36-K39) we observed that compounds K26, with a o-diphenylcarbonyl, had low micromolar A3R binding affinities (Supplementary Table 1) and according to the MD simulations of K26 in complex with A3R ( Supplementary Fig. 11) had a similar binding pose to that of K18 (Fig. 4). However, in our functional assays, K26 (and K34, which also had a o-diphenylcarbonyl and low micromolar binding) showed weak antagonistic potency below the concentration of 1 µM ( Supplementary Fig. 3) suggesting a more complex binding mode is present. We observed that the p-substitution in compounds K29 and K36-38 was not favorable for binding at all since this led to a loss of the van der Waals interaction with the hydrophobic area of the A3R towards TM5 and TM6; as was demonstrated in MD simulations for K36 ( Supplementary Fig. 11).
Finally, we also examined how the activity was affected when the 4-thiazolyl in the mid-range antagonist K17 was changed to 2-,3-or 4-pyridyl in compounds K32, K10, K11 which bind to A3R (Table 1). We found antagonistic activity only for compounds K32 and K10; compared to K11, in compounds K32 and K10 the pyridine nitrogen can interact with N250 6.55 due to their proximate positions in binding conformation (see Fig. 4B for K17). This interaction appears to be preserved with both the 2-,3-pyridyl groups but lost when the nitrogen is in the 4-position. Thus, while we have been able to identify compounds that are structurally able to mimic some of the features of compound K18 (and its derivaties K17 and K5), K18 remains the most potent antagonist present within this study.

Experimental evaluation of the binding mode of K18 at A3R
The potential binding site of our lead A3R selective antagonist, K18, was investigated through the use of point mutations as an experimental approach to give insight into structure-function relationships. The determination of critical residues for antagonist binding becomes particularly difficult in the case of competitive antagonists whereby important amino acids are likely overlapping with those for agonist binding. Through performing Schild analysis, whereby the pA2 is independent of agonist, we were able to experimentally determine the effect of receptor mutation on antagonist binding. Whereas an increase in the pA2 for a particular mutant when compared to WT suggested the antagonist was more potent, a decrease indicated a reduced potency. Of the identified residues predicted to mediate an interaction between K18 and the A3R, the ones which showed (according to the MD simulations) the most frequent and the most important contacts were chosen for investigation and included amino acids L90 3.32 , F168 5.29 , V169 5.30 , M177 5.40 , L246 6.51 , I249 6.54 , N250 6.55 , L264 7.34 and I268 7.39 (Fig.  4). Site-directed mutagenesis was performed replacing each residue with an alanine and expressed then in the Flp-In-CHO™ cells lines. Each mutant was then screened for their ability to supress forskolin-induced cAMP accumulation in response to NECA/IB MECA stimulation in the presence and absence of K18.
Mutation of residues F168 5.29 , L246 6.51 , N250 6.55 and I268 7.39 abolished agonist induced suppression of forskolin-induced cAMP accumulation and were discontinued in this study (Stamatis et al., 2019, in preparation). Both L90A 3.32 and M177A 5.40 showed a significantly decreased NECA and IB-MECA potency. L264A 7.34 showed a slight decrease in IB-MECA potency whereas the potency of NECA was similar to WT. Whereas the NECA stimulated cAMP inhibition in V169A 5.30 or I249A 6.54 expressing Flp-In CHOs was comparable to WT, the IB-MECA stimulated cAMP inhibition was enhanced in potency ( Table 5). Mutation of V169 5.30 to glutamate, the amino acid present in the remaining three AR subtypes, enhanced both NECA and IB-MECA potency.

Schild analysis of K18 at WT and mutant A3R
The pA2 values obtained through conducting Schild analysis of K18 at WT and mutant A3R were compared in order to determine the potential antagonist binding site (Fig. 5, Table 5). The pA2 value for I249A 6.54 A3R was similar to WT, whereas M177A 5.40 and V169A 5.30 were significantly smaller. Interestingly we found an increase in the pA2 for L90A 3.32 and L264A 7.34 when compared to WT, suggesting an enhanced ability of K18 to act as an antagonist. Our confidence in the obtained pA2 values for K18 was enhanced by testing with NECA and IB-MECA at an A3R mutant which caused enhanced activity (L90A 3.32 ). As would be expected, the pA2 values for this mutant was not significantly different between agonists, confirming agonist independence ( Supplementary Fig. 12). These experimental findings are reflected in our final binding pose of K18 at the WT A3R (Fig. 4).

Kinetics of A3R antagonists determined through BRET
BRET techniques have been successfully used to determine the real time kinetics of ligand binding to GPCRs (Stoddart et al., 2018). In BRET ligand-binding experiments, we investigated the ability of the selective A3R antagonist MRS 1220, K5, K17 or K18 to inhibit specific binding of the fluorescent A3R antagonist CA200645 to Nluc-A3R. The kinetic parameters for CA200645 at Nluc-A3R were initially determined as Kon (k1) = 2.86 ± 0.89 x 10 7 M -1 , Koff (k2) = 0.4397 ± 0.014 min -1 with a KD of 17.92 ± 4.45 nM. (Supplementary Fig 13). Our MRS 1220 kinetic data was fit with the original 'kinetic of competitive binding' model (Motulsky and Mahan, 1984; built into GraphPad Prism 8.0) with a determined Kon (k3) and Koff (k4) rate of 3.25 ± 0.28 x 10 8 M -1 min -1 and 0.0248 ± 0.005 min -1 , respectively. This gave a residence time (RT) (RT = 1/Koff) of 40.32 min. It was noticed in the analysis for K5, K17 and K18 that the fit in some cases was ambiguous (Regression with Prism 8: "Ambiguous", 2019) and/or the fitted value of the compound dissociation rate constant was high (k4 > 1 min -1 , corresponding to a dissociation t1/2 of < 42 sec). In order to determine the reliability of the fitted k4 value, data were also analysed using an equation that assumes compound dissociation is too rapid for the dissociation rate constant to be determined reliably and the fits to the two equations compared ("Kinetics of competitive binding, rapid competitor dissociation", derived in the Appendix I). This model allowed estimate of the equilibrium binding affinity of the compound (Ki) but not the binding kinetics of K5, K17 and K18 (Supplementary Fig. 14 and Table 4). These pKi values were found to be similar to those calculated through fitting the Cheng-Prusoff equation (Cheng and Prusoff, 1973) and notably, the order of affinity for K5, K17 and K18 reflected that determined through Schild analysis and previously published radioligand binding (Table 4).

DISCUSSION
In silico structure-based drug design efforts in ligand discovery, using molecular docking calculations, have proven to be highly successful (Meng et al., 2011). The broad and similar orthosteric binding site of ARs makes the determination of AR subtype selective compound a challenging task. Indeed, given the similarity between ARs orthosteric binding sites, the search for an AR subtype specific compound often leads to compounds active at more than one of the AR subtypes (Kolb et al., 2012). Given that AR subtypes play distinct roles throughout the body, obtaining highly specific receptor antagonists and agonists is crucial. Here, we presented the pharmacological characterisation of eight A3R antagonists determined though virtual screening. Of these eight compounds, K10, K17, K18, K20, K23 and K32 were determined to be competitive. Whereas K20 and K23 were antagonists at both the A1R and A3R, K10, K17, K18 and K25 were A3R selective antagonists. Indeed, we found no functional activity, or indeed binding affinity (< 30 µM), at the other AR subtypes.
K1, K20 and K23 showed weak antagonism of the A2AR and K32 was the only tested antagonist which showed A2BR antagonisitc potency. These selectivity findings were in agreement with our radioligand binding data (presented here and in Lagarias et al., 2018 for K1-25, K28 and K35). However, a number of compounds previously determined to have micromolar binding affinity for A3R (K5, K9, K21, K22, K24, K26, K27 and K34), showed no antagonistic potency in our initial functional screen. Further testing confirmed that these compounds were low potency antagonists and, although supporting the previously published radioligand binding data, confirmed the need for functional testing: not all compounds with binding affinity showed high functional potency.
We showed the A3R, when expressed in Flp-In TM -CHO cells, displays constitutive activity. Compounds which preferably bind to the inactive (R) state, decreasing the level of constitutive activity (Giraldo et al., 2007) and in the case of a Gi/o -coupled GPCR leading to an elevated cAMP, are referred to as inverse agonists. All eight characterised A3R antagonists and both characterised A1R antagonists (K20 and K23) were found to act as inverse agonists. We also reported an elevation in cAMP accumulation when cells were stimulated with DMSO, which was concentration-dependent. Given that even low concentrations of DMSO has been reported to interfere with important cellular processes (Tunçer et al., 2018), the interpretation of this data should be made with caution. The initial virtual screening described in Lagarias et al., 2018 was carried out using a combination of a ligand-based and structure-based strategy and the experimental structure of A2AR in complex with the antagonist ZM241385 (PDB ID 3EML) (Jaakola et al., 2008), described as A2AR selective antagonist and inverse agonist (Lebon et al., 2011). Our high hit rate for A3R selective antagonist appears counter-intuitive since the ligand-based virtual screening tool Rapid Overlay of Chemical Structures (ROCS) was used to predict structures similar to ZM241385 (Lagarias et al., 2018). Indeed, ZM241385 has little affinity for A3R and 500-to 1000-fold selectivity for A2AR over A1R. However, as has been previously reported, the search for an AR subtype specific compound often leads to compounds active at multiple AR subtypes (Kolb et al., 2012), likely due to their similar binding site.
We hypothesize that the presence of a chloro substituent in the phenyl ring of 3-phenylisoxazole favors A3R affinity and activity, as following 0Cl < 1Cl < 2Cl i.e. K5 < K17 < K18. This theory is supported by both our radioligand binding, NanoBRET ligand-binding and functional data which determine the relative potency and affinity of the three related compounds K5, K17 and K18 as K5 < K17 < K18. The MD simulations showed that these compounds adopted a similar binding mode at the A3R orthosteric binding site, but the freeenergy calculations showed that the two chlorine atoms in K18 increases its lipophilicity, thus allowing it to more efficiently leave the solution state and enter the highly lipophilic binding area.
For the first time, we demonstrate the utilisation of a new model which expands on the 'Kinetic of competitive binding' model (Motulsky and Mahan, 1984; built into Prism) for fitting fast kinetics data obtained from NanoBRET experiments and assumes the unlabelled ligand rapidly equilibrates with the free receptor. Very rapid competitor dissociation can lead to failure of the fit, eliciting either an ambiguous fit (Regression with Prism 8: "Ambiguous", 2019) or unrealistically large K3 and K4 values. Whereas we were able to successfully fit the MRS 1220 kinetic data with the Motulsky and Mahan model due to its slow dissociation, fitting of K5, K17 and K18 kinetic data with this model often resulted in an ambiguous fit. Our new model, assuming fast compound dissociation, successfully fit the data and allowed the determination of binding affinity. In the cases where the data was able to fit the Motulsky and Mahan model, the dissociation constant was higher (of the order of 1 min-1), indicating rapid dissociation. Although we found nearly a 10-fold differences in determined binding affinity for MRS 1220, K5, K17 and K18 between BRET ligand binding and radioligand binding assays, we demonstrate the order of affinity remains consistent. Indeed, this was seen across all three experimental approached: Schild analysis, NanoBRET ligand-binding assay and radioligand binding.
Our MD simulations showed the potential binding site of K18, our most potent and selective A3R antagonist, within the A3R orthosteric binding area (Fig. 4A). Here, K18 is stabilised through hydrogen bonding interactions between the amino group and thiazole ring of the ligand and the amide side chain of N250 6.55 . In addition, the dichloro-phenyl ring can be oriented to the unique lipophilic area of A3R including V169 5.30 , M177 5.40 , I249 6.54 and L264 7.34 stabilized in that cleft through attractive van der Waals interactions; K18 is further stabilized through π-π aromatic stacking interactions between isoxazole ring and the phenyl group of F168 5.29 and the thiazole group is oriented deeper into the receptor favoring interactions with L246 6.51 and L90 3.32 and possibly with I268 7.39 . In combination with our mutagenesis data, the final binding pose of K18 appears to be within the orthosteric binding site, involving residues previously described to be involved in binding of A3R compounds (Arruda et al., 2017). We reported no detectable Gi/o response following co-stimulation with forskolin and NECA or IB-MECA for A3R mutants F168A 5.29 , L246A 6.51 , N250A 6.55 and I268A 7.39 (Stamatis et al., 2019, in preparation). These findings are in line with previous mutagenesis studies investigating residues important for agonist and antagonist binding at the human A3R (Gao et al., 2002, May et al., 2012. L90A 3.32 , V169A 5.30 , M177A 5.40 , I249A 6.54 and L264A 7.34 A3R all showed a detectable Gi/o response when stimulated with agonists (Stamatis et al., 2019).
Through performing Schild analysis, we were able to experimentally determine the effect of receptor mutation on antagonist affinity for L90A 3.32 , V169A/E 5.30 , M177A 5.40 , I249A 6.54 and L264A 7.34 A3R. The pA2 value for I249A 6.54 A3R was similar to WT, whereas M177A 5.40 and V169A 5.30 were significantly smaller suggesting these residues appear to be involved in K18 binding. Interestingly we found an increase in the pA2 for L90A 3.32 and L264A 7.34 when compared to WT, suggesting an enhanced ability of K18 to act as an antagonist. Further evidence was provided by the MM-PBSA calculations which were in agreement, based on the proposed binding model, between the calculated binding free energy by congeners of K18 having one or no chlorine atoms, i.e. compounds K17 and K5, and binding affinities and antagonistic potency. Importantly, substitution of the 1,3-thiazole ring in K17 with either a 2pyridyl ring (K32) or a 3-pyridyl ring (K10) but not a 4--pyridyl ring (K11) maintained A3R antagonistic potency. Although we have not directly determined the effects of similar pyridyl ring subsitutions on the higher affinity antagonist K18, we suspect there would be no significant increase in the potency of K18 given the small changes we observed for K17.
In conclusion, through pharmacological characterisation of a number of potential A3R antagonists, this study has determined K18 as a specific (<1 µM) A3R competitive antagonist. Our mutagenic studies, supported by MD simulations, identified the residues important for K18 binding are located within the orthosteric site of the A3R. Importantly, the absence of a chloro substituent, as is the case in K5, led to affinity loss. We suggest that the high affinity subtype selectivity of K18 makes it a molecule to begin detailed SAR and represents a useful tool compound that warrants further assessment for its therapeutic potential.

Cell culture and Transfection
Cell lines were maintained using standard subculturing routines as guided by the European Collection of Cell Culture (ECACC) and checked annually for mycoplasma infection using an EZ-PCR mycoplasma test kit from Biological Industries (Kibbutz Beit-Haemek, Israel). All procedures were performed in a sterile tissue culture hood using aseptic technique and solutions used in the propagation of each cell line were sterile and pre-warmed to 37 o C. All cells were maintained at 37°C with 5% CO2, in a humidified atmosphere. CHO-K1-A1R or CHO-K1-A3R cells were routinely cultured in Hams F-12 nutrient mix (21765029, Thermo Fisher Scientific) supplemented with 10% Foetal bovine serum (FBS) (F9665, Sigma-Aldrich). Flp-In-CHO cells purchased from Thermo Fisher Scientific (R75807) were maintained in Hams F-12 nutrient mix supplemented with 10% FBS containing 100 μg/mL Zeocin TM Selection Antibiotic (Thermo Fisher Scientific).
Stable Flp-In-CHO cell lines were generated through co-transfection of the pcDNA5/FRT expression vector (Thermo Fisher Scientific) containing the gene of interest and the Flp recombinase expressing plasmid, pOG44 (Thermo Fisher Scientific). Transfection of cells seeded in a T25 flask at a confluency of ≥80% was performed using TransIT®-CHO Transfection Kit (MIR 2174, Mirus Bio), in accordance with the manufacturer's instructions. Here, a total of 6 μg of DNA (receptor to pOG44 ratio of 1:9) was transfected per flask at a DNA:Mirus reagent ratio of 1:3 (w/v). 48 hours post-transfection, selection using 600 μg/mL hygromycin B (Thermo Fisher Scientific), concentration determined through preforming a kill curve, was performed for two days prior to transferring the cells into a fresh T25 flask. Stable Flp-In-CHO cell lines expressing the receptor of interest were selected using 600 μg/mL hygromycin B whereby the media was changed every two days. Successful mutant cell line generation for non-signalling mutants were confirmed by Zeocin TM sensitivity (100 μg/mL).

Constructs
The human A3R originally in pcDNA3.1+ (ADRA3000000, cdna.org) was cloned into the pcDNA5/FRT expression vector and co-transfected with pOG44 to generate a stable Flp-In-CHO cell line. Mutations within the A3R were made using the QuikChange Lightening Site-Directed Mutagenesis Kit (Agilent Technologies) in accordance with the manufacturer's instructions. All oligonucleotides used for mutagenesis were designed using the online Agilent Genomics 'QuikChange Primer Design' tool (Supplementary Table 7) and purchased from Merck. All constructs were confirmed by in-house Sanger sequencing. were purchased from Sigma-Aldrich and dissolved in dimethyl-sulphoxide (DMSO). PMA was purchased from Sigma-Aldrich. Compounds under investigation were purchased from emolecules and dissolved in DMSO.

cAMP accumulation assay
For cAMP accumulation (A2AR and A2BR) or inhibition (A1R or A3R) experiments, cells were harvested and re-suspended in stimulation buffer (PBS containing 0.1% BSA and 25 μM rolipram) and seeded at a density of 2,000 cells per well of a white 384-well Optiplate and stimulated for 30 minutes with a range of agonist concentrations. In order to allow the A1R/A3R mediated Gi/o response to be determined, co-stimulation with forskolin, an activator of AC (Zhang et al., 1997), at the indicated concentration (depending on cell line) was performed. For testing of potential antagonists, cells received a co-stimulation stimulated with forskolin, agonist and compound/DMSO control. cAMP levels were then determined using a LANCE® cAMP kit as described previously (Knight et al, 2016). In order to reduce evaporation of small volumes, the plate was sealed with a ThermalSeal® film (EXCEL Scientific) at all stages.
Phospho-ERK assay ERK1/2 phosphorylation was measured using the homogeneous time resolved fluorescence (HTRF)® Phospho-ERK (T202/Y204) Cellular Assay Kit (Cisbio Bioassays, Codolet, France) two-plate format in accordance with the manufacturer's instructions. A3R expressing Flp-In-CHO were seeded at a density of 2,000 cells per well of a white 384-well Optiplate and stimulated with agonist and test compounds for 5 minutes at 37°C. Plate reading was conducted using a Mithras LB 940 (Berthold technology). All results were normalised to 5 minutes stimulation with 1 μM PMA, a direct protein kinase C (PKC) activator (Jiang and Fleet, 2012). To determine if the measured pERK1/2 level was Gi-mediated, we treated cells with Pertussis toxin (PTX) (Tocris Biosciences) for 16 hours at 100 ng/mL prior to pERK assay.

Radioligand Binding
All pharmacological methods followed the procedures as described in the literature (Klotz et al., 1998). In brief, membranes for radioligand binding were prepared from CHO cells stably transfected with hAR subtypes in a two-step procedure. In the first step, cell fragments and nuclei were removed at 1000 x g and then the crude membrane fraction was sedimented from the supernatant at 100000 x g. The membrane pellet was resuspended in the buffer used for the respective binding experiments and it was frozen in liquid nitrogen and stored at -80ºC. For radioligand binding at the A1R, 1 nM [ 3 H]CCPA was used, for A2AR 10 nM [ 3 H]NECA and

for A3R 1 nM [ 3 H]HEMADO. Non-specific binding of [ 3 H]CCPA was determined in the presence of 1 mM theophylline and in the case of [ 3 H]NECA and [ 3 H]HEMADO 100 μM R-PIA was used.
Ki values from competition experiments were calculated using Prism (GraphPad Software, La Jolla, CA, U.S.A.) assuming competitive interaction with a single binding site. The curve fitting results (see Fig. 8 in Lagarias et al. 2018) showed R 2 values ≥ 0.99 for all compounds and receptors, indicating that the used one-site competition model assuming a Hill slope of n=1 was appropriate.

Determining Kon and Koff rates of A3R antagonists
Through the use of NanoBRET, real-time quantitative pharmacology of ligand-receptor interactions can be investigated in living cells. CA200645, a high affinity AR xanthine amine congener (XAC) derivative containing a polyamide linker connected to the BY630 fluorophore, acts as a fluorescent antagonist at both A1R and A3R with a slow off-rate (Stoddart et al., 2012). Using an N-terminally NanoLuc (Nluc)-tagged A3R expressing cell line, competition binding assays were conducted. The kinetic data was fitted with the 'kinetic of competitive binding' model (Motulsky and Mahan, 1984; built into Prism) to determine affinity (pKi) values and the association rate constant (Kon) and dissociation rates (Koff) for unlabelled A3R antagonists. This model resulted in several cases in an ambiguous fit (Regression with Prism 8: "Ambiguous", 2019). We developed a new model which expands on the 'kinetic of competitive binding' model to accommodate very rapid competitor dissociation, assuming the unlabelled ligand rapidly equilibrates with the free receptor. This method allows determination of compound affinity (pKi) from the kinetic data.
Filtered light emission at 450 nm and > 610 nm (640-685 nm band pass filter) was measured using a Mithras LB 940 and the raw BRET ratio calculated by dividing the 610 nm emission with the 450 nm emission. Here, Nluc on the N-terminus of A3R acted as the BRET donor (luciferase oxidizing its substrate) and CA200645 acted as the fluorescent acceptor. CA200645 was used at 25 nM, as previously reported (Stoddart et al., 2015). BRET was measured following the addition of the Nluc substrate, furimazine (0.1 µM). Nonspecific binding was determined using a high concentration of unlabelled antagonist, MRS 1220 (10 nM), for Nluc-A3R. Data were also analysed using an equation that assumes compound dissociation is too rapid for the dissociation rate constant to be determined reliably and the fits to the two equations compared ("Kinetics of competitive binding, rapid competitor dissociation", derived in the Appendix I, Supplementary material). This equation assumes rapid equilibration between compound and receptor and consequently provides an estimate of the equilibrium binding affinity of the compound (Ki) but not the binding kinetics of the compound. The equation is,

Receptor binding kinetics data analysis
where rI is fractional occupancy of receptors not bound by L: and kobs,+ I is the observed association rate of tracer in the presence of competitor, defined as, ( HIJ,LG = [#]( ) (1 − F G ) + ( @ The fits to the two equations were compared statistically using a partial F-test in Prism 8 (Motulsky, 2019).

Data and Statistical analysis
All in vitro assay data was analysed using Prism 8.0 (GraphPad software, San Diego, CA), with all dose-inhibition or response curves being fitted using a 3-parameter logistic equation to calculate response range or Emax and IC/EC50. Dose-inhibition/dose-response curves were normalised to either forskolin response or forskolin inhibition (A1R and A3R), relative to NECA/IB-MECA. In the case of pERK1/2 response, normalisation was performed to PMA.
Schild analysis was performed to obtain pA2 values (the negative logarithm to base 10 of the molar concentration of an antagonist that makes it necessary to double the concentration of the agonist to elicit the original submaximal response obtained by agonist alone (Schild, 1947)) for the potential antagonists. In cases where the Schild slope did not differ significantly from unity, the slope was constrained to unity giving an estimate of antagonist affinity (pK B). pA2 and pK B coincide when the slope is exactly unity. The pA2 values obtained through conducting Schild analysis of K18 at WT and mutant A3R were compared in order to indicate important residues involved in K18 binding.
The data and statistical analysis comply with the recommendations on experimental design and analysis in pharmacology (Curtis et al., 2018). Statistical significance (*, p< 0.05; **, p<0.01; ***, p<0.001; ****, p<0.0001) was calculated using a one-way ANOVA with a Dunnett's post-test for multiple comparisons or Students' t-test, as appropriate. Compounds taken forwards for further experiments after initial screening were identified as having the highest statistical significance (P value of 0.001 (***) or <0.0001 (****)). All statistical analysis was performed using Prism 8.0 on data which were acquired from experiments performed a minimum of five times, conducted in duplicate.

MD simulations
Preparation of the complexes between A3R and K5, K17 or K18 was based on a homology model of A2AR (see Appendix II in Supplementary material). Each ligand−protein complex was embedded in hydrated POPE bilayers. A simulation box of the protein-ligand complexes in POPE lipids, water and ions was built using the System Builder utility of Desmond (Desmond Molecular Dynamics System, version 3.0; D.E. Shaw Res. New York, 2011;Maest. Interoperability Tools, 3.1;Schrodinger Res. New York, 2012.). A buffered orthorhombic system in 10 Å distance from the solute atoms with periodic boundary conditions was constructed for all the complexes. The MD simulations were performed with Amber14 and each complex-bilayer system was processed by the LEaP module in AmberTools14 under the AMBER14 software package (Case et al., 2014). Amber ff14SB force field parameters (Maier et al., 2015) were applied to the protein, lipid14 to the lipids (Dickson et al., 2014), GAFF to the ligands (Wang et al., 2004) andTIP3P (Jorgensen et al., 1983) to the water molecules for the calculation of bonded, vdW parameters and electrostatic interactions. Atomic charges were computed according to the RESP procedure (Bayly et al., 1993) using Gaussian03 (Frisch et al., 2003 and antechamber of AmberTools14 (Case et al., 2014). The temperature of 310 K was used in MD simulations in order to ensure that the membrane state is above the main phase transition temperature of 298 K for POPE bilayers (Koynova and Caffrey, 1998). In the production phase, the relaxed systems were simulated in the NPT ensemble conditions for 100 ns. The visualization of produced trajectories and structures was performed using the programs Chimera (Pettersen et al., 2004) and VMD (Humphrey et al., 1996). All the MD simulations were run on GTX 1060 GPUs in lab workstations or on the ARIS Supercomputer.

MM-PBSA calculations
Relative binding free energies of the complexes between K5, K17 and K18 and A3R was estimated by the 1-trajectory MM-PBSA approach (Massova and Kollman, 2000). Effective binding energies (ΔGeff) were computed considering the gas phase energy and solvation free energy contributions to binding. For this, structural ensembles were extracted in intervals of 50 ps from the last 50 ns of the production simulations for each complex. Prior to the calculations all water molecules, ions, and lipids were removed, and the structures were positioned such that the geometric center of each complex was located at the coordinate origin. The polar part of the solvation free energy was determined by calculations using Poisson-Boltzmann (PB) calculations (Homeyer and Gohike, 2013). In these calculations, a dielectric constant of εsolute = 1 was assigned to the binding area and εsolute = 80 for water. Using an implicit solvent representation for the calculation of the effective binding energy is an approximation to reduce the computational cost of the calculations. The binding free energy for each complex was calculated using equation (1) ΔGeff = ΔEMM + ΔGsol (1) In equation (1) ΔGeff is the binding free energy for each calculated complex neglecting the effect of entropic contributions or assuming to be similar for the complexes studied. ΔEMM defines the interaction energy between the complex, the protein and the ligand as calculated by molecular mechanics in the gas phase. ΔGsol is the desolvation free energy for transferring the ligand from water in the binding area calculated using the PBSA model. The terms for each complex ΔEMM and ΔGsol are calculated using equations (2)    one of the remaining AR subtypes were exposed to forskolin in the case of Gi-coupled A1R and A3R (1 μM or 10 μM, respectively) or DMSO control in the case of Gs-coupled A2AR and A2BR, NECA and test compound (10 μM) for 30 min and cAMP accumulation detected. All values are mean ± SEM expressed as percentage forskolin inhibition (A1R and A3R) or stimulation (A2AR and A2BR), relative to NECA. n ≥ 3 independent experimental repeats, conducted in duplicate.   Table 3. IB-MECA stimulated cAMP inhibition at WT A3R: activity of MRS 1220 and potential antagonists. Forskolin stimulated cAMP inhibition was measured in Flp-In-CHO stably expressing A3R following stimulation with 10 µM forskolin, compound at the indicated concentration and varying concentrations of IB-MECA.  Side chains of critical residues for binding indicated from the MD simulations are shown in sticks. Residues L90 3.32 , V169 5.30 , M177 5.40 , I249 6.54 and L264 7.34 , in which carbon atoms are shown in grey, were confirmed experimentally; in residues F168 5.29 , L246 6.51 , I268 7.39 and N250 6.55 carbon atoms are shown in magenta; nitrogen, oxygen and sulfur atoms are shown in blue, red and yellow respectively. --47.0 ± 2.4 -8.8 ± 2.7 29.8 ± 2.9 -25.9 ± 3.6 6.35 6.33 ± 0.03 5.38 K18 -46.3 ± 2.9 -7.5 ± 2.4 26.9 ± 3.1 -26.9 ± 2.7 7.20 6.92 ± 0.10 6.05 a vdW energy of binding calculated using molecular mechanics b Electrostatic energy of binding calculated using molecular mechanics c Difference in solvation energy between the complex, the protein and the ligand, i.e. Gcomplex, solv -(Gprotein, solv + Gligand, solv) d Effective binding free energy calculated as ΔGeff = ΔEΜΜ + ΔGsol; in Table 4, ΔEΜΜ = ΕvdW + EEL (see Materials and Methods) e Equilibrium dissociation constant of MRS 1220, K5, K17 and K18 as determined through three independent experimental approaches: Schild analysis (pKB), NanoBRET (pKi) or radioligand binding (pKi) f pKB obtained through Schild analysis in A3R stably expressing Flp-In CHO cells g pKi (mean ± SEM) obtained in NanoBRET binding assays using Nluc-A3R stably expressing HEK 293 cells and determined through fitting our "Kinetics of competitive binding, rapid competitor dissociation" model or in the case of MRS 1220 through fitting with the 'Kinetics of competitive binding' model with a determined Kon (k3) and Koff (k4) rate of 3.25 ± 0.28 x 10 8 M -1 min -1 and 0.0248 ± 0.005 min -1 , respectively h pKi values previously published for K5, K17 and K18 (Lagarias et al., 2018) or MRS 1220 (Stoddart et al., 2015) through radioligand binding assays.