New Insights into Key Determinants for Adenosine 1 Receptor Antagonists Selectivity Using Supervised Molecular Dynamics Simulations

Abstract

Adenosine receptors (ARs), like many otherGprotein-coupledreceptors (GPCRs), are targets of primary interest indrug design. However, one of the main limits for the development of drugs for this class of GPCRs is the complex selectivity profile usually displayed by ligands. Numerous efforts have been madefor clarifying the selectivity of ARs, leading to the development of many ligand-based models. The structure of the AR subtype A₁ (A₁AR) has been recently solved, providing important structural insights. In the present work, we rationalized the selectivity profile of two selective A₁AR and A_2AAR antagonists, investigating their recognition trajectories obtained by Supervised Molecular Dynamics from an unbound state and monitoring the role of the water molecules in the binding site.

1. Introduction

Adenosine receptors (ARs) are class A G protein-coupled receptors (GPCRs) that bind the endogenous agonist adenosine.ARs are composed of foursubtypes: A₁, A_2A, A_2B, A₃. While A₁ and A₃ARs (which share 49% of a sequence identity) are preferentially coupled to Gα_iproteinsand therefore inhibit the adenylatecyclase, A_2A and A_2BARs (sharing 59% of a sequence identity) stimulatethis enzyme, as being coupled to Gα_s proteins [1]. Several ARs antagonists are in clinical trials for various diseases. With regards to A_2AAR, is tradefylline has been recently approved for Parkinson’sdisease (NCT02610231) [2], PBF-509 is in phase I/II trials for non-small cell lung cancer (NCT02403193), and CPI-144 is in phase I trials for various cancer types (NTC02655822). PBF-680, on the other hand, is the only A₁ AR antagonist inclinical phase II, for the treatment of Asthma (NTC02635945) [3].

One of the difficulties during the development of ARs agonists and antagonists as therapeutic agents is the poor selectivity between differentreceptorssubtypes [4]. For this reason, many efforts have been madeto elucidate the molecular basis of ARs ligands selectivity, and several structure-activity relationship (SAR) models have been developed for selective ligandsof all the four subtypes [1,2,5,6,7,8,9,10,11,12,13,14].

With the increasing availability of structural information (mainly from mutagenesis, X-ray, and cryo-EM approaches [15]) in the last few years, light has been shed on the origin of selectivity onARs.Recently, the A₁AR inactive (PDB, Protein Data Bank, code 5UEN [16] and 5N2S [17]) and active (PDB code 6D9H [18]) structures have been solved. Interestingly, in [17], Cooke and colleagues obtained the X-ray crystal structure of A₁ and A_2A Ars, in a complex with the same xanthine ligand PSB36, providing insight about the selectivity. There are several structural differences between A₁AR and A_2AAR. The second extracellular loop (ECL2), in particular, is more folded in A₁AR and orients perpendicularly to the plane of the membrane, while in A_2A ARit forms a longer helix, which is parallel to the lipid bilayer. This difference is probably due to the presence of two disulfide bondsuniquely present in A₂AR. Indeed, the bond between Cys71 and Cys159 anchorsECL2 to ECL1,while the bond between Cys74 and Cys146 tethers TM3 to ECL2. The class A conserved disulfide bonds between Cys80 and Cys169 is present in both the two subtypes. These differences in the disulfide bonds likely contribute to the outward movement of the top of transmembrane helix 2 (TM) in the inactive A₁AR (Figure 1).

Further divergence involves TM7, which is shifted outward compared to A_2AAR, due to the shorter ECL3, andTM6 slightly shifted inward in A₁AR. These rearrangements, in turn, affect the orthosteric site of A₁AR, which is wider than A_2AAR. Interestingly, the key residues in the orthosteric site of the two receptors are conserved and drive the same binding mode of the antagonist PSB36 (Figure 1). More precisely, the xanthine scaffold forms two hydrogen bonds with Asn254 (A₁AR, Asn253 in A_2A AR) and a π- π stacking with Phe171 (A₁AR, Phe168 in A_2AAR). Nevertheless, Asn254 (A₁AR) is located in the binding site deeper than Asn253 in A_2AAR, and the xanthine ligand is consequently positioned deeply in the orthosteric site (Figure 1).

