Stephanie
Federico
*a,
Enrico
Margiotta
b,
Silvia
Paoletta‡
b,
Sonja
Kachler
c,
Karl-Norbert
Klotz
c,
Kenneth A.
Jacobson
d,
Giorgia
Pastorin
e,
Stefano
Moro
b and
Giampiero
Spalluto
a
aDepartment of Chemical and Pharmaceutical Sciences, University of Trieste, Via Licio Giorgeri 1, 34127 Trieste, Italy. E-mail: sfederico@units.it
bMolecular Modeling Section (MMS), Dipartimento di Scienze del Farmaco, Università degli Studi di Padova, Via F. Marzolo 5, 35131 Padova, Italy
cInstitut für Pharmakologie und Toxicologie, Universität Würzburg, Versbacher Straße 9, 97078 Würzburg, Germany
dLaboratory of Bioorganic Chemistry, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, Maryland 20892-0810, USA
eDepartment of Pharmacy, National University of Singapore, 3 Science Drive 3, 117543 Singapore
First published on 18th February 2019
A series of adenosine receptor antagonists bearing a reactive linker was developed. Functionalization of these derivatives is useful to easily obtain multi-target ligands, receptor probes, drug delivery systems, and diagnostic or theranostic systems. The pyrazolo[4,3-e][1,2,4]triazolo[1,5-c]pyrimidine scaffold was chosen as a pharmacophore for the adenosine receptors. It was substituted at the 5 position with reactive linkers of different lengths. Then, these compounds were used to synthesise probes for the adenosine receptors by functionalization with a fluorescent moiety. Both series of compounds were evaluated for their binding at the four adenosine receptor subtypes. Different affinity and selectivity profiles were observed towards hA1, hA2A and hA3 adenosine receptors. In particular, fluorescent compounds behave as dual hA2A/hA3 ligands. Computational studies suggested different binding modes for developed compounds at the three receptors. Both molecular docking and supervised molecular dynamics (SuMD) simulations confirmed that the preferred binding mode at the single receptor was driven by the substitution present at the 5 position. Obtained results rationalized the compounds' binding profile at the adenosine receptors and pave the way for the development of more potent conjugable and conjugated ligands targeting these membrane receptors.
Adenosine receptors (ARs) are family A G protein-coupled receptors (GPCRs) involved in a number of physiological and pathological conditions that make them promising therapeutic targets for several diseases. In particular, four distinct AR subtypes have been identified so far: A1, A2A, A2B and A3 ARs.10 Each subtype has a specific distribution in the human body, which confers to their ligands particular potential therapeutic applications, such as in: pain conditions (A1 and A3 agonists),11 myocardial perfusion imaging (A2A agonists),12 cardiac and renal conditions (A1 antagonists),13 neurodegenerative diseases (A2A antagonists),14 cancer immunotherapy (A2A antagonists),14 inflammatory bowel disease (A2B antagonists)15 and glaucoma (A3 antagonists).16 It is clear that a pharmacophore, with tunable affinity towards AR subtypes might be useful in order to easily obtain theranostic agents, drug delivery systems or multivalent ligands for these receptors. The pyrazolo[4,3-e][1,2,4]triazolo[1,5-c]pyrimidine (PTP) is a well-known scaffold (Fig. 1) from which many potent AR antagonists have been derived. Among these, preladenant (1) is an A2AAR antagonist under clinical development to treat neurodegeneration and other pathological conditions.17 As depicted in Fig. 1, the PTP scaffold could be simultaneously substituted at different positions, which determines a desired specific affinity and selectivity profiles for the targets. In addition, we have extensive experience on this scaffold.18–21 For these reasons, in this work we have developed conjugable AR antagonists bearing a PTP scaffold (2–19). In order to understand the structural requirements crucial to retain the affinity and selectivity towards ARs, we have investigated these antagonists both from a pharmacological and computational perspective. In a previous study, these compounds were functionalized with a fluorophore leading to fluorescent receptor probes (20–25), which we have reported here, with their full characterization, as an example of application of our strategy.22,23
Fig. 1 General structure of the pyrazolo[4,3-e][1,2,4]triazolo[1,5-c]pyrimidine (PTP) nucleus and structure of the hA2AAR antagonist preladenant (1). |
In order to investigate the effect on the AR affinity, we decided to design nine homologous derivatives bearing diamino linkers of different lengths, from a hydrazino linker to longer diaminopolyglycolic chains. These chains are versatile because one amino group could be reacted with the AR ligand and the other with the detection moiety. In fact, there are several groups that could react with amines, such as carboxylic acids or their activated derivatives, isocyanates and isothiocyanates.
Compd | X | R | hA1a (Ki nM) | hA2Ab (Ki nM) | hA2Bc (Ki nM) | hA3d (Ki nM) |
---|---|---|---|---|---|---|
a Displacement of specific [3H]-CCPA binding at hA1ARs expressed in CHO cells (n = 3–6). b Displacement of specific [3H]-NECA binding at hA2AARs expressed in CHO cells (n = 3–6). c Inhibition of the NECA-stimulated adenylyl cyclase in CHO cells (n = 3–6) expressing hA2BARs. d Displacement of specific [3H]-HEMADO binding at hA3ARs expressed in CHO cells (n = 3–6). Data are expressed as geometric means with 95% confidence limits. | ||||||
2, WS86 | — | Boc | 7820 (6510–9380) | 3900 (3630–4200) | >10000 | 173 (124–241) |
3 | –CH2CH2– | Boc | 2550 (1740–3760) | 624 (467–834) | >10000 | 553 (470–651) |
4, WS88 | –CH2CH2CH2– | Boc | 1770 (1430–2200) | 195 (129–295) | >3000 | 32.1 (21–49.1) |
5 | –CH2CH2CH2CH2– | Boc | 1750 (1310–2340) | 313 (256–383) | 11500 (8590–15300) | 84.4 (79.5–89.7) |
6 | –CH2CH2CH2CH2CH2– | Boc | 2580 (1690–3940) | 346 (299–401) | 3000 | 129 (98.1–171) |
7, WS91 | –CH2CH2OCH2CH2OCH2CH2– | Boc | 12200 (10700–13800) | 1260 (995–1590) | >10000 | 585 (438–781) |
8 | –CH2CH2CH2O(CH2CH2O)2CH2CH2CH2– | Boc | 8480 (5960–12060) | 674 (568–799) | >10000 | 264 (136–515) |
9 | — | Boc | 4070 (3470–4770) | 2170 (1980–2370) | ∼10000 | 856 (696–1053) |
10 | — | Boc | >100000 | 21400 (18900–24300) | >10000 | 737 (505–1077) |
11 | — | H × HCl | 7100 (4870–10300) | 233 (128–422) | >10000 | 174 (126–241) |
12 | –CH2CH2– | H × HCl | 2500 (1480–4230) | 756 (454–1260) | >10000 | 2400 (1938–2971) |
13 | –CH2CH2CH2– | H × HCl | 1440 (1290–1610) | 253 (186–344) | >10000 | 3784 (3031–4723) |
14 | –CH2CH2CH2CH2– | H × HCl | 1220 (1080–1370) | 169 (108–263) | >3000 | 4205 (3715–4760) |
15, WS98 | –CH2CH2CH2CH2CH2– | H × HCl | 778 (608–994) | 52.6 (45.1–61.4) | >3000 | 2246 (2139–2359) |
16 | –CH2CH2OCH2CH2OCH2CH2– | H × HCl | 5110 (4830–5400) | 314 (256–387) | >10000 | 2662 (2461–2880) |
17 | –CH2CH2CH2O(CH2CH2O)2CH2CH2CH2– | H × HCl | 1660 (1350–2050) | 348 (269–450) | >10000 | 312 (264–368) |
18 | — | H × HCl | 5520 (3530–8620) | 682 (610–762) | >10000 | 2914 (1714–4955) |
19 | — | H × HCl | >100000 | 12900 (10400–15900) | >10000 | 8646 (7293–10250) |
In particular, in the Boc-protected series (2–10) we observed low affinity towards hA1AR: long or bulky linkers (7–10) were not tolerated at this receptor, while simple alkyl chains lead to Ki values in the range of 1.7–2.6 μM (i.e.7, hA1AR Ki = 12200 nM vs.4, hA1AR Ki = 1770 nM). Better results were obtained at the hA2AAR, where the affinity follows the same trend observed at the hA1AR. Nevertheless, most of the compounds displayed affinities in the submicromolar range, with a propyl linker providing the best Ki value (4, hA2AAR Ki = 195 nM). On the other hand, Boc-protected derivatives were quite good A3AR antagonists (2–10, hA3AR Ki = 32.1–856 nM). In particular, a clear hA3AR affinity enhancement was observed comparing derivatives 3 and 4 (WS88), which bear an ethyl and a propyl chain, respectively (3, hA3AR Ki = 553 nM vs.4, hA3AR Ki = 32.1 nM). Elongation of the linker caused a progressive decrease of hA3AR affinity (i.e.5, hA3AR Ki = 84.4 nM, 7, hA3AR Ki = 585 nM) as also observed when constrained in a piperazine ring (10, hA3AR Ki = 737 nM). With a Ki value of 32.1 nM, compound 4 was the most potent derivative towards the hA3AR, and it showed 55- and 6-fold selectivity against hA1 and hA2A ARs, respectively. Thus, it was the analogue that displayed the best balance between affinity and selectivity. The most selective compound for the hA3AR was the piperazine derivative 10, which displayed a >136- and 29-fold higher affinity value at the hA3AR, with respect to hA1 and hA2A ARs, respectively, but it showed an affinity in the high nanomolar range towards hA3AR (Ki = 737 nM). The other compounds of this subseries displayed a common selectivity profile at the ARs that was very low against hA2AAR (hA2A/hA3 = 1.1–3.7) and higher towards hA1AR (hA1/hA3 = 20–55). A different behaviour was observed when the shorter hydrazine linker was present at the 5 position of PTP. In fact, first of all, with compound 2 (WS86) we gained hA3AR affinity with respect to the ethyl compound 3 (3, hA3AR Ki = 553 nM vs.2, hA3AR Ki = 173 nM). Moreover, this compound exhibited one of the best selectivity profiles against both hA2AAR (22.5-fold) and hA1AR (45-fold).
However, when the Boc group was removed (leading to the free amino series 11–19) affinity at the hA1AR was not particularly affected (e.g.9, hA1AR Ki = 4070 nM vs.18, hA1AR Ki = 5520 nM). Conversely, at the hA3AR we observed a pronounced detrimental effect (e.g.9, hA3AR Ki = 856 nM vs.18, hA3AR Ki = 2914 nM) and a simultaneous increase of hA2AAR affinity (e.g.9, hA2AAR Ki = 2170 nM vs.18, hA2AAR Ki = 682 nM). As a consequence, in general, free amino compounds could be considered hA2AAR antagonists. Among these, an exception was seen with compounds 11 and 17, which showed similar Ki values towards hA2A and hA3 ARs (11, hA2AAR Ki = 233 nM, hA3AR Ki = 174 nM; 17, hA2AAR Ki = 348 nM, hA3AR Ki = 312 nM). Thus, they could be considered as dual hA2A–hA3 AR antagonists. Interestingly, it seems that for hA3AR, with the hydrazine linker compounds (2, 11) the presence of the Boc group did not influence the affinity at this receptor (2, hA3AR Ki = 173 nM vs.11, hA3AR Ki = 174 nM), while it had a strong effect on the hA2AAR (2, hA2AAR Ki = 3900 nM vs.11, hA2AAR Ki = 233 nM). Notably, observing data for compounds bearing oxygen-containing chains (16, 17), the longer {[(propoxy)-ethoxy]-ethoxy}-propyl linker conferred higher hA3AR affinity than the diethoxyethyl linker (17, hA3AR Ki = 312 nM vs.16, hA3AR Ki = 2662 nM). A similar behaviour was observed, although to a less extent, in Boc-protected compounds (8, hA3AR Ki = 264 nM vs.7, hA3AR Ki = 585 nM). Concerning the structure–activity relationship (SAR) at the hA2AAR, an affinity enhancement was observed upon chain elongation from 2 to 5 methylene units that decreased when longer or bulkier linkers were present at the 5 position. Compound 15 displayed the best hA2AAR affinity of the series, with a Ki value of 52.6 nM and a 15- and 43-fold selectivity against hA1 and hA3 ARs, respectively.
Longer linkers, which should extend into the extracellular part of the receptor, gave especially low affinity compounds and variable selectivity profiles. This observation suggested that the linker had a very crucial role in the recognition process. In order to better understand this behaviour, we performed computational studies on the ligand–receptor recognition process.
Scheme 2 Synthesis of fluorescent probes 20–25.22 Reagents and conditions: a) FITC, TEA, DMF, rt, 48 h. |
Compd | X | hA1b (Ki nM) | hA2Ac (Ki nM) | hA3d (Ki nM) |
---|---|---|---|---|
a All data are from ref. 22. b Displacement of specific [3H]-RPIA binding at hA1ARs expressed in CHO cells (n = 3–5). c Displacement of specific [3H]-CGS21680 binding at hA2AARs expressed in HEK293 cells (n = 3–5). d Displacement of specific [125I]I-AB-MECA binding at hA3AR expressed in CHO cells (n = 3–5). Duplicate determination. Data are expressed as geometric means with 95% confidence limits. | ||||
20 | — | >10000 | 2480 (1640–3320) | 2390 (1070–3710) |
21 | –CH2CH2CH2– | >10000 | 460 (340–580) | 720 (580–860) |
22 | –CH2CH2CH2CH2– | >10000 | 1750 (810–2690) | 5660 (4440–6880) |
23 | –CH2CH2CH2CH2CH25– | 9500d | 920 (700–1140) | 780d |
24 | –CH2CH2OCH2CH2OCH2CH2– | >10000 | >10000 | 9390 (5930–12850) |
25 | –CH2CH2CH2O(CH2CH2O)2CH2CH2CH2– | >10000 | 6010 (5170–6850) | 7090 (6350–7830) |
The low affinity of these compounds highlights again the crucial role of the linker in the recognition process, suggesting that the PTP core alone is not sufficient to provide the desired affinity and selectivity at the hARs. The introduction between the PTP core and the linker of a substituent able to achieve a specific affinity for one of the receptors could enhance the affinity and selectivity of the conjugated compounds.
As previously reported, also in this case, regardless of the nature of the receptor subtype considered, docking simulations suggested, for all ligands belonging to the PTP class, three different and alternative binding modes (BM) in their orthosteric binding sites.20 With a specific reference to the position assumed by the furyl moiety, the binding modes differ for the orientation of the ring along the main axis of the binding cavity (see details in Fig. 3): in particular, BM1 and BM3 orient the furyl ring down compared to the extracellular environment (for BM1, alternative up orientation is also retrieved), while BM2 orients it towards the TM2. Interestingly, the three binding modes guarantee the stabilizing interaction with Asn (6.55) and with other hydrophobic residues, such as Phe (EL2). The relative BM sampling frequency strongly depends critically on the nature of the substituent at the 5-position of the PTP scaffold, suggesting that the specific substituents there can modulate the receptor potency and/or receptor subtype selectivity. Considering this specific PTP series, docking studies suggested that BM2 is the most frequently sampled for any receptor, even if BM3 is its most relevant alternative on the hA3AR. In Fig. 4–6 are shown the interaction energy fingerprint (IEF) maps of both electrostatic and hydrophobic contributions, calculated for each sampled pose on all three receptor subtypes.
Interestingly, regardless of the nature of the receptor subtype considered, almost any analyzed PTP analogue can be accommodated in the orthosteric site by the canonical binding pattern, already annotated for other AR antagonists: these include the hydrogen bonds and the aromatic hydrophobic π–π interactions with Asn (6.55) and Phe (EL2), respectively. Results can be related to the structural and sequence diversities between receptors. Indeed, favorable electrostatics are detected along the EL1 and TM2 for ligands docked to the hA2AAR, while, for the hA1AR, some repulsion is achieved at the TM3, i.e. non-conserved residues play a major role.
To analyze in more detail the interactive scheme of some representative ligands belonging to this class of antagonists, per-residue interaction energy profiles (electrostatic, van der Waals and hydrophobic contributions) were calculated for compounds 4 (WS88), 7 (WS91), 15 (WS98) and 23, and shown in Fig. 7–10. Conserved and non-conserved residues were selected for each receptor subtype to better appreciate their contributions in the stabilizing ligand–receptor binding state. Moreover, docking poses are also collected in Videos S1–S3 (ESI†), together with electrostatic and hydrophobic contributions of peculiar residues.
Fig. 7 Interaction energy histograms of compound 4 (WS88) for conserved and non-conserved residues of any receptor subtype. |
Fig. 8 Interaction energy histograms of compound 7 (WS91) for conserved and non-conserved residues of any receptor subtype. |
Fig. 9 Interaction energy histograms of compound 15 (WS98) for conserved and non-conserved residues of any receptor subtype. |
Fig. 10 Interaction energy histograms of compound 23 for conserved and non-conserved residues of any receptor subtype. |
The first analyzed antagonist was derivative 4, the most potent compound of the series towards the hA3AR (Ki = 32.1 nM), despite its modest selectivity (hA1/hA3 = 55.3 and hA2A/hA3 = 6.1). Interestingly, derivative 4 preferably recognizes the orthosteric binding site of the hA1AR using the interactive scheme BM1up, the hA2A subtype using the interactive scheme BM3, and the hA3 using the interactive scheme BM1down. As shown in Fig. 7, these three different binding modes are characterized by three different interaction energy profiles and the analysis shows that the most favorable interactive scheme is for subtype hA3 engaging strong hydrophobic contacts with the conserved Phe168 (EL2), with the non-conserved Leu90 (3.32), and an electrostatically favorable interaction with Glu19 (1.39), without showing relevant repulsive van der Waals contributions. Conversely, an unfavorable van der Waals term is detected when derivative 4 is docked to both hA1 and hA2A ARs, mediated by Asn254 and Asn253, respectively (6.55).
The second analyzed antagonist is derivative 7, a moderately selective compound for the hA3AR, especially when compared to the hA1AR (hA1/hA3 = 135.7 and hA2A/hA3 = 2.1). This antagonist preferably recognizes the orthosteric binding site of the hA1AR using the interactive scheme BM2, and the hA2A and hA3 using the interactive scheme BM3. Interestingly, unless for Glu172 (EL2) and Val87 (3.32), the analysis of the interaction energy profiles indicates a very modest stabilizing interaction network compared to the hA3AR. At the hA2AAR, the ligand interacts unfavorably with His250 (6.52). The most favorable profile is observed at the hA3AR.
Derivative 15 is the highest affinity compound of the series at the hA2AAR (Ki = 52.6 nM). Nevertheless, it shows affinity constants in high nanomolar (Ki = 778 nM) and micromolar ranges (Ki = 2246 nM), for the hA1 and hA3 ARs, respectively. Derivative 15 preferably recognizes the orthosteric binding site of the hA1 and hA2AARs using the interactive scheme BM2, but interactive scheme BM3 is used at the hA3AR. For this ligand, at the hA1AR, several stabilizing interactions were detectable, in particular with Phe171 (EL2), Glu172 (EL2) and Asn254 (6.55), as shown in Fig. 9. This stabilizing network of interactions is also present at the hA2A AR, even if the intensity of each interaction energy term is lower than that measured at the hA1AR. To appropriately interpret the lower affinity determined experimentally against the hA3 subtype, it is possible to invoke the crucial repulsive interactions detected for Val169 (EL2) and Phe168 (EL2), together with highly unfavorable van der Waals contributions by Asn250 (6.55). Involvement of these residues is likely to be critical for affinity.
The final analyzed antagonist is derivative 23, the most representative compound among the fluorescent ligands of the series, and the most selective towards hA3AR. As expected, BM2 is the only interactive scheme sampled for each receptor subtypes, allowing a better comparison (Fig. 10). As already described, this specific binding mode allows us to orient the substituent that carries the fluorescent probe towards the extracellular environment. Unfortunately, in this case, considering the relevant volumetric properties of these ligands and the fact that their receptor recognition involves not only the orthosteric binding sites but also relevant portions of the extracellular loops (in particular, the EL2), the analysis of the interaction energy profiles is not very helpful to interpret the experimental selectivity profile. However, it is emphasized that all the most crucial interactions with the conserved Phe (EL2) and Asn (6.55) are preserved. Moreover, the amino acid sequences and the architectures of EL2 play an important role in controlling both the selectivity profiles and ligand affinities of the fluorescent analogues against all receptor subtypes.
Finally, to support or negate the hypothesis suggested by docking studies concerning the presence of different and alternative ligand binding modes (BM) inside the receptor orthosteric binding site, we decided to carry out a supervised molecular dynamics (SuMD) study with the aim to explore the possible ligand–receptor recognition pathways (from an unbound state to the orthosteric bound state) and to verify if the final states reached by supervised trajectories converge to the most populated poses obtained by molecular docking. For this specific scope, we decided to use the smaller PTP analogue of this series, derivative 2 (WS86), which showed a moderate affinity against the hA3AR (Ki = 173 nM) and a reasonable selective profile (hA1/hA3 = 45.2 and hA2A/hA3 = 22.5). Docking simulations suggested that derivative 2 preferably recognizes the orthosteric binding site of the hA1AR using the interactive scheme BM1up; the hA2AAR used the interactive scheme BM3, whereas the hA3 used two distinct interactive schemes BM1down and BM2. As already reported for all previously analyzed antagonists, these four different binding modes are characterized by four different interaction energy profiles and the analysis shows that the most favorable interactive scheme for subtype hA3AR is the one in which the ligand preferably adopts the interactive scheme BM1down, establishing the already described stabilizing interaction with the conserved Phe168 (EL2) and Asn250 (6.55) (Fig. 11).
Fig. 11 Interaction energy histograms of compound 2 (WS86) for conserved and non-conserved residues of any receptor subtype. BM1down and BM2 are reported for the hA3AR subtype. |
Considering the appreciable hA3AR affinity, all SuMD simulations have been performed using this specific receptor subtype. Five independent supervised simulations have been carried out starting from an unbound state in which the ligand has been randomly positioned 40 Å far away from the center of the orthosteric hA3AR binding cavity.
As expected, and already reported in previous studies,26 each SuMD trajectory describes a specific ligand–receptor recognition pathway (see Fig. 12, Panel A) and allows the ligand to sit in different regions of the orthosteric binding site (bound state) and adopt different binding poses (see Fig. 12, Panel B). Even if the different bound states obtained by SuMD simulations are not perfectly superimposable to those suggested by docking experiments, we emphasize that the possibility to have a multiple binding mode for a single ligand is clearly shown in the results of the SuMD trajectory analysis. Particularly interesting is the result extrapolated from Replica 5, in which the observed bound state is close to the one described as BM2, as shown in Video S4.† As previously reported, also for derivative 2, the different SuMD trajectories described the presence of several meta-binding sites (energetically stable binding sites spatially distinct from the orthosteric site) with an energy profile (Fig. 12, Panel A) that smoothly decreases from an average value 〈E1〉 of −30 kcal mol−1 (for distances >10 Å) to an average value 〈E2〉 of −50 kcal mol−1 (〈E2〉 ≈ 2 〈E1〉), in a distance window between 2 and 10 Å. This behavior of the interaction energy profile is typical of a ligand–receptor recognition process characterized by a reasonably high kon value.
As already anticipated, different bound states obtained by SuMD simulations are not perfectly superimposable to the poses suggested by docking experiments, even if in Replica 5 the SuMD bound state is close to the one described as BM2. Starting from the Replica 5 final bound state, a non-biased molecular dynamics simulation has been performed to analyse the time-evolution (stability) of this specific bound state. Surprisingly, after 250 ns of simulation, the ligand flips its original position assuming a new configuration very similar to the BM3 detected in docking simulations, as shown in Fig. 13, and more energetically favored. This evidence indicates that the accommodation of the furyl ring down, compared to the extracellular side, is the most favored one.
The analysis of both supervised and non-biased MD trajectories showed that the adaptive positioning of the ligand inside the orthosteric binding site is strictly related to the simultaneous conformational variation of the side chains of the two amino acids Phe168 (EL2) and Trp243 (6.48) or, in other words: the final ligand pose depends on the reshaping of the binding cavity (induce fit) caused by the side chains of the two amino acids Phe168 (EL2) and Trp243 (6.48). Moreover, considering that all SuMD simulations have been carried out with a sodium ion and its first hydration shell in the hA3AR structure, it is interesting to note that both the sodium ion and its coordinated water molecules are also involved in this transition (Video S4†).
In conclusion, the combination of molecular docking and dynamics (supervised and non-biased) simulations can be considered an interesting computational strategy to improve the descriptive capacity of the interaction modality of this new series of AR antagonists. In particular, SuMD can facilitate the interpretation at a molecular level of the possible ligand–receptor recognition pathways bridging the stability of the canonical orthosteric binding (thermodynamic binding propensity) with the ease of accessing the orthostatic site (kinetic binding propensity).
The hA2A and hA1 AR coordinates were retrieved from the Protein Data Bank39 using the crystal structures coded as 3PWH40 and 5UEN,41 respectively. Given the lack of any crystallographic structures for hA3AR, a homology model was used. The model was built using 5UEN as a template, as reported in our recently submitted work,42 after a docking/structural based model assessment. According to our recently published work,25 all molecular docking studies have been carried out with a sodium ion and its first hydration shell in all adenosine receptor structures.
The Ballesteros–Weinstein43 numbering system is used sometimes to indicate conserved residues. IE electrostatic and van der Waals energy contributions to the binding energy were calculated by MOE, together with per residue electrostatic and hydrophobic interactions. Per-residue information was reported in the so-called “Interaction Energy Fingerprints”: they consist of heat maps reporting the strength of the interaction of each residue (x-axis) and each ligand (y-axis) according to a colorimetric scale going from blue to red for negative to positive values in the case of electrostatic contributions and from white to dark green for low to high values for hydrophobic contributions. Interactions of most relevant residues were also reported by histograms, whose height is proportional to the strength of the interaction. Plots were generated using Gnuplot 4.6.44 2D chemical structures in figures were drawn with Marvin Sketch.45 Molecular graphics were performed with the UCSF Chimera package,46 except Fig. 12B, which was obtained using the VMD program (version 1.9.3).47 The in-house MMsDocking video maker tool was exploited to produce videos showing the docking poses, per residue hydrophobic and electrostatic contributions for selected residues, experimental binding data, and scoring values. Representations of docking poses were produced using the UCSF Chimera package, 2D depictions were constructed by the cheminformatics toolkit RDKit,48 and the heat maps were obtained by Gnuplot 4.6; in the end, videos were mounted using MEncoder.
Ligand 3 force field parameters for MD simulations were initially retrieved from the Paramchem web service and then deeply optimized in concordance with CGenFF, at the MP2/6-31G* level of the theory49 by using Gaussian 09 (ref. 50) and RESP partial charges. Systems were embedded in a 1-palmitoyl-2oleyl-sn-glycerol-3-phospho-choline (POPC) lipid bilayer, according to the pre-orientation provided by the Orientations Proteins in Membrane (OPM) database51 and by using the VMD membrane builder plugin. Lipids within 0.4 Å from the protein were removed and TIP3P52 model water molecules were added to solvate the system by means of Solvate1.0.53 System charge neutrality was reached by adding Na+/Cl− counterions to a final concentration of 0.154 M. Equilibration was performed through a three-step procedure. In the first one, 1500 conjugate-gradient minimization steps were applied in order to reduce the clashes between protein and lipids. Then, a 5 ns long MD simulation was performed in the NPT ensemble, with a positional constraint of 1 kcal mol−1 Å−2 on ligand, protein, and lipid phosphorus atoms. During the second stage, 10 ns MD simulation in the NPT ensemble was performed constraining all the protein and ligand atoms but leaving the POPC residues free to diffuse in the bilayer. In the last equilibration stage, positional constraints were applied only to the ligand and protein backbone alpha carbons for further 5 ns MD simulation. All the MD simulations were performed using: (1) an integration time step of 2 fs; (2) the Berendsen barostat54 in order to maintain the system pressure at 1 atm; (3) the Langevin thermostat55 to maintain the temperature at 310 K with a low dumping of 1 ps−1; (4) the M-SHAKE algorithm56 to constrain the bond lengths involving hydrogen atoms; and (5) a long-range cutoff of 10 Å. According to the SuMD approach,57,58 the timescale needed to reproduce complete intermolecular complex formations is in the range of nanoseconds, instead of hundreds of nanoseconds or microseconds usually necessary with unsupervised MD. Sampling is gained by applying a tabu-like algorithm to monitor the distance between the centers of mass of the ligand and the binding site, during short unbiased MD simulations. SuMD considers the ligand atoms and the atoms of user-defined protein residues to monitor the distance between the centers of mass of the binder and the binding site. A series of 600 ps unbiased MD simulations are performed and after each simulation, the distance points collected at regular time intervals are fitted into a linear function. If the resulting slope is negative the next simulation step starts from the last set of coordinates and velocities produced, otherwise the simulation is restarted by randomly assigning the atomic velocities. Short simulations are perpetuated under the supervision until the distance between the ligand and the binding site goes below 5 Å, then, the supervision is disabled and a classical MD simulation is performed. For the orthosteric center of mass, we considered hA3AR residues N250, F168, H272, and S247.
NAMD energy Plugin 1.4 (ref. 59) was used to calculate ligand–protein interaction energies.
Footnotes |
† Electronic supplementary information (ESI) available: Videos concerning docking poses, electrostatic and hydrophobic contributions of peculiar residues and results of SuMD for synthesised compounds. See DOI: 10.1039/c9md00014c |
‡ Current affiliation: Sygnature Discovery, Nottingham, UK. |
This journal is © The Royal Society of Chemistry 2019 |