Structure-Based Design of Potent and Selective Ligands at the Four Adenosine Receptors

molecules Article Structure-Based Design of Potent and Selective Ligands at the Four Adenosine Receptors Willem Jespers 1, Ana Oliveira 1, Rubén Prieto-Díaz 2, María Majellaro 2, Johan Åqvist 1, Eddy Sotelo 2ID and Hugo Gutiérrez-de-Terán1,*ID 1Department of Cell and Molecular Biology, Uppsala University, Biomedical Centre (BMC), BOX 596, SE-751 24 Uppsala, Sweden; [email protected] (W.J.); [email protected] (A.O.); [email protected] (J.Å.) 2Centro Singular Investigación Quimica Biologica e Materiales Moleculares (CIQUS), Departamento de Quimica Orgánica, Facultade de Farmacia, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain; [email protected] (R.P.-D.); [email protected] (M.M.); [email protected] (E.S.) *Correspondence: [email protected]; Tel.: +46(0)18-471-5056 Received: 20 October 2017; Accepted: 8 November 2017; Published: 10 November 2017 Abstract: The four receptors that signal for adenosine, A 1 , A 2A , A 2B and A 3 ARs, belong to the superfamily of G protein-coupled receptors (GPCRs). They mediate a number of (patho)physiological functions and have attracted the interest of the biopharmaceutical sector for decades as potential drug targets. The many crystal structures of the A 2A , and lately the A 1 ARs, allow for the use of advanced computational, structure-based ligand design methodologies. Over the last decade, we have assessed the efficient synthesis of novel ligands specifically addressed to each of the four ARs. We herein review and update the results of this program with particular focus on molecular dynamics (MD) and free energy perturbation (FEP) protocols. The first in silico mutagenesis on the A 1 AR here reported allows understanding the specificity and high affinity of the xanthine-antagonist 8-Cyclopentyl-1,3-dipropylxanthine (DPCPX). On the A 2A AR, we demonstrate how FEP simulations can distinguish the conformational selectivity of a recent series of partial agonists. These novel results are complemented with the revision of the first series of enantiospecific antagonists on the A 2B AR, and the use of FEP as a tool for bioisosteric design on the A3AR. Keywords: free energy perturbation (FEP); G protein-coupled receptors (GPCRs); molecular dynamics (MD) simulations; structure-based drug design (SBDD) 1. Introduction The superfamily of G protein-coupled receptors (GPCRs), which encompasses targets for more than 30% of marketed drugs [ 1 ], was traditionally relegated from the field of structure-based drug design (SBDD) due to the inherent difficulty of solving the structure of these membrane receptors. However, recent advancements in membrane crystallography in the past decade have led to an explosive growth in available crystal structures, currently comprising 196 receptor–ligand complexes [ 2 ]. A particularly privileged family in this sense is the adenosine receptors (ARs), with several complexes of the A 2A and lately for the A 1 ARs deposited in the PDB. These structures can be used to model the remaining A 2B and A 3 receptor subtypes in the family, thus providing a full pallet of atomistic models of the AR family. ARs play an essential role in many physiological processes in the cardiovascular and central nervous systems, and control anti-inflammatory and immunosuppressive responses. Consequently, this family of receptors is of outstanding interest as drug targets for cardiovascular, neurodegenerative and autoimmune diseases, as well as cancer [3]. Molecules 2017,22, 1945; doi:10.3390/molecules22111945 www.mdpi.com/journal/molecules Molecules 2017,22, 1945 2 of 17 Over the last decades, a large number of agonists and antagonists have been synthetized and pharmacologically characterized, accumulating to 25,000 entries in the ChEMBL database, which have been source of a number of computational screening and SAR analyses [ 4 – 7 ]. The exponential growth of experimental GPCR structures [ 8 ] proved extremely beneficial for the structural biology and structure-based ligand design of ARs [ 9 , 10 ]. The current repertoire, consisting of several A 2A AR structures in complex with antagonists and agonists, was recently complemented with fully active conformation of this receptor in ternary complex with a G protein mimic [ 11 ], and two antagonist-bound (inactive) A 1 AR structures [ 12 , 13 ]. The effects of point mutations on ligand binding constitute another important information resource for the characterization of receptor–ligand interactions and receptor activation. Thus, site-directed mutagenesis data further complement the structural and chemical information, and over 2500 data points are deposited in the GPCRdb [ 2 ]. This combination of structural, pharmacological and chemical information allowed the pharmaceutical SBDD of 1,2,4 triazines as A 2A AR antagonists, one of which reached clinical studies for the treatment of Parkinson’s disease [9]. The family of ARs, and in particular the A 2A AR, have become targets for several computational studies to evaluate methods and protocols proposed in the field of GPCR ligand design [ 10 ]. One of the first topics investigated after the release of the first A 2A AR structure was the actual role of water molecules in ligand binding, studied through molecular dynamics (MD) sampling of the solvent in the binding cavity [ 14 ]. This allowed the identification of water clusters whose displacement should favor the binding affinity of prospective ligands [ 15 ], or of structural waters that should be considered for antagonist docking and virtual screening (VS) on this receptor [ 16 ]. The detailed structural characterization of the A 2A AR also allowed for the estimation of binding free energies, based on MD sampling. The low affinity of the micromolar antagonist caffeine was rationalized using the molecular mechanics/Poison-Boltzman surface area (MM/PBSA) method [ 17 ], introducing the idea of a multiple binding mode for this molecule, which was partially supported by recent crystal structures with xanthine-like antagonists [ 13 , 18 ]. In this area, our group developed a protocol based on free energy perturbation (FEP) for the estimation of relative binding affinities upon point mutations [ 19 ]. The results of this “in-silico” site directed mutagenesis can be directly compared to experimentally determined ligand affinity ratios between wild-type (WT) and mutant receptor variants, coming from the numerous site-directed mutagenesis studies as we will discuss here for the case of the A 2A AR [ 20 , 21 ]. The method, which was initially developed and later applied in other GPCR systems [ 19 , 22 ], overcomes the convergence problems associated with large perturbations and allows reproducing the effect of point mutations on ligand binding with high precision and convergence. During the last years, our labs have combined FEP and other computational methodologies with high-throughput synthetic methodologies and pharmacological characterization, to develop various series of structurally simple, novel AR ligands [ 23 ]. Our SBDD protocol includes homology modeling, protein–ligand docking, 3D-QSAR, MD and FEP simulations, and was used to assist the design of novel ligands as well as to predict or rationalize their pharmacological profile [ 24 – 32 ]. The different series obtained with this strategy, represented in Figure 1, have been optimized to yield potent and selective antagonists of the different ARs. The first scaffold reported, a 2,4-diarylpyrimidine, was recently followed by a bioisosteric series of 2,4-diarylpyridines, both of them revealing high affinity and selectivity as A 3 AR antagonists [ 24 , 25 , 30 ]. For the A 2B receptor, we have documented the first non-planar, enantioselective antagonists, and disclosed the reasons of their enantiospecificity [ 31 , 32 ]. The A 2A AR emerges as the most attractive target for SBDD within the AR family, due to the existence of several crystal structures complemented with broad mutagenesis and SAR data. Consequently, we used this system to train our protocol for in-silico mutagenesis, based on well-converged FEP calculations [20,21] . Moreover, we went one step further in the design of new ligands and isolated, from a series of antagonists, a group of compounds that presented a pharmacological profile of partial agonists, where we identified a prolinol moiety as a replacement of the ribose group of classical agonists [ 28 ]. Finally, recently published crystal structures of the A 1 AR allow for the characterization Molecules 2017,22, 1945 3 of 17 of molecular determinants of high affinity, which we here exemplify with the calculation and interpretation of the effect of point mutations on the binding affinity of the reference xanthine-like antagonist DPCPX (see Figure 1). In the following sections, we review our recent results on the design and characterization of ligand binding for each of the four ARs, and further complement these studies with new analysis and calculations. We put the results of this project in the broad perspective of ARs ligand design, and discuss the methodological developments in our SBDD pipeline, with an impact in the design of novel ligands for this and other families of receptors. Molecules 2017, 22, 1945 3 of 17 In the following sections, we review our recent results on the design and characterization of ligand binding for each of the four ARs, and further complement these studies with new analysis and calculations. We put the results of this project in the broad perspective of ARs ligand design, and discuss the methodological developments in our SBDD pipeline, with an impact in the design of novel ligands for this and other families of receptors. Figure 1. Adenosine receptor antagonists discussed in this work. 2. Results 2.1. A 2A AR The A 2A AR is one of the best-characterized GPCRs, with various entries in the PDB including inactive structures, in complex with eight different antagonists, four active-like structures with agonists, and one fully-activated A 2A AR, in a ternary complex with an agonist and a mimic of the intracellular G protein [11]. In addition, this receptor accumulates 42% of the mutagenesis data deposited in the GPCRdb for the ARs, while 31% of the 26.274 entries in ChEMBL [33] (release 23) report affinities for this receptor subtype. The vast amount of chemical and structural information corresponds to the biopharmaceutical interest in this receptor, where antagonists have reached clinical trials as a nondopaminergic therapy for Parkinson’s disease [34], while the agonist Regadenoson has been approved by the FDA as a coronary vasodilator [35]. Given this unique combination of mutagenesis and structural information, we have used this system to further develop our in silico mutagenesis approach based on FEP [20,21]. The thermodynamic cycle underlying this approach was originally developed by Kollman and co-workers to characterize the effect of a single mutation on ligand binding free energies (and catalysis) [36], and consists on the perturbation of the amino acid sidechain from WT to mutant in both the complex and the apo versions of the receptor (which acts therefore as the reference state). Starting from this concept, we developed a dedicated FEP protocol to routinely evaluate affinity shifts for a given ligand to mutant versions of the receptor, as compared to the WT. Although standard FEP implementations might work well enough for simple cases (e.g., Serine to Alanine), more drastic mutations (e.g., Tryptophan to Alanine) are subject to convergence and sampling problems, usually leading to high hysteresis values (defined as the difference in the forward and backwards pathway of the sidechain transformation). In our implementation, described in detail in the methods section and applied here to the original characterization of mutagenesis experiments on the A 1 AR (see below), the sidechain of the residue to mutate is gradually annihilated to alanine in a stepwise fashion. With this strategy, we examined the effect of 18 alanine mutations reported for the A 2A AR for agonist and antagonist binding with high precision (s.e.m. = 0.7 kcal/mol) and convergence values (average hysteresis of 0.2 kcal/mol) [21]. This analysis revealed effects that were not evident by the sole analysis of crystal structures, such as the disruption of water-mediated interactions causing decreased antagonist Figure 1. Adenosine receptor antagonists discussed in this work. 2. Results 2.1. A2AAR The A 2A AR is one of the best-characterized GPCRs, with various entries in the PDB including inactive structures, in complex with eight different antagonists, four active-like structures with agonists, and one fully-activated A 2A AR, in a ternary complex with an agonist and a mimic of the intracellular G protein [ 11 ]. In addition, this receptor accumulates 42% of the mutagenesis data deposited in the GPCRdb for the ARs, while 31% of the 26.274 entries in ChEMBL [ 33 ] (release 23) report affinities for this receptor subtype. The vast amount of chemical and structural information corresponds to the biopharmaceutical interest in this receptor, where antagonists have reached clinical trials as a nondopaminergic therapy for Parkinson’s disease [ 34 ], while the agonist Regadenoson has been approved by the FDA as a coronary vasodilator [35]. Given this unique combination of mutagenesis and structural information, we have used this system to further develop our in silico mutagenesis approach based on FEP [ 20 , 21 ]. The thermodynamic cycle underlying this approach was originally developed by Kollman and co-workers to characterize the effect of a single mutation on ligand binding free energies (and catalysis) [ 36 ], and consists on the perturbation of the amino acid sidechain from WT to mutant in both the complex and the apo versions of the receptor (which acts therefore as the reference state). Starting from this concept, we developed a dedicated FEP protocol to routinely evaluate affinity shifts for a given ligand to mutant versions of the receptor, as compared to the WT. Although standard FEP implementations might work well enough for simple cases (e.g., Serine to Alanine), more drastic mutations (e.g., Tryptophan to Alanine) are subject to convergence and sampling problems, usually leading to high hysteresis values (defined as the difference in the forward and backwards pathway of the sidechain transformation). In our implementation, described in detail in the methods section and applied here to the original characterization of mutagenesis experiments on the A 1 AR (see below), the sidechain of the residue to mutate is gradually annihilated to alanine in a stepwise fashion. With this strategy, we examined the Molecules 2017,22, 1945 4 of 17 effect of 18 alanine mutations reported for the A 2A AR for agonist and antagonist binding with high precision (s.e.m. = 0.7 kcal/mol) and convergence values (average hysteresis of 0.2 kcal/mol) [ 21 ]. This analysis revealed effects that were not evident by the sole analysis of crystal structures, such as the disruption of water-mediated interactions causing decreased antagonist binding observed in mutants V84 3.32 A, L294 6.51 A or Ile274 7.39 A (Ballesteros–Weinstein nomenclature for GPCRs is indicated in superscripts as X.YY, with X indicating the helix and YY the correlative position relative to the most conserved residue in the helix, when YY = 50 [ 37 ]). Another outcome of such analysis was deciphering the detrimental effect on agonist binding of the M177 5.38 A mutation due to an increased solvation of the binding site that perturbed agonist-receptor interactions. The methodology was then extended to account for non-alanine mutations, by joining two thermodynamic cycles describing the reduction of a sidechain, i.e., from WT and mutant, respectively, to a common fragment (e.g., Alanine). This allowed us to further characterize the effect of 17 non-alanine mutations on agonist binding, again in excellent agreement with the experimental data [20]. The most traditional and extensive use of FEP calculations on ligand design, however, is to predict the relative binding affinity between pairs of ligands. In this case, the application involves a thermodynamic cycle where ligand 1 is transformed into ligand 2 in both the bound and the reference (i.e., water solvated) states. An example of the application of this technique on antagonist design can be found in the section covering the A 3 AR (see below). We herein decided to mix the idea of a ligand-perturbation with a receptor-comparison in a new thermodynamic cycle, which we apply to identify the preferred receptor conformational state for different pharmacological classes of ligands. In theory, the capacity of a ligand to activate the receptor (i.e., agonist potency) can be related to its specific preference for the active conformation. Therefore, the relative affinities between two ligands for the inactive and active conformations of the receptor can be qualitatively related to their relative difference in potency. To illustrate this point, we here perform such calculations on our recently reported series of thiazolo[5,4-d]pirimidines, designed as potential agonists of the A 2A AR as bioisosteres of the purine ring on the basis of molecular docking and superposition with the classical agonists (e.g., adenosine, NECA) [ 28 ]. The molecules in this series that contain a 7-prolinol substitution were characterized as a new class of non-nucleoside A 2A AR partial agonists, confirming the hypothesis of the design where the 2-hydroxymethylene moiety of the prolinol would mimic the interactions of the 2 0 /3 0 hydroxyl groups of the ribose of agonists, observed in all agonist bound A 2A AR crystal structures. To understand the molecular determinants of the functional behavior of these molecules, we designed a thermodynamic cycle where the partial agonist ligand 10n (containing prolinol, relative efficacy as compared to NECA > 50%) was transformed in its antagonist analog lacking the 2-hydroxymethylene group ( 10m , with a relative efficacy value < 10%), both in the inactive and active-like A 2A AR structures as shown in Figure 2. The results indicate that the loss of the 2-hydroxymethylene “agonist tag” is more unfavorable in the active-like than in the inactive receptor structure ( ∆∆G10m→10n = 1.74 ±0.29 kcal/mol ), which is in qualitative agreement with the observed reduction in potency between this pair of compounds. Interestingly, the main difference in the simulations between the inactive and active structures are observed in the interactions of the 2-hydroxymethylene moiety, where the interaction with S277 7.42 in the (unfavored) inactive conformation is replaced by a hydrogen bond to H278 7.43 in the simulations of the active-like structure. This might explain why there are several partial agonists that exert a higher potency on the S277A 7.42 mutant receptor [ 38 ]. If followed a similar evaluation of a second pair of compounds within this series with similar results: The partial agonist denoted as 10d (40% efficacy as compared to NECA) contains a different substitution pattern on the ring exposed to the extracellular loops (i.e., 4-Cl-phenyl instead of 3,4,5-tris-OMe-phenyl in 10m / 10n , see Figure 1), and is transformed into the corresponding antagonist ( 10c , 10% efficacy) upon removal of the prolinol moiety. The calculated effect in this pair of ligands also indicates that the affinity of the partial agonist for the active form of the receptor is favored as compared to the antagonist. Moreover, the reduced energetic barrier as compared to the 10m / 10n Molecules 2017,22, 1945 5 of 17 pair ( ∆∆ G 10d→10c = 0.65 ± 0.26 kcal/mol) agrees with the lower experimental difference in the efficacy between the 10d/10c pair of compounds. Molecules 2017, 22, 1945 5 of 17 Figure 2. Thermodynamic cycle showing the relationship between the potency of the ligand pair 10m/10n and the relative free energies of perturbing 10n to 10m in the active-like (blue) and inactive (orange) structures. The grey sticks represent the conformation of the residues in the active-like structure as compared to the inactive structure, and vice versa. A ligand that activates the receptor (partial or full agonist) should in theory have a higher affinity for the active state of that receptor. Thus, perturbing the (partial) agonist (10n) into an antagonist (10m) should result in an unfavorable change in binding free energies, as shown here for the 10n→10m transformation. 2.2. A1AR The adenosine A1 receptor, which plays a significant role in the regulation of neural, renal and cardiac systems [3], was the first AR to be characterized [39]. A1AR antagonists have been predominantly derivatives of the xanthine scaffold, such as the reference antagonist DPCPX (see Figure 1), allowing alkyl substitutions on the 1 and 3 positions while more bulky, cyclic substituents on the 8 positions proved essential for subtype selectivity. Initial mutagenesis studies indicated an essential role of T2707.35 (Met in A2AAR and A2BAR, Leu in A3AR) in the accommodation of substituents at position 8. This was confirmed by two recently solved crystal structures in complex with selective xanthine derivatives [12,13]. Besides revealing T2707.35 as a selectivity determinant, the structures showed a slightly widened area as compared to the A2AAR in the bottom of the binding pocket, selectively Figure 2. Thermodynamic cycle showing the relationship between the potency of the ligand pair 10m / 10n and the relative free energies of perturbing 10n to 10m in the active-like (blue) and inactive (orange) structures. The grey sticks represent the conformation of the residues in the active-like structure as compared to the inactive structure, and vice versa. A ligand that activates the receptor (partial or full agonist) should in theory have a higher affinity for the active state of that receptor. Thus, perturbing the (partial) agonist ( 10n ) into an antagonist ( 10m ) should result in an unfavorable change in binding free energies, as shown here for the 10n→10m transformation. 2.2. A1AR The adenosine A 1 receptor, which plays a significant role in the regulation of neural, renal and cardiac systems [ 3 ], was the first AR to be characterized [ 39 ]. A 1 AR antagonists have been predominantly derivatives of the xanthine scaffold, such as the reference antagonist DPCPX (see Figure 1), allowing alkyl substitutions on the 1 and 3 positions while more bulky, cyclic substituents on the 8 positions proved essential for subtype selectivity. Initial mutagenesis studies indicated an essential role of T270 7.35 (Met in A 2A AR and A 2B AR, Leu in A 3 AR) in the accommodation of Molecules 2017,22, 1945 6 of 17 substituents at position 8. This was confirmed by two recently solved crystal structures in complex with selective xanthine derivatives [ 12 , 13 ]. Besides revealing T270 7.35 as a selectivity determinant, the structures showed a slightly widened area as compared to the A 2A AR in the bottom of the binding pocket, selectively accommodating the substituents at position 1. Finally, the binding pocket of the A 1 AR was shown to be unexpectedly large. All these differences deserve a detailed exploration on the basis of MD and free energy calculations, to understand the specific requirements of ligand binding to this receptor. We here apply our in silico mutagenesis approach, described below for the case of A 2A AR co-crystallized ligands [ 20 , 21 ], to understand the structural requirements of binding of the selective xanthine DPCPX to the A 1 AR. The first stage was to manually dock DPCPX onto the A 1 AR crystal structure 5N2S, based on the structural similarity of the co-crystalized ligand PSB36 and DPCPX (see Figure 3and Methods). A double hydrogen bond with N254 6.55 formed the main anchor point between the core scaffold and the receptor, similar to many other crystal structure complexes of A 1 and A2A ARs. Molecules 2017, 22, 1945 6 of 17 accommodating the substituents at position 1. Finally, the binding pocket of the A1AR was shown to be unexpectedly large. All these differences deserve a detailed exploration on the basis of MD and free energy calculations, to understand the specific requirements of ligand binding to this receptor. We here apply our in silico mutagenesis approach, described below for the case of A2AAR co-crystallized ligands [20,21], to understand the structural requirements of binding of the selective xanthine DPCPX to the A1AR. The first stage was to manually dock DPCPX onto the A1AR crystal structure 5N2S, based on the structural similarity of the co-crystalized ligand PSB36 and DPCPX (see Figure 3 and Methods). A double hydrogen bond with N254 6.55 formed the main anchor point between the core scaffold and the receptor, similar to many other crystal structure complexes of A1 and A2A ARs. (A) (B) Figure 3. (A) Average structure from MD simulations of the DPCPX-A1AR complex (blue) superimposed with the crystal structure of the same receptor in complex with PSB-036 (PDB 5N2S); and (B) close look in the binding site, highlighting selected residues for in silico mutagenesis studies (see Table 1). Although no affinity data is available for DPCPX for the N254A6.55 mutant A1AR, the fact that many other ligands show severely reduced affinities for N6.55 mutant ARs (37 receptor–ligand pairs in the GPCRdb) indicate that this residue is essential for ligand binding. To investigate this hypothesis, we mutated this residue to alanine. As shown in Table 1, the mutation led to unfavorable binding affinities for the mutant receptor. This effect is of similar magnitude as for the F171AEL2 mutation, which was shown in literature to result in completely abolished binding [40]. Mutational effects of other residues within the binding site were also correctly reproduced (Table 1). Particularly interesting is T277A7.41, the only mutation resulting in favorable binding affinities as compared to WT A1AR. This mutation has been thoroughly studied in ligand binding, and is included as a thermostabilizing mutant in the crystal structure 5N2S [13]. Although we correctly capture the favorable effect on ligand binding, there is an overestimation of the effect of about 2 kcal/mol. Therefore, we ran the same mutation on another crystal structure (5UEN) of the A1AR, which does not include this thermostabilizing mutation. This time the calculated effect is milder, in agreement with the experimental data (Table 1). Figure 3. ( A ) Average structure from MD simulations of the DPCPX-A 1 AR complex (blue) superimposed with the crystal structure of the same receptor in complex with PSB-036 (PDB 5N2S); and ( B ) close look in the binding site, highlighting selected residues for in silico mutagenesis studies (see Table 1). Although no affinity data is available for DPCPX for the N254A 6.55 mutant A 1 AR, the fact that many other ligands show severely reduced affinities for N 6.55 mutant ARs (37 receptor–ligand pairs in the GPCRdb) indicate that this residue is essential for ligand binding. To investigate this hypothesis, we mutated this residue to alanine. As shown in Table 1, the mutation led to unfavorable binding affinities for the mutant receptor. This effect is of similar magnitude as for the F171A EL2 mutation, which was shown in literature to result in completely abolished binding [ 40 ]. Mutational effects of other residues within the binding site were also correctly reproduced (Table 1). Particularly interesting is T277A 7.41 , the only mutation resulting in favorable binding affinities as compared to WT A 1 AR. This mutation has been thoroughly studied in ligand binding, and is included as a thermostabilizing mutant in the crystal structure 5N2S [ 13 ]. Although we correctly capture the favorable effect on ligand binding, there is an overestimation of the effect of about 2 kcal/mol. Therefore, we ran the same mutation on another crystal structure (5UEN) of the A 1 AR, which does not include this thermostabilizing mutation. This time the calculated effect is milder, in agreement with the experimental data (Table 1). Molecules 2017,22, 1945 7 of 17 Table 1. Mutational effects on DPCPX affinity for A 1 AR. The experimental affinities were retrieved from the GPCRdb and the reference(s) are given in brackets. The binding affinity (pK D ) was converted to ∆∆G following ∆∆G = RT·ln (KDMutant/KDWT). Mutation ∆∆G (kcal/mol) In Vitro In Silico F171AEL2 [40]4.32 a6.15 ±0.69 I175AEL2 [40]0.98 0.26 ±0.39 M177A5.37 [40]1.06 1.84 ±0.46 N254A6.55 ND b4.14 ±0.67 T270A7.34 [40]0.46 0.40 ±0.48 T277A7.41 [41–43]−0.32 ±0.19 c−2.77 ±0.55 −0.42 ±0.54 d a No detectable binding, the value represents the detection threshold of the experiment; b No experimental value determined in literature; c An average value and associated s.e.m. were calculated based on the reported values from literature (n= 3); dCalculations performed on the 5UEN crystal structure. 2.3. A2BAR The A 2B AR is a receptor with low affinity for the natural agonist. It remains silent under physiological conditions to be rapidly activated during chronic highly oxidative stress conditions (e.g., hyperglycemia or mast cell activation). For these reasons, it is an attractive target for inflammatory processes like asthma or colitis, as well as for other diseases like diabetic retinopathy or cancer [ 1 ]. However, because of its poor expression levels and low affinity for standard ligands, the pharmacological characterization of this receptor is scarce as compared with the rest of the members of the ARs family. We are consequently engaged in the development of potent and selective antagonists as novel pharmacological tools to aid in this process. The hypothesis underlying this project was the idea that ligands with non-planar, stereodiverse topologies would provide additional selective tools, due to more specific stereoselective ligand-target interactions. We initially designed the 3,4-dihydropyrimidin-2(1H)-one scaffold, which is easily assembled via the Biginelli multicomponent reaction, an advantage when it comes to develop series of compounds and evaluate their SAR [ 29 ]. The promising affinity of some compounds in the initial series, measured as a racemic mixture, was explored by means of molecular docking on our A 2B AR homology model. The model revealed an enantiospecific binding orientation, suggesting the potential stereoselectivity of this scaffold and identifying the potentially active isomer (Figure 4). These promising results led to a further exploration of non-planar derived scaffolds, by structural diversification at face 2 of the diazinone scaffold. Three series of analogues were developed, namely 3-deaza[pyridin-2(1H)-ones] and bicyclic and tricyclic derivatives fused at positions 2,3 or 5,6 of the heterocyclic framework, identifying selective and potent A 2B AR ligands (Figure 1) where the hypothesis of the stereoselectivity was maintained [ 32 ]. In a later modification of the original scaffold, we created a series of 2-cyanoimino-4-substituted6-methyl-1,2,3,4-tetrahydro-pyrimidine-5-carboxylates [ 31 ]. This time the two enantiomers of the most attractive ligand ( 16b ) were separated by chiral HPLC and their absolute configurations established by circular dichroism (Figure 4). The biological evaluation of the two enantiomers demonstrated that the affinity is exclusively due to the (S)- 16b enantiomer, validating the prediction from our molecular modeling studies. Molecules 2017,22, 1945 8 of 17 Molecules 2017, 22, 1945 8 of 17 Figure 4. Synthesis and enantiospecific binding characterization of compound 16b as antagonist of the A 2B AR: (A) the Biginelli multicomponent reaction allows the assembly of the 3,4-dihydropyrimidine scaffold bearing different diversity points; (B) the chiral resolution identified the S isomer as the one responsible of the biological affinity of the racemic mixture; and (C) the molecular modeling confirmed the binding mode previously hypothesized for isomer (S)-16b. The water-density maps (black mesh) are calculated here by means of MD exploration (grid spacing 1 Å, occupancy > 80%), and are overlaid with the position of structural waters previously considered during the docking of this series (red spheres). The docking studies of these series of non-planar scaffolds were complemented by MD simulations, which allowed characterizing the role of water molecules in ligand binding. In a first stage, the selected complex of A 2B AR was equilibrated with our PyMemDyn protocol, as implemented in the GPCR-ModSim web server [44,45] (see Methods). From this analysis, we identified a water molecule that mediated the interaction between the NH in position 1 of the parent 3,4-dihydropyrimidin-2(1H)- ones and Glu169 EL2 , and a second water molecule in the region between the alkoxy substituent and the furan/thienyl ring, which is actually equivalent to a conserved water molecule observed in the crystal structures of A 2A AR with antagonist ZM241385 [23]. These pair of water molecules were maintained in successive docking runs for the expanded series of compounds, verifying their potential role in mediating this inter and intra molecular interactions [29,31,32]. We here provide further support for the existence of these structural water molecules, by means of an MD exploration of the solvent in the binding site (see methods). The idea is that the structural water molecules, which might favor the binding of ligands, is reflected in a higher occupancy during the MD trajectories, as we recently showed for the A 3 AR (see below and reference [30]). Figure 4 shows the calculated water density maps at 80% occupancy, supporting the position of the two water molecules discussed above, mediating interactions between the ligand and E170 EL2 as well as the intramolecular interaction between the 2-furan and esther groups of the ligand. 2.4. A 3 AR The identification of A 3 AR chemical modulators would help in the development of novel drugs for the pathologies in which this receptor is involved, such as glaucoma, inflammation, asthma and COPD, as well as several types of cancer. The medicinal chemistry developed around this receptor also responds to growing demand for pharmacological tools to study the (pato)physiology of the A 3 AR, and includes the development of agonists and antagonists, as well as radiolabeled versions of these [46,47]. Our initial efforts in this area resulted in the report of several potent and selective antagonists with a monocyclic pyrimidine scaffold (compound ISVY130, Figure 1) [24]. The pyrimidine core, being part of the heterocyclic moiety of the endogenous adenosine, is a recurrent substructural motif within biand tricyclic AR antagonists. The functionalized pyrimidine template could be decorated in a divergent and parallel fashion with one alkilamino substituent (L1) and two symmetric aromatic substitution points (L2 and L3, see Figure 5), by means of the well-established Suzuki–Miyaura crosscoupling reaction to a collection of commercially available boronic acids. The initial series of structures Figure 4. Synthesis and enantiospecific binding characterization of compound 16b as antagonist of the A 2B AR: ( A ) the Biginelli multicomponent reaction allows the assembly of the 3,4-dihydropyrimidine scaffold bearing different diversity points; ( B ) the chiral resolution identified the Sisomer as the one responsible of the biological affinity of the racemic mixture; and ( C ) the molecular modeling confirmed the binding mode previously hypothesized for isomer (S)- 16b . The water-density maps (black mesh) are calculated here by means of MD exploration (grid spacing 1 Å, occupancy > 80%), and are overlaid with the position of structural waters previously considered during the docking of this series (red spheres). The docking studies of these series of non-planar scaffolds were complemented by MD simulations, which allowed characterizing the role of water molecules in ligand binding. In a first stage, the selected complex of A 2B AR was equilibrated with our PyMemDyn protocol, as implemented in the GPCR-ModSim web server [ 44 , 45 ] (see Methods). From this analysis, we identified a water molecule that mediated the interaction between the NH in position 1 of the parent 3,4-dihydropyrimidin-2(1H)-ones and Glu169 EL2 , and a second water molecule in the region between the alkoxy substituent and the furan/thienyl ring, which is actually equivalent to a conserved water molecule observed in the crystal structures of A 2A AR with antagonist ZM241385 [ 23 ]. These pair of water molecules were maintained in successive docking runs for the expanded series of compounds, verifying their potential role in mediating this inter and intra molecular interactions [ 29 , 31 , 32 ]. We here provide further support for the existence of these structural water molecules, by means of an MD exploration of the solvent in the binding site (see methods). The idea is that the structural water molecules, which might favor the binding of ligands, is reflected in a higher occupancy during the MD trajectories, as we recently showed for the A 3 AR (see below and reference [ 30 ]). Figure 4shows the calculated water density maps at 80% occupancy, supporting the position of the two water molecules discussed above, mediating interactions between the ligand and E170 EL2 as well as the intramolecular interaction between the 2-furan and esther groups of the ligand. 2.4. A3AR The identification of A 3 AR chemical modulators would help in the development of novel drugs for the pathologies in which this receptor is involved, such as glaucoma, inflammation, asthma and COPD, as well as several types of cancer. The medicinal chemistry developed around this receptor also responds to growing demand for pharmacological tools to study the (pato)physiology of the A 3 AR, and includes the development of agonists and antagonists, as well as radiolabeled versions of these [46,47]. Our initial efforts in this area resulted in the report of several potent and selective antagonists with a monocyclic pyrimidine scaffold (compound ISVY130, Figure 1) [ 24 ]. The pyrimidine core, being part of the heterocyclic moiety of the endogenous adenosine, is a recurrent substructural motif within biand tricyclic AR antagonists. The functionalized pyrimidine template could be decorated in a divergent and parallel fashion with one alkilamino substituent (L1) and two symmetric aromatic Molecules 2017,22, 1945 9 of 17 substitution points (L2 and L3, see Figure 5), by means of the well-established Suzuki–Miyaura cross-coupling reaction to a collection of commercially available boronic acids. The initial series of structures combined structural simplicity with a low molecular weight (MW < 350), allowing for optimization during the subsequent hit-to-lead process [ 24 ]. This process was assessed by systematic docking on our A 3 AR homology model and the consensus binding mode obtained featured interactions with well-conserved residues in the ARs: the double hydrogen bond of the key residue Asn250 6.55 with both the exocyclic amido group and the closest nitrogen in the pyrimidine scaffold, and a π -stacking interaction with Phe168 EL2 complemented by hydrophobic interactions with Leu246 6.51 (Figure 5). In addition, this binding pose identified selectivity hotspots on the extracellular area around the L1 region: Leu 7.35 in TM7 (Thr/Met/Met on A 1 /A 2A /A 2B , respectively) Ile 6.58 (Val on the remaining ARs) and Val 5.30 in EL2 (otherwise Glu in the family). The binding mode and interactions identified were later supported by a crystal structure of the A 2A AR in complex with another monocyclic scaffold, i.e., 1,2,4-triazine, also bearing aryl substitutions [ 34 ]. The biological superposition obtained from the consensus docking was the starting point for a 3D-QSAR model with the second generation of GRid INdependent Descriptors (GRIND-2), used to rationalize and further optimize the three points of diversity of this scaffold. The first model was built on the initial series of 62 compounds, and showed excellent correlation with the experimental data (r 2 = 0.86) and good predictive power ( q2= 0.62 ). Based on this model, we investigated the role on ligand binding of different methoxyphenyl fragments of the L2/L3 diaryl substituents of the 4-amidopyrimidine scaffold. The increased selectivity of these series confirmed a second A 3 AR selective hotspot created by serines at positions 6.52 (His in the other ARs) and 5.42 (Asn in the other ARs), located deeper at the TM5-TM6 interface, which create a cavity that accommodates the metaand/or para-methoxy groups of the most potent and selective compounds [25]. The exploration of the pyrimidine scaffold included the characterization of the two isomeric series, i.e., 2-amidopyrimidines and 4-amidopyrimidines, differing in the position of the second nitrogen (N1) in the pyrimidine ring. The 3D-QSAR model suggested a role of this nitrogen in the affinity of these compounds, presumably by interaction with a potentially conserved structural water molecule within ARs linked to Ser271 7.42 , for which evidence in the hA 3 AR homology model was suggested by a GRID analysis with a water probe [ 24 ]. To further explore the complex effect of a “necessary nitrogen” [ 48 ] in the structure–activity (SAR) of these heterocycles, we designed a new series of bioisosteres derived from the 4-acetamidopyrimidine scaffold, i.e., lacking this nitrogen atom [ 30 ] (Figure 5). The most potent pyridines obtained displayed affinity values in the low nanomolar range, comparable to the best pyrimidines of the former series, maintaining the excellent selectivity profile toward the remaining ARs. The affinity data of the two series remained within the same nanomolar range, which confirmed our design of the bioisosteric replacement. However, a closer look into the data revealed that the nitrogen replacement by CH led to a small to moderate decrease in binding affinity in most compounds within the pyridine series. This would be indeed in line with our previous hypothesis that the second nitrogen would stabilize a water molecule (Figure 5) [ 24 ], which in analogy with the latest high resolution crystal structures of the A 2A AR, could be part of a water network in the deep cavity of the binding site. A comparative MD sampling of the solvent in the two situations (i.e., pyridine and pyrimidine bound states), revealed a specific increase in water densities around the second nitrogen of the pyrimidine, as compared to the A 3 AR-pyridine complex, thus supporting the idea of a water network stabilized by the binding of a molecule bearing the pyrimidine scaffold [ 30 ]. We therefore performed FEP simulations, which allowed to quantitatively estimating an understanding the origin of the relative binding affinities between analogous pairs of pyridine/pyrimidine compounds. The calculations, which essentially consisted of the N → CH perturbation in the ring, were in remarkable agreement with the experimental data, reproducing the small effect (8 fold loss in binding affinity) in the non-selective diphenyl substituted pair of compounds, as well as the moderate (30 fold) loss of affinity observed for the most potent compound in the parent pyrimidine series (4-metoxyphenyl derivative, ISVY130, and the corresponding pyridine ISVY177, see Figure 5) [ 30 ]. An extreme case was the complete Molecules 2017,22, 1945 16 of 17 39. Fredholm, B.B.; IJzerman, A.P.; Jacobson, K.A.; Linden, J.; Müller, C.E. International union of basic and clinical pharmacology. LXXXI. Nomenclature and classification of adenosine receptors—An update. Pharmacol. Rev. 2011,63, 1–34. [CrossRef] [PubMed] 40. Nguyen, A.T.N.; Baltos, J.A.; Thomas, T.; Nguyen, T.D.; Munoz, L.L.; Gregory, K.J.; White, P.J.; Sexton, P.M.; Christopoulos, A.; May, L.T. Extracellular loop 2 of the adenosine A 1 receptor has a key role in orthosteric ligand affinity and agonist efficacy. Mol. Pharmacol. 2016,90, 703–714. [CrossRef] [PubMed] 41. Dalpiaz, A.; Townsend-Nicholson, A.; Beukers, M.W.; Schofield, P.R.; Ijzerman, A.P. Thermodynamics of full agonist, partial agonist, and antagonist binding to wild-type and mutant adenosine A 1 receptors. Biochem. Pharmacol. 1998,56, 1437–1445. [CrossRef] 42. Palaniappan, K.K.; Gao, Z.G.; Ivanov, A.A.; Greaves, R.; Adachi, H.; Besada, P.; Hea, O.K.; Ae, Y.K.; Seung, A.C.; Lak, S.J.; et al. Probing the binding site of the A 1 adenosine receptor reengineered for orthogonal recognition by tailored nucleosides. Biochemistry 2007,46, 7437–7448. [CrossRef] [PubMed] 43. Townsend-Nicholson, A.; Schofield, P.R.A. Threonine residue in the seventh transmembrane domain of the human A 1 adenosine receptor mediates specific agonist binding. J. Biol. Chem. 1994 ,269, 2373–2376. [PubMed] 44. Gutiérrez-de-Terán, H.; Bello, X.; Rodríguez, D. Characterization of the dynamic events of GPCRs by automated computational simulations. Biochem. Soc. Trans. 2013,41, 205–212. [CrossRef] [PubMed] 45. Esguerra, M.; Siretskiy, A.; Bello, X.; Sallander, J.; Gutiérrez-de-Terán, H. GPCR-ModSim: A comprehensive web based solution for modeling G-protein coupled receptors. Nucleic Acids Res. 2016 ,44, W455–W462. [CrossRef] [PubMed] 46. Baraldi, P.G.; Preti, D.; Borea, P.A.; Varani, K. Medicinal chemistry of A 3 adenosine receptor modulators: Pharmacological activities and therapeutic implications. J. Med. Chem. 2012 ,55, 5676–5703. [CrossRef] [PubMed] 47. Ciancetta, A.; Jacobson, K. Structural probing and molecular modeling of the A 3 adenosine receptor: A focus on agonist binding. Molecules 2017,22, 449. [CrossRef] [PubMed] 48. Pennington, L.D.; Moustakas, D.T. The necessary nitrogen atom: A versatile high-impact design element for multiparameter optimization. J. Med. Chem. 2017,60, 3552–3579. [CrossRef] [PubMed] 49. Liu, W.; Chun, E.; Thompson, A.A.; Chubukov, P.; Xu, F.; Katritch, V.; Han, G.W.; Roth, C.B.; Heitman, L.H.; IJzerman, A.P.; et al. Structural basis for allosteric regulation of GPCRs by sodium ions. Science 2012 ,337, 232–236. [CrossRef] [PubMed] 50. Lebon, G.; Warne, T.; Edwards, P.C.; Bennett, K.; Langmead, C.J.; Leslie, A.G.W.; Tate, C.G. Agonist-bound adenosine A 2A receptor structures reveal common features of GPCR activation. Nature 2011 ,474, 521–525. [CrossRef] [PubMed] 51. Jaakola, V.-P.; Griffith, M.T.; Hanson, M.A.; Cherezov, V.; Chien, E.Y.T.; Lane, J.R.; Ijzerman, A.P.; Stevens, R.C. The 2.6 angstrom crystal structure of a human A 2A adenosine receptor bound to an antagonist. Science 2008 , 322, 1211–1217. [CrossRef] [PubMed] 52. Šali, A.; Blundell, T.L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993 , 234, 779–815. [CrossRef] [PubMed] 53. Fiser, A.; Šali, A. MODELLER: Generation and refinement of homology-based protein structure models. Methods Enzymol. 2003,374, 461–491. [PubMed] 54. Madhavi Sastry, G.; Adzhigirey, M.; Day, T.; Annabhimoju, R.; Sherman, W. Protein and ligand preparation: Parameters, protocols, and influence on virtual screening enrichments. J. Comput. Aided Mol. Des. 2013 ,27, 221–234. [CrossRef] [PubMed] 55. Schrödinger Release 2014-3: Maestro, version 9.9; Schrödinger, LLC: New York, NY, USA, 2014. 56. Verdonk, M.L.; Cole, J.C.; Hartshorn, M.J.; Murray, C.W.; Taylor, R.D. Improved protein-ligand docking using GOLD. Proteins Struct. Funct. Bioinform. 2003,52, 609–632. [CrossRef] [PubMed] 57. Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008 ,4, 435–447. [CrossRef] [PubMed] 58. Marelius, J.; Kolmodin, K.; Feierberg, I.; Åqvist, J.; Aqvist, J. Q: A molecular dynamics program for free energy calculations and empirical valence bond simulations in biomolecular systems. J. Mol. Graph. Model. 1998,16, 213–225. [CrossRef] 59. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 33–38. [CrossRef] Molecules 2017,22, 1945 17 of 17 60. King, G.; Warshel, A. A surface constrained all-atom solvent model for effective simulations of polar solutions. J. Chem. Phys. 1989,91, 3647. [CrossRef] 61. Lee, F.S.; Warshel, A. A local reaction field method for fast evaluation of long-range electrostatic interactions in molecular simulations. J. Chem. Phys. 1992,97, 3100. [CrossRef] 62. Ryckaert, J.-P.J.; Ciccotti, G.; Berendsen, H.J.H. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977,23, 327–341. [CrossRef] 63. Robertson, M.J.; Tirado-Rives, J.; Jorgensen, W.L. Improved peptide and protein torsional energetics with the OPLS-AA force field. J. Chem. Theory Comput. 2015,11, 3499–3509. [CrossRef] [PubMed] 64. Zwanzig, R.W. High-temperature equation of state by a perturbation method. I. nonpolar gases. J. Chem. Phys. 1954,22, 1420. [CrossRef] Sample Availability: Calculations of the affinity of compounds DPCPX for A1AR, and 10m and 10n on A2AAR are available from the authors. © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).