Despite the huge help provided from high-resolution structural biology techniques, certain selectivity profiles cannot be only rationalized by the mere coordinates of bound state or “final” state. Ligand recognition is an articulated mechanism in which many variables may play a relevant role and over the last few years, there has been rising attention in the understanding of binding kinetics at GPCRs and its determinant role to successfully target this class of proteins [19].

In the present study, we used supervised molecular dynamics (SuMD) simulations toshed light on the molecular basis of the selectivity of three different ligands to A₁AR and A_2AAR, not only considering the bound states, but also the possible different recognition mechanism preceding the final orthosteric site and the role of the solvent. We focus our attention on three antagonists: the A_2AAR selective antagonists Z48 (Ki 16.9 nM in A_2AAR and 1345.7 nM in A₁AR) [20]; the A₁AR selective antagonist LC4 (Ki 16,800 nM in A_2AAR and 89 nM in A₁AR) [21]; and the nonselective antagonist caffeine (Figure 2). SuMD [22,23] is a molecular dynamics (MD) approach that allows for the study of molecular recognition processes in a fully atomistic way, in the nanosecond timescale, without introducing any energetic biases.

2. Materials and Methods

2.1. System Setup

The crystal structures of the two receptors were retrieved from PDB (the PDB code is 5N2S for A₁-AR and 5NM4 for A_2A-AR). Systems preparation was performed using a Molecular Operating Environment (MOE)) suite (Chemical Computing Group ULC, Molecular Operating Environment (MOE), 2019.01. 1010 Sherbrooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2019) [24] for protein preparation (removal of crystallographic water molecules, ions, and other solvent molecules, selection of the highest occupancy for each residue, assignment of the correct protonation state at pH 7.4). Systems preparation for the molecular dynamics simulations was carried out using VMD [25]. The protein was explicitly solvated in a water box with the borders placed at a distance of 15 Å from any protein atom, the water model used was TIP3P [26]. The system charge was neutralized to a concentration of 0.154 M using Na⁺/Cl⁻. The lipid bilayer consisted of phosphatidylcholine (POPC) units.

The sodium ionwithin the TMD allosteric site of A_2A-AR was retained, andit was also placed by superposition in A₁-AR.

2.2. Equilibration of the System

All the simulations were performed with a CHARMM36 force field [27] and using ACEMD2 [28]. Ligands parameters were retrieved from Paramchem [29], a web interface for the assignment of parameters based on the CGenFF [30] force field.

The system energy was minimized in 1500 steps using the conjugate-gradient method, then the equilibrationof the system was done in four steps. The first one consistedof 5 ns of NPT simulation with harmonic positional constraints of 1 kcal mol⁻¹Å⁻² on each atom of the protein and the lipid bilayer. The second step consisted in 10 ns of NPT simulation with harmonic positional constraints of 1 kcal mol⁻¹Å⁻² only on eachprotein atom and on the phosphorus atom of the POPC units, the third step consisted in 5 ns of NPT simulation with harmonic positional constraints of 1 kcal mol⁻¹Å⁻² only on the alpha carbons of the protein, and the last step consisted in 50ns of NVT simulation without any constraints.

For the productive simulations, the temperature wasmaintained at 310 K using the Langevin thermostat, with a low dumping of 1 ps⁻¹. The pressure was set at 1 atm using the Berendsen barostat [31]. The particle-mesh Ewald (PME) method was used to calculate the electrostatic interactions with a 1 Å grid [32]. A 9.0 Å cutoff was applied for long-term interactions.The M-SHAKE algorithm was applied to constrain the bond lengths involving hydrogen atoms.

