Understanding Ligand Binding Selectivity in a Prototypical GPCR Family

Adenosine receptors are involved in many pathological conditions and are thus promising drug targets. However, developing drugs that target this GPCR subfamily is a challenging task. A number of drug candidates fail due to lack of selectivity which results in unwanted side effects. The extensive structural similarity of adenosine receptors complicates the design of selective ligands. The problem of selective targeting is a general concern in GPCRs, and in this respect adenosine receptors are a prototypical example. Here we use enhanced sampling simulations to decipher the determinants of selectivity of ligands in A2a and A1 adenosine receptors. Our model shows how small differences in the binding pocket and in the water network around the ligand can be leveraged to achieve selectivity.


■ INTRODUCTION
Adenosine receptors are Class A G-protein-coupled receptors (GPCRs) that are widely expressed in neurons, immune system cells, and vascular system. 1 The adenosine signaling has long been shown to regulate cytoprotective and immune response in tissues. 2 Of the four subtypes, A 1 R and A 2a R have received most attention from the scientific community, mainly due to their therapeutic potential for the treatment of reperfusion injury and neuropathic pain (A 1 R), cancer, attention deficit hyperactivity disorder (ADHD), and Parkinson's Disease (PD) (A 2a R). 3,4 Notwithstanding significant efforts, the success in developing molecules with high potency and selectivity has been limited. 5 The proteins have a sequence identity of 37%, 6 with limited structural differences and remarkably similar orthosteric binding pockets 7−9 (Figure 1), as all key pocket-lining residues are conserved with the exception of the position 7.35 (Ballesteros-Weinstein numbering scheme, 10 where the first number refers to the transmembrane helix and the second to the position relative to the most conserved residue of the helix, which is arbitrarily set to 50). The lack of differences in the pockets of the two receptors poses a real challenge for the design of selective ligands. 5 Thus, understanding how ZM241385, an A 2a R-selective inverse agonist, achieves its selectivity is of great importance. To this aim we have recently solved and analyzed a new crystal structure of the complex. 11 However, the highly dynamical nature of the receptors is not fully captured by the crystal structures, 11−14 leaving a number of questions unanswered. For instance, while ZM241385 shows greatly reduced affinity for A 1 R, its dehydroxy derivative LUF5452 15 binds with the same affinity to the two proteins ( Table 1). The two ligands are very similar (Figure 2), and the crystallographic binding mode of ZM241385 to A 2a R shows that the only part of the ligand differing to LUF5452 does not form significant direct interactions with the receptor. 11 The challenges of explaining the differential binding affinities of ligands based on static poses hints at an important role of the binding mechanisms and perhaps different flexibility of the receptors. Thus, to fully understand the determinants of selectivity, the crystal structures need to be complemented by a method that is able to reconstruct the binding process in atomistic detail and provide information on the different dynamical properties of the receptors. To this end, here we use molecular dynamics simulations (MD) combined with enhanced sampling algorithms, that have been successfully used to model ligand binding in GPCRs. 11,16,17 ■ RESULTS AND DISCUSSION Interaction of the Ligands with the Binding Sites. To elucidate the binding pose and the role of the interfacial water, 1 μs-long unbiased MD were run for human A 1 R and A 2a R in complex with either ZM241385 or LUF5452 starting with the ligand bound to the binding site. An important difference in the binding sites of the two receptors is constituted by the 7.35 position, T270 7.35 in A 1 R and M270 7.35 in A 2a R. Its role in determining the selectivity of ligands with bulky substituents has been proposed in previous studies. 9 Indeed, in our simulations we observe a poorer stacking of the ligand cores in A 1 R due to the smaller side chain of the threonine residue with respect to the methionine of A 2a R. Thus, while in A 2a R the ligands remained close to the initial pose (with two relatively similar conformations), in A 1 R a significant reorientation of the (4-hydroxy)phenyl-ethyl tail was observed.
Of the two conformations adopted by the ligands in A 2a R, one is close to the starting crystallographic structure (PDB 5IU4 11 ), with the (4-hydroxy)phenyl-ethyl tail pointing toward the extracellular side of the protein, and one is very similar to the crystallographic structure of the thermostabilized receptor (PDB 3PWH), 13 with the tail pointing toward TM1. The various conformations can be monitored using the RMSD of the common heavy atoms of the two ligands ( Figure 3). While ZM241385 preferentially adopts the upright (5IU4-like) position in the simulations, LUF5452 tends to explore 3PWH-like conformations. In this conformation, ZM241385 is able to form polar interactions with Y9 1.35 , E13 1.39 , and A63 2.61 . The interconversion of the poses has significant effect on the hydration of the tail, with partial desolvation happening for 3PWH-like states and allowing LUF5452 to minimize unfavorable interactions between the solvent and the phenyl ring, as seen in Figure 4. A hydrogen bond network involving the ligands and water molecules stabilizes the molecules in the pocket of A 2a R ( Figure 3, conformations a 1 and b 1 ), as is also evident from crystallographic structures. 11,12 Although the cavity is wide enough for observing much movement of solvent molecules, there is a tendency to adopt ordered patterns of interaction, that are disrupted in 3PWH-like states. The desolvation of the phenyl ring leads LUF5452 to frequently adopt such conformations and leads to a poor arrangement of water molecules, contributing to the lower affinity for the receptor with respect to ZM241385.
In A 1 R, the ligands are much more mobile. ZM241385 was observed rotating the tail to a 3PWH-like conformation ( Figure 3, conformation c 2 ) and temporarily interacting with an hydrophobic cavity of the extracellular vestibule (ECV) formed by ECL2 and the ends of the transmembrane helices (TM) 2 and 3 (See Figure S1a). LUF5452 instead undergoes a more significant conformational rearrangement, rotating its core away from the upright orientation and projecting the tail toward the same pocket without reverting to the starting pose ( Figure 3, conformation d 2 ). In this conformation the hydrogen bond between the exocyclic primary amine of the ligand and N254 6.55 is broken to allow the phenyl ring of the molecule to interact with the pocket. The distance between the amine nitrogen of LUF5452 and the carbonyl oxygen of N254 6.55 increases from 3 to 6 Å. This results in a significant rearrangement of the interfacial water molecules (see Figure  S2). Similar to A 2a R, adopting a 3PWH-like conformation results in a partial desolvation (0.30 > ligand RMSD to 5IU4 > 0.45 nm in Figure 4). Moreover, the wide volume of the pocket of A 1 R allows for the presence of a large number of solvent molecules. While in this receptor the interactions of the tail of ZM241385 with the solvent are stabilized by the hydroxy group, LUF5452 is instead able to bury the phenyl ring in the hydrophobic ECV pocket, leading to more effective desolvation than in A 2a R ( Figure 4B).
Binding Poses of the Ligands. To better quantify the relative free energy of the different poses, as well as to compute the full free energy landscape associated with the binding of the ligands, we run parallel tempering metadynamics simulations (PT-metaD) 19,20 on the four complexes. For A 2a R, the chosen collective variables (CVs) were the projection of the vector connecting the binding pocket and the ligand onto the Z-axis   15 and Guo et al. 21 and was converted with the relation ΔG = −RT ln(K a ). Superscripts refer to the source organism: (a) rat or (b) human. Data in kcal/mol.

