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, A1R and A2aR have received most attention from the scientific community, mainly due to their therapeutic potential for the treatment of reperfusion injury and neuropathic pain (A1R), cancer, attention deficit hyperactivity disorder (ADHD), and Parkinson’s Disease (PD) (A2aR).[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 ), 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]

Figure 1

ZM241385 bound to the binding sites of (a) A2aR and (b) A1R in the enhanced-sampling simulations. (c) The salt bridge hindering the pocket entrance is shown as spheres. A zoomed view of the ligands in the binding sites and the key residues is presented in the boxes.

Thus, understanding how ZM241385, an A2aR-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 A1R, its dehydroxy derivative LUF5452[15] binds with the same affinity to the two proteins (Table ). The two ligands are very similar (Figure ), and the crystallographic binding mode of ZM241385 to A2aR shows that the only part of the ligand differing to LUF5452 does not form significant direct interactions with the receptor.[11]

Table 1

Experimental and Calculated Binding Free Energy of ZM241385 and LUF5452 to A1R and A2aRa

	A1R		A2aR
	ΔGexp	ΔGcalc	ΔGexp	ΔGcalc
ZM241385	–9.05a,–7.81b	–8.23 ± 0.73	–12.01a,–12.90b	–11.87 ± 0.43
LUF5452	–11.19a	–10.34 ± 0.27	–11.35a	–10.57 ± 0.53

Experimental affinity data was obtained from de Zwart et al.[15] and Guo et al.[21] and was converted with the relation ΔG = −RT ln(Ka). Superscripts refer to the source organism: (a) rat or (b) human. Data in kcal/mol.

Figure 2

Ligands in this study, ZM241385 and LUF5452.

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 A1R and A2aR 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, T2707.35 in A1R and M2707.35 in A2aR. 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 A1R due to the smaller side chain of the threonine residue with respect to the methionine of A2aR. Thus, while in A2aR the ligands remained close to the initial pose (with two relatively similar conformations), in A1R a significant reorientation of the (4-hydroxy)phenyl-ethyl tail was observed.

Of the two conformations adopted by the ligands in A2aR, 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 ). 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 Y91.35, E131.39, and A632.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 . A hydrogen bond network involving the ligands and water molecules stabilizes the molecules in the pocket of A2aR (Figure , conformations a and b), 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.

Figure 3

Analysis of the conformations adopted by ZM241385 and LUF5452 in the binding pocket of A2aR and A1R 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.

Figure 4

Solvation of the 4-hydroxyphenyl or phenyl groups of the ligands over the unbiased MD simulations, represented by the number of water oxygen atoms within 4 Å from the moieties’ heavy atoms. The plots illustrate the average solvation and its standard deviation as a function of the RMSD of the ligand to an upright, 5IU4-like conformation (see Figure ). Conformations A and B show the displacement of water molecules in the hydrophobic ECV pocket of A1R by LUF5452.

In A1R, the ligands are much more mobile. ZM241385 was observed rotating the tail to a 3PWH-like conformation (Figure , conformation c) 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 , conformation d). In this conformation the hydrogen bond between the exocyclic primary amine of the ligand and N2546.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 N2546.55 increases from 3 to 6 Å. This results in a significant rearrangement of the interfacial water molecules (see Figure S2). Similar to A2aR, adopting a 3PWH-like conformation results in a partial desolvation (0.30 > ligand RMSD to 5IU4 > 0.45 nm in Figure ). Moreover, the wide volume of the pocket of A1R 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 A2aR (Figure B).

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 A2aR, the chosen collective variables (CVs) were the projection of the vector connecting the binding pocket and the ligand onto the Z-axis (Z-projection) and onto the XY-plane (XY-projection) parallel to the membrane. In A1R, 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 Z-projection, 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 A2aR and closely resemble it in A1R (Figure ). 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 A2aR. The binding free energy associated with the minima in all four systems agrees well with experimental values (Table ).

Figure 5

Projection of the free energy landscapes onto Z-projection and XY-projection by reweighting.[18] Boxes: Representative conformations of the ligands bound to the main binding site of the receptors, LUF5452 interacting with the hydrophobic ECV pocket in A1R, and the two molecules breaking the salt bridge at the entrance of the orthosteric site of the proteins.

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 A1R (Figure ). 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 N2546.55 (Figure S1a). The region is lined by L652.60, A662.61, I692.64, V833.28, A843.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 I692.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 A2aR the interaction is formed between the residues E169 and H264; in A1R 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 A1R with respect to its counterpart in A2aR (Figure S4a) and was reflected in the need to bias the distance between the side chains in our metadynamics simulations of A1R. We also find that the strong salt bridge in A1R hinders the binding and unbinding of the ligand (Figure ). On the contrary, the weaker salt bridge in A2aR 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 (A2aR) or E172 (A1R) 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 A2aR, stacking of the core and the side chain of H264 is also observed (Figure ).

The importance of the salt bridge of A2aR for the binding mechanism is supported by mutagenesis and structural data[11,25] and was previously suggested for A1R 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 (A2aR) and the equivalent F171 and E172 in A1R, a longer helical segment extending further into the bulk solvent differs in length and conformation (see Figure ). In the binding path, we observe frequent interactions between the ligands and the loops, especially with A151, A155, A158, and K173 in A1R and K150 or K153 in A2aR. 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 ), similar to what was observed by Dror et al.[27] for β-adrenergic receptors.

Figure 6

Representative binding or unbinding paths of the ligands, extracted from the PT-metaD simulations. The ligands need to break the salt bridge at the entrance of the binding site in order to unbind. Extensive interactions with ECL2 are found before complete solvation of the ligands in the fully unbound states.

Conclusions

Our results suggest that a combination of three structural differences of A1R and A2aR 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 A1R and A2aR 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 A1R allows LUF5452 to adopt an alternative pose, displacing weakly interacting water molecules in the region. Third, the desolvation of the 4-hydroxyphenyl 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 A2aR this results in poor stabilization of LUF5452, the hydrophobic pocket of A1R 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 A1R and A2aR 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 A1 and A2a receptors were downloaded from the Protein Data Bank (A1R: 5N2S,[7] A2aR: 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 MODELLER 9.19[28] and the proteins were embedded in a pre-equilibrated 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 GROMACS 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 Z-projection 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]