At the end of the equilibration, several parameterswere calculated to assess the stability of the system:theroot mean square deviation (RMSD) of the alpha carbons of the protein, the root mean square fluctuation (RMSF) of each protein residue, the volume of the cell (which should tend to a plateau in the NPT ensemble), and the area per lipid (APL) for each membrane layer (calculated using GridMAT-MD [33]). We also computed the volume of the orthostericsiteduring the equilibration using POVME [34]. Figure S1 (A_2A-AR) and Figure S2 (A₁-AR) report the analysis performed during the equilibration of the system. In both the two systems, the protein reached a stable conformation (RMSD of the protein Cα (panel A, Figures S1 and S2) stably below 2 Å for A_2A-AR and below 3 Å for A₁-AR). The volume of the orthosteric site reached a plateau in both cases (Panel B, Figures S1 and S2). The most flexible parts of the protein, as expected, are the loops. Indeed, the RMSF of these regions is higher than the TMs (panels C and D, Figures S1 and S2).

As shown in Figure S5, the orthosteric site of A₁-AR appears to be deeper than A_2A-AR, due to a cleft between TM5 and TM6. The APL and the volumetric analysis are reported in Figure S3 (A₁-AR) and Figure S4 (A_2A-AR). For both systems, the cell volume reached stable values.

2.3. Supervised Molecular Dynamics Simulations

SuMD [22,23] is a molecular dynamics (MD) approach that allows for the investigation of molecular recognition processes in a fully atomistic way, in the nanosecond timescale, without introducing any energetic bias (Figure 3). Ligands were placed 35 Åaway from the protein. EachSuMD step was set to 600 ps. During each SuMD step, the distance between the center of mass of the binding site (defined by a series of residues) and the center of mass of the ligand is monitored. These dataare thenfitted, and if the slope of the interpolating linear function is negative, then the coordinates and the velocities are used for the successive time window, otherwise, the last time window is simulated again reassigning the velocities (this reassignment of the velocities is intrinsic in the use of Langevin thermostat). If the condition fails 30 consecutivetimes, then the simulation is stopped. Otherwise, the algorithm continues until the distance between the two centers of mass is below the threshold of 5 Å; at this point, the supervision is turned off and 30 short classical MD simulations are performed, switching on the supervision if the distance between the centroids becomesgreater than 5 Å. At the end of the SuMD process, the trajectoryis prolonged for 25 ns of classical MD simulation.

Twenty simulations were performed for each of the six systems (Z48/A₁AR, Z48/A_2AAR, LC4/A₁AR, LC4/A_2AAR, Caffeine/A₁AR, Caffeine/A_2AAR). Only the simulationsthat sampled the ligand reaching the orthosteric site and interacting with the classic fingerprint of these class of ligands (Figure S6) were here reported (e.g., one replica for each system, excepted LC4/A₁AR for which two replicas were analyzed). For each system, the reasonsfor failure are similar. In most cases, the ligands interact strongly with the residues of the ECLs, and do not reach the binding site. In a few cases, the ligands only partially reach the orthosteric site. Finally, in other rare cases, the ligands get stuck between the protein and the membrane.

2.4. Trajectories Analysis

The SuMD trajectories were analyzed using anin-house python tool (described in [35]) that provides information on the geometry and the energetic of the system. The output consists of a per-residue analysis of the electrostatic and van Der Waals contributions to the protein-ligand interactions; a representation of the distance between the center of mass of the ligand and the center of mass of the binding site as a function of time; a global energetic evaluation of the system as a function of the aforementioned distance. This analysis allowscomparing the energetic profile of two systemsboth in a general way and, through the per-residue analysis (i.e., it is possible to evaluate which ligand has better interaction with some specific amino acids of interestover the time).

Three replicas of 50 ns were performed on the apo form of A₁-AR and A_2A-AR, after the equilibration stage described before. These three replicas were merged and analyzed by AquaMMapS [36].

3. Results

3.1. SuMD Binding of the A₁AR Nonselective Antagonist Caffeine