Journal of Chemical Information and Modeling
Article (Z-projection) and onto the XY-plane (XY-projection) parallel to the membrane. In A 1 R, extensive tests showed that including the distance between the salt bridge-forming residues, E172 and K265 speeds the convergence of the free energy. Thus, for ZM241385, we used the distance and the Zprojection, while for LUF5452 all three CVs were used (see the Supporting Information). The deepest free energy minimum for the two ligands corresponds to the crystallographic binding pose of ZM241385 11,12 in A 2a R and closely resemble it in A 1 R ( Figure  5). As expected, the presence of the same key residues in the pocket allows the molecules to form analogous interactions and therefore adopt similar binding poses with respect to A 2a R. The binding free energy associated with the minima in all four systems agrees well with experimental values (Table 1).
LUF5452 Interacts with the Accessory Hydrophobic Site. As observed in the unbiased MD simulations, LUF5452 also shows an equally stable secondary minimum in which it makes significant contacts with the hydrophobic ECV pocket of A 1 R ( Figure 5). The interaction allows the ligand to rotate away from the expected pose, breaking the hydrogen bond between the exocyclic amine and the side chain of N254 6.55 ( Figure S1a). The region is lined by L65 2.60 , A66 2.61 , I69 2.64 , V83 3.28 , A84 3.29 , C169, E170, and F171. In this conformation, the phenyl group of the ligand is nearly parallel to TM2 and is positioned between the side chains of I69 2.64 and F171. The ring of F171 is involved in a stacking interaction with the core of the ligand and a T-shaped (quadrupolar) one with its phenyl group. The presence of the aromatic ring of LUF5452 in the pocket displaces a number of water molecules ( Figure S1a).
Analysis of the pocket with the GRID software 22 confirms the highly hydrophobic nature of its surface, with overlap between the hotspot identified with the C1 (lipophilic, carbon sp2 probe) and the phenyl ring of LUF5452 ( Figure  S3).
As also observed in the unbiased MD, the hydrophobic nature of the region makes the interaction with the solvent unfavorable and displacing the "unhappy waters" increases the strength of the interaction with the ligand stabilizing the secondary minimum. The importance of exploiting the presence of "unhappy waters" has been previously shown in other GPCRs. 23,24 Salt Bridge Influences Binding. In both proteins a salt bridge is present at the entrance of the binding site, hindering the binding and unbinding of the ligands. While in A 2a R the interaction is formed between the residues E169 and H264; in A 1 R a lysine residue (K265) forms a stronger contact with E172. The different nature of the salt bridge in the two receptors results in a much higher stability of the E172-K265 contact in A 1 R with respect to its counterpart in A 2a R ( Figure  S4a) and was reflected in the need to bias the distance between the side chains in our metadynamics simulations of A 1 R. We also find that the strong salt bridge in A 1 R hinders the binding and unbinding of the ligand ( Figure 5). On the contrary, the Figure 3. Analysis of the conformations adopted by ZM241385 and LUF5452 in the binding pocket of A 2a R and A 1 R over the unbiased MD simulation. Probability over the space defined by the RMSD of the common heavy atoms of the ligands to the conformation adopted in PDB 5IU4 or 3PWH, after alignment to the backbone of the proteins. TM7 not shown for clarity.

Journal of Chemical Information and Modeling
Article weaker salt bridge in A 2a R has a smaller influence on the binding dynamics.
The restriction of the cross-sectional area of the entrance of the binding site determined by the helical turn of ECL2 hosting E169 (A 2a R) or E172 (A 1 R) and ECL3 forces the ligand to approach the salt bridge at a right angle between the common core plane and the bridge axis. Generally, the molecules initially interact with the residues with either the furane ring or with the bicyclic region of the cores. In A 2a R, stacking of the core and the side chain of H264 is also observed ( Figure 5).
The importance of the salt bridge of A 2a R for the binding mechanism is supported by mutagenesis and structural data 11,25 and was previously suggested for A 1 R in the literature. 26 Dominant Binding Pathway of the Ligands. The free energy landscapes of the four sets of simulations also indicate a significant role of the extracellular vestibule in driving the binding process. ECL2, in particular, is often the first point of contact between the ligand and the receptors, as reported for β-adrenergic receptors. 27 The conformation of ECL2 is different in the two receptors. While a short helical turn is present in both proteins toward the orthosteric binding pocket, hosting F168 and E169 (A 2a R) and the equivalent F171 and E172 in A 1 R, a longer helical segment extending further into the bulk solvent differs in length and conformation (see Figure  1). In the binding path, we observe frequent interactions between the ligands and the loops, especially with A151, A155, A158, and K173 in A 1 R and K150 or K153 in A 2a R. The contact between ECL2 and the ligands likely helps preorient the molecules toward the binding site, forcing the cores of the molecules to lay on the loop surface, and funnelling the ligands to the salt bridge at the entrance of the site (Figure 6), similar