According to the SuMD simulations (Figures S7 and S8; Video S1 for A₁AR and Video S2 for A_A2AAR in the Supplementary Information), the nonselective antagonist caffeineestablishes intermediate interactions with the extracellular vestibule of A₁AR and A_2AAR, before reaching the orthosteric site. The most stable bound configurations sampled on A₁AR and A_2AAR differed for theorientation of the xanthine ring (Figure S9). On A_2AAR, caffeine pointed the N7-methyl toward Asn253, while in A₁AR it was rotated by 180°, with the N3 methylinthe proximity of Asn254. Notably, both these two conformations have been experimentally observed in X-ray crystal structuresofA_2AAR [17]. Fluctuations of the ligand in both orthosteric sites and the transient nature of the interactions are easily depicted in the energy interaction landscapes (Panels A,S7,S8) in which the points are particularly scattered over the distance between the centers of mass of the ligand and the orthosteric site.

3.2. SuMD Binding of the A_2AAR Selective Antagonist Z48

Figures S10 and S11 reportthe analysis of Z48/A₁ARand Z48/A_2AAR, including 25 ns of classic MD simulations performed at the end of each SuMD simulation (Videos S3 and S4 in the Supplementary Information).

For what concerns the SuMD trajectory of Z48 on A₁AR, notable electrostatic repulsion took place between the ligand and Lys168, Lys173, and Lys265, before the ligand reached the orthosteric site, as clearly observable in the per-residue electrostatic interaction energy plot (Panel C,S10). These residues are positioned on the ECLs (Lys168 and Lys173 on ECL2, andLys265 on ECL3). As a result, in A₁AR, Z48 did not adopt the binding fingerprint of the ARs antagonists (e.g., only one hydrogen bond with Asn254 out of two was formed). Moreover, this binding mode wasunstable over the 25 ns of classic MD simulations.Interestingly, the terminal amine group of the ligand strongly interacts with Glu170 (Leu167 in A_2AAR). Figure 4 reports the binding modes of Z48 at the end of the SuMD simulation on A₁AR (the binding mode at the end of the 25 ns of classic MD simulation is shown in Figure S12) and A_2A-AR (at the end of the 25 ns of classic MD simulation), respectively.

Moving to the binding simulations of Z48 toA_2AAR, two SuMDreplicas led to the classic binding mode of the ARs antagonists. The π-π stacking with Phe168 was present, along with the hydrogen bonds with Asn253 and with Glu169. Both binding modeswere stable during 25 ns of classic MD simulations (Figure S11, Panel C,D, Figure 4). No significantprotein-ligand electrostatic repulsionswereobserved. The narrow funnel-like interaction energy profile of the ligand in the orthostatic also suggested a good complementarity of the ligand in the pocket, and a rapid reaching of a stable bound state.

3.3. SuMD Binding of the A₁AR Selective Antagonist LC4

With regards to LC4, two SuMDreplicas led to different binding modes of LC4 into the A₁AR orthosteric site. During Replica 1, the ligand formed a complex characterized by the xanthine scaffold positioned in the orthosteric binding site, and the N8-substituent pointed outward the receptor (Figure 5). A hydrogen bond with Asn254 and hydrophobic contacts with Phe171 occurred (Figure 6A). Interestingly, the LC4 oxygen atom in the N8 substituent interacted with and stabilized water molecules in the proximity of the Phe171 backbone, a hydrated spot on ECL3, according to the AquaMMapS analysis (Figure 6A).

In the case of A₂AR, the xanthine scaffold reached the orthostericsite rotated by 180° compared to the binding mode adopted in A₁AR (Figure 6B) and with an unfavorable geometry for hydrogen bonding with Asn253. In both of these two binding modes (Figures S13 and S14), LC4 did not interactwithconserved glutamate Glu172 (A₁AR) or Glu169 (A_2AAR).

During SuMDReplica 2 on A₁A, LC4 experienced a two-step binding (Figure 6). First, the antagonist entered the orthosteric site pointing the methylphenyl group into a hydrophobic pocket located between TM2 and TM3. This cryptic pocket is not present in A₂AR in light of bulkier residues and a higher degree of packing between the helixes (Figure 6B). From this metastable configuration, LC4 moved deeper into the orthosteric site and engagedAsn254 and Glu172 in hydrogen bonds (Figure 6). In this binding mode, the ligand inserted the N8-substituent inside a further cryptic pocket between TM5 and TM6 (Figure 6C,D), which is delimited by the “toggle switch” residue W247 [15], and the residues Ile95 and Phe253, being part of the conserved class A structural motif PIF [37]. Notably, this hydrophobic sub-pocket was occupied by likely “unhappy” water molecules during simulations of the apo-A₁R (Figure 6A).

The video of these three simulations can be found in Supplementary Materials Video S5 and Video S6 for the two simulations of LC4/A₁AR, and Video S7 for the simulation of LC4/A_A2AAR.

4. Discussion

Here we present results from SuMD simulations performed to shed light on the selectivity displayed by LC4 and Z48, two antagonists of A₁AR and A_2AAR, respectively. The nonselective antagonist caffeine was also dynamically docked to the two ARs subtype.Caffeine, which is in a weak binder(micromolar range [38]) of all the ARs, experienced more than one binding mode, in line with our previous simulations [22] and experimental observation [17].

SuMD binding of the selective A_2AAR antagonistZ48 on A₁-AR and A_2AAR suggested different interaction patterns along the pathways. Z48 experienced unfavorable electrostatic interactions between positively charged A₁AR residues Lys168, Lys173 (ECL2), and Lys265 (ECL3), and the ligand charged amine on the N8-substituent. Such transitory states did not take place during binding to A_2AAR, as no significant electrostatic repulsions were computed at the extracellular vestibule. Interestingly, A_2AAR bears Ala265 instead of Lys265 on the ECL3, while Lys168 and Lys173 (ECL2) are farther from the binding site, compared to A₁AR. This is consistent with mutagenesis studies that demonstrated the importance of these lysine residues forthe binding of several A₁AR ligands [39,40,41]. On A_2AAR, Z48 reachedthe orthosteric site producing the classic interactions fingerprint of the ARs antagonists. On A₁AR, on the other hand, the ligand sampled a different binding mode. Moreover, the A₁AR residue Glu170 (Val167 in A_2AAR) strongly interacted with the charged terminal amine of the ligand, stabilizing this alternative binding mode. Z48 was proposed to bind to A₂AR overcoming low enthalpy transition state(s) [40]. From this standpoint, the unfavorable electrostatic interactions with ECL2 could implicate a slower binding toA₁AR, and therefore a kinetic selectivity for A₂AR.

SuMD simulations of the antagonist LC4 proposed two possible binding modes that could drive selectivity. It is possible that the ligand binds A₁AR and A_2AAR with the same conformation, but differently interacting with water molecules nearby ECL3. From this standpoint, and in analogy with [41], we propose these “happy” water molecules contribute to the ligand stabilization in A₁AR (Figure S16), but not in A₂AR (Figure S17). An alternative and unique LC4 binding mechanism was sampled only on A₁AR, with a metastable state before the final complex formation. Along this pathway, two cryptic hydrophobic pockets (between TM2 and TM3 and between TM5 and TM6) allowed the N8-substituent of the ligand to correctly orient first, and then engage key residues for the receptor activation (the “toggle switch” Trp247, Ile95, and Phe253, which are part of the conserved class A motif PIF). Notably, the cryptic pocket between TM2 and TM3 has recently been proposed as a determinant for A₁AR selectivity displayed by the triazolotriazine antagonist LUF5452 [42,43].

5. Conclusions

Understanding the selectivity of GPCRs ligands is an important task in drug design. This study supports the emerging idea that selectivity is driven by a plethora of phenomena, other than the protein-ligand interactions in the bound state. Receptor-ligands recognitions are multistep events modulated by intermediate interaction along withthe (un)binding paths. This picture may be further complicated by the presence of stable water molecules, whichcan have a tremendous impact on stabilizing or destabilizing an orthosteric complex. To consider different aspects that may affect the selectivity on A₁AR and A_2AAR, we used SuMD simulations to investigate the recognition of thee different antagonists. Overall, our results suggest that kinetic selectivity may favor the binding of Z48 to A_2AAR over LC4. A different scenario was observed for A₁AR, the recognition trajectories highlighted the key role of water molecules in the binding mode of LC4, which is favored by two hidden sub pockets within A₁AR.