Journal of Chemical Information and Modeling
Article to what was observed by Dror et al. 27 for β-adrenergic receptors.

■ CONCLUSIONS
Our results suggest that a combination of three structural differences of A 1 R and A 2a R concur to the observed selectivity of ZM241385, while leading to similar binding affinities of LUF5452. First, the salt bridge at the opening of the orthosteric binding site of A 1 R and A 2a R significantly influences the dynamical nature of the interaction between the receptor and ligands. Second, albeit ZM241385 and LUF5452 possess a similar binding mode, the presence of a hydrophobic pocket in the binding site of A 1 R allows LUF5452 to adopt an alternative pose, displacing weakly interacting water molecules in the region. Third, the desolvation of the 4hydroxyphenyl and phenyl groups of the two ligands leads to the exploration of alternative conformations at the expense of the water network that stabilizes them in the main binding pose. While in A 2a R this results in poor stabilization of LUF5452, the hydrophobic pocket of A 1 R allows an effective desolvation of the group. Taken together, our findings confirm how the differential behavior of GPCR ligands arises from a complex interplay of small but relevant structural and dynamical differences and from different solvation patterns. The new insights on the molecular determinants of ligand selectivity between A 1 R and A 2a R provide a clear indication for the structure-based rational design of selective drugs for these receptors and possibly for other GPCRs.

■ METHODS
The crystal structures of the human A 1 and A 2a receptors were downloaded from the Protein Data Bank (A 1 R: 5N2S, 7 A 2a R: 5IU4 11 ). After correcting mutations and removing the apocytochrome b562 (bRIL) inserted in the ICL3 loop or at the N-terminus, missing residues were added with MOD-ELLER 9.19 28 and the proteins were embedded in a preequilibrated POPC membrane. The membrane was then aligned to the XY plane, resulting in the main axis of the GPCRs and the Z axis being close to parallel.
The complexes were solvated with TIP3P water, and chloride ions were added to balance the net charge. Ligands were parametrized with the generalized AMBER force field (GAFF) 29 and charges were calculated using Gaussian09 30 with a 6-31G* basis set at the Hartree−Fock level. Protein, water, and ion parameters were generated with AMBER14SB force field, 31,32 and the phospholipid topology was downloaded from LipidBook. 33 The simulations were run using GRO-MACS 5.1.4 34 with the PLUMED 2.3.1 plugin 35 in the NPT ensemble.
The metadynamics simulations were run using a parallel tempering scheme in the 300−310 K temperature range. During the production metadynamics runs, Gaussian hills were deposited every 2 ps in the well-tempered scheme 36 with a bias factor of 15. The Gaussian width was set to 0.1 nm for the Zprojection and XY-projection and to 0.033 nm for the salt bridge distance. In the simulations where two CVs were biased, the initial height was 1 kJ/mol, whereas in that biasing three CVs it was set to 1.5 kJ/mol.
The exploration of the bulk water region by the ligand was restrained with the use of a funnel-like restraint in order to aid the convergence of the simulations. 17,37 ■ ASSOCIATED CONTENT

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.9b00298.
Comparison of rat and human receptor sequences; details of the setup of the molecular dynamics and

Journal of Chemical Information and Modeling
Article parallel-tempering metadynamics simulations; and further information on the hydrophobic pocket of A 1 R, the role of water in the orthosteric pocket of A 1 R, and on the results of the metadynamics simulations (PDF)