Binding of the Antagonist Caffeine to the Human Adenosine Receptor hA2AR in Nearly Physiological Conditions.

Lipid composition may significantly affect membrane proteins function, yet its impact on the protein structural determinants is not well understood. Here we present a comparative molecular dynamics (MD) study of the human adenosine receptor type 2A (hA(2A)R) in complex with caffeine--a system of high neuro-pharmacological relevance--within different membrane types. These are POPC, mixed POPC/POPE and cholesterol-rich membranes. 0.8-μs MD simulations unambiguously show that the helical folding of the amphipathic helix 8 depends on membrane contents. Most importantly, the distinct cholesterol binding into the cleft between helix 1 and 2 stabilizes a specific caffeine-binding pose against others visited during the simulation. Hence, cholesterol presence (~33%-50% in synaptic membrane in central nervous system), often neglected in X-ray determination of membrane proteins, affects the population of the ligand binding poses. We conclude that including a correct description of neuronal membranes may be very important for computer-aided design of ligands targeting hA(2A)R and possibly other GPCRs.

Introduction

Increasingly emerging experiments point to an important role of membrane lipid composition in structure/function relationships of G-protein coupled receptors (GPCRs)—the largest membrane protein family in mammals [1-4]. For instance, in rhodopsin, adding (1-palmitoyl-2-oleoylphosphatidylethanolamine) POPE lipids into POPC (1-palmitoyl-2-oleoylphosphatidylcholine) lipid bilayers affect the equilibrium between two sub-states of the rhodopsin functional cycle [5,6]. Another example is given by the class A GPCR serotonin1A receptor: its binding to the agonist 8-OH-DPAT decreases with cholesterol concentration [7]. This issue is crucial for pharmacological applications, as more than one quarter of FDA-approved drugs target GPCRs [8].

Molecular Dynamics (MD) simulations are being instrumental to assess the effect of native membrane environment on conformational properties of GPCRs. Indeed, MD studies of class A GPCR, beta-2 adrenergic receptor, using four different membrane types with or without cholesterol, have suggested that the stability of the receptor “ionic lock” [9] varies with different lipid compositions [4]. In addition, 1.6 μs MD simulations of rhodopsin embedded in (1-stearoyl-2-docosahexaenoyl-phosphatidylcholine) SDPC/SDPE (1-stearoyl-2 docosahexaenoyl-phosphatidylethanolamine) lipid bilayer, in the presence of cholesterol, have lead to the conclusion that specific cholesterol-rhodopsin binding modulates the TM1-TM2-TM7 helices/helix 8 interactions, essential for the receptor’s activation [10]. Finally, MD simulations indicates that the stabilization of the amphipathic helix 8 (H8) of the class C GPCR mGluR2 receptor increases with cholesterol concentration and that such stabilization depends also on membrane thickness [3].

Here we use MD simulations to elucidate the effect of membrane composition on an antagonist-bound neuronal GPCR. This is the class A GPCR human adenosine receptor type 2A (hA2AR) [11]—highly localized in the so-called “striatum” of the brain [12]—in complex with the antagonist caffeine (CFF, chemical formula in S1 Fig). CFF binding to hA2AR may lead to neuroprotection [13-17]. It prevents apoptotic cell death in a Parkinson’s' model [18]. The effect on passing from artificial lipid bilayers to conditions near to real neuronal membranes, where cholesterol content varies from 33% to 50% [19], is investigated by performing 0.8 μs-long MD simulations on three systems, named in sequence as I, II and III (see Table 1). I is composed by pure POPC lipid bilayer (see S1 Fig), which is commonly used for MD studies on adenosine receptor [20-25]. II is a bilayer of equally mixed POPC and POPE lipids (S1 Fig). III resembles the synaptic membrane (42% POPC, 34% POPE and 25% of cholesterol molecules, S1 Fig), where hA2AR is expressed [26] and the binding to CFF exerts its beneficial neuroprotection effects [16]. Notably, III markedly differs from the artificial membrane mimics (detergent n-nonyl-β-D-glucopyranoside), where the CFF/hA2AR complex is embedded for X-ray structure determination [27].

Table 1

Composition of the three systems simulated here.

System	Lipid ratio	#POPC molecules	#POPE molecules	#Cholesterol molecules	#Water molecules	Ions	#Total atoms
I	POPC	495	-	-	26839	152 Na+; 165 Cl-	152,414
II	POPC/POPE (0.50:0.50)	247	249	-	27244	152 Na+; 165 Cl-	151,522
III	POPC/POPE /Cholesterol (0.42:0.34:0.24)	248	204	141	25703	148 Na+; 161 Cl-	151,834

Results and Discussion

Convergence Analysis

The degree of convergence of the 800 ns long MD simulations of systems I, II and III is here investigated by using the so-called ‘all-to-all RMSD analysis’ [28]. This assembles pairs of Cα atoms' RMSDs in matrices along an MD simulation. In the last 400 ns, the matrices of the three systems (S2 Fig) show a "leave-and-return pattern", a converged-alike feature according to [28], see supporting information. On longer time-scales, this feature is not observed. This suggests that the protein conformations of the three systems might have reached a fair convergence in the last 400 ns timescale here. Consistently with this fact, the Cα RMSDs oscillate around average values in such time scale (Panel A in S3 Fig). The properties presented here are therefore calculated for the last 400 ns.

Overall Fold

The typical seven transmembrane helices are maintained across the three systems during the overall trajectory (Panel B in S3 Fig). However, the fold of H8, located at the membrane-cytoplasm interface, which was not solved in the X-ray structure [27], shows membrane-sensitive conformations. Indeed residues 292–317 of H8 preserve a helical conformation in II and III only (Fig 1). Yet, the helical content of H8 decreases in I: residues 305–317 unfolded into flexible loop after 250 ns MD (Fig 1). This is associated with two features, which are present only in I: a very large increase of the Cα RMSD (Panel A in S3 Fig), and the presence of two non-overlapping blocks in the matrix calculated with the all-to-all RMSD analysis (S2 Fig). The latter is a signature of a significant, irreversible transition between two distinct conformations [28]. The different stability of H8 across the three systems is likely to be related to an increase of membrane thickness on passing from I to II, and, more, to III, observed here (see Section ‘Membrane Structure’ below for details). As results, H8 in II and III is only half-exposed to the solvent, being the other half immersed in the membrane. This stabilizes H8, because this helix is amphipathic (Panel B in S4 Fig). Such stabilization is not present for system I, where H8 is more solvent-exposed because of the thinner thickness of the bilayer. H8 is a key structural element for hA2AR function, as it connects the transmembrane helices interacting with ligands with the cytoplasmic C terminus coupling with alpha-actinin (type 2), dopamine receptors (types 2 and 3), glutamate mGlu5 receptors and other regulatory GPCRs [29,30]. Hence, our simulations point to the importance of using proper membrane environment to study this neurotransmitter receptor. Our findings share similarities with an NMR study of the structurally-related class A GPCR human β2 adrenergic receptor, where H8 is helical in DMSO and disordered in water [31]. Also for the class C GPCR mGluR2 receptor [3], mixed POPC/cholesterol membrane was shown to stabilizes the helical structure of the H8, whereas the pure POPC membrane induces a disruption of H8.

Fig 1

Membrane-sensitive folding of H8.

A) The cartoon representations of receptor’s backbone of systems I–III are shown in blue to red according to residues’ increased flexibility, as emerging from the so-called PAD index values [32]; the inserted panel shows the location of H8 (in yellow cartoon) in each membrane; POPC, POPE and cholesterol molecules are shown in red, blue, green lines respectively. The phosphorous atoms are shown as violet Van de Waals spheres. B) Secondary structure content of H8-including C segment (res 292 to 329) (292REFRQTFRKIIRSHVLRQQEPFKAAAAHHHHHHHHHHH329) is reported as a function of the simulated time. β sheet, α helix, coil and bend, and turn are shown in red, blue, white, green and yellow, respectively.

The flexibility of the three systems, here analyzed in terms of the so-called PAD index for backbone atoms [32], is similar for the three systems, with the exceptions of the N-term of helix H1 and the second extracellular loop ECL2 (S5 Fig), which are significantly more flexible in I.

Binding Site

CFF exhibits multiple binding poses (A-D in Fig 2) in receptor binding cavity across the three systems, comparable to what found for adenosine in this receptor [33]. Most of the identified binding poses are similar to those found in the 0.069-μs MD study of a H8-truncated CFF/hA2AR complex [34]. Our poses yet differ from those in the 0.005-μs MD study [35], possibly because of the large difference (more than two orders of magnitude) between our time scale and theirs. The population of the CFF poses depends on the type of membrane environment. A (38%), B (29%), and C (94%) are the most populated poses for I, II, III, respectively (see Fig 2 and Table 2). Notably, the pose in the X-ray structure (D) [27] (Fig 2), is not the most populated one in any of the three systems (see S1 Text). C is almost the only pose assumed by the antagonist in the physiologically relevant system III. It is stabilized by hydrophobic contacts between the CFF/C5 methyl group and ILE66 and SER67 side chains on the extracellular side of H2 (panel A in Fig 3). This stabilizing interaction may be triggered by the diffusion of a cholesterol molecule, already after 0.3 μs, to the cleft between H1 and H2 (panel B in Fig 3). Indeed, this specific cholesterol binding induces conformational rearrangements of VAL57, LEU58, ILE66 and SER67 (S6 Fig), which in turn result in the enhanced stabilization of the hydrophobic interaction between CFF and H2 residues. We conclude that the cholesterol very likely drives specific pose for CFF. The calculated lateral diffusion coefficient of cholesterol molecules around the receptor is not too dissimilar from that of cholesterol molecules in the proximity of the lipids (5*10−8 cm2 s-1, 8*10−8 cm2 s-1 respectively, S2 Text and S7 Fig), suggesting that the observed cholesterol binding event is not strongly dependent on cholesterol’s starting location. Notably, in the absence of cholesterol molecules (systems I and II), one POPC molecule replaces cholesterol in the binding cleft. Hence, this receptor’s binding cleft seems to act like an anti-diffusion trap for lipids and, more, for cholesterol molecules, showing higher specificity for the latter. Interestingly, in the μs-MD simulations of apo hA2AR, Lyman et al. also detected the specific cholesterol presence between helices H1 and H2 [23].

Fig 2

CFF’s most populated binding poses in systems I-III (A-D).

For each binding pose, the upper panel shows the protein backbone in yellow cartoon, CFF and residues interacting with CFF in thick and thin sticks, respectively. Water molecules forming H-bonds with CFF and residues are represented as red sphere; the lower panel shows the corresponding 2-d chart.

Table 2

Populations of CFF binding poses (%) detected across systems I-III over the last 400 ns of MD simulated time.

BP Index	I	II	III
A	38.1		-
B	-	29.1	-
C	31.9	5.8	92.4
D	-	-	3.8

Fig 3

Specific cholesterol binding to hA2AR A) Cartoon showing cholesterol-binding pose in H1/H2 cleft in system III.

The receptor is shown in yellow cartoon; the cholesterol molecule is shown as green sticks surrounded by its solvent accessible surface; CFF, cholesterol-interacting residues, VAL57, LEU58, as well as CFF-interacting residues ILE66, SER67 are shown as green sticks with oxygen and nitrogen atoms colored in red and blue, respectively. B) The diffusion of cholesterol into of the H1/H2 cleft enhances hydrophobic contacts between CFF and H2. The minimum distances between the specific cholesterol molecule and H1 (residues 5–34), between cholesterol and H2 (residues 41–67), between C5@CFF and heavy atoms of ILE66 and SER67 side chains, are shown in black, red, blue and green, respectively.

Next, we investigated the mobility of CFF, i.e. the different extent of roto-translation of the ligand inside the cavity across the three systems over the last 400 ns. This mobility is measured in terms of the orientational flipping angle (defined in the method section, see S8 Fig) and the CFF center of mass (Fig 4), sampled along the simulation time. Not surprisingly, III features the smallest fluctuations of both quantities, as the ligand is mostly in the C conformation. II exhibits the largest fluctuations whilst I features intermediate values (see Fig 4 and S8 Fig).

Fig 4

The distribution of CFF center of mass within the ligand-binding cavity of hA2AR across systems I-III.

The receptor and CFF are shown as yellow cartoon and red sticks, respectively. CFF center of mass at each collected frame over the last 400 ns of MD simulated time is depicted as one black dot.

The MD-averaged number of water molecules in the binding cavity of II (32 molecules) is much larger than those of I and, more, of III (21 and 16 molecules, respectively, see S9 Fig). Hence, this study corroborates the plausible hypothesis that enhanced hydration of binding cavity increases ligand dynamics in II [33].

Membrane Structure

The primary amine group of POPE (in II and III) forms intra and intermolecular hydrogen bonding with the lipids' phosphate groups (S2 Table) [36]. Water- and receptor-lipid hydrogen bonds are comparable across the three systems (S2 Table).

The different membrane composition affects its thickness. The latter increases on passing from I to II, and from II to III. The latter observation is consistent with the experimental observation that the presence of cholesterol causes an increase of thickness of lipid bilayers [37]. These features are shown by a plot of phosphate groups’ density distributions (panel A in S4 Fig). The area per lipid decreases from 0.61 nm2 (system I) to 0.56 nm2 (system II) and 0.50 nm2 (system III), see S2 Table. This is consistent with the fact that the average area per lipid of pure POPC is greater than that of pure POPE membrane [38]. Hence, the area per lipid is anti-correlated with the thickness of the membrane, consistently with what already found in ref. [39].

The POPC and headgroups’ dipole moments turn out to be oriented differently on passing from I to the other two systems. Let us define the PN vector (from the phosphorus atom to the nitrogen atom of one lipid headgroup, see S10 Fig) and the angle (ΦPN) between the PN vector and the axis perpendicular to the lipid bilayer surface (z-axis). The POPC headgroups’ dipoles turn out to be perpendicular to z-axis in system I, with a large standard deviation (ΦPN ~ 90(36) degrees, see S11 Fig). Instead, in systems II and III, these headgroups show a bivariate distribution with two shallow peaks around ΦPN ~ 65 degree and ΦPN ~ 115. Interestingly, the POPE headgroups are oriented perpendicular to the z-axis (ΦPN ~ 90(28) degrees and 90(29) for systems II and III, respectively).

In system III, the POPC/POPE ratio is 1.2:1. This differs from that of system II, which features a POPC/POPE ratio of 1:1. To test whether this difference in POPC/POPE ratio plays a role for protein structural, we performed an additional simulation where we replaced 24 POPE molecules of system II with 24 POPC molecules. The results are rather similar to the ones of system II and are reported in SI (S12 Fig).

Conclusion

Both hA2AR fold and CFF binding dynamics are sensitive to the lipid environment where hA2AR is embedded. The lipid polar headgroups also exhibit varied dipole orientations in different membrane environments. Most importantly, the presence of cholesterol in the membrane is shown to drastically affect CFF binding pose population and mobility. X-ray studies commonly crystallize ligand/receptor complexes in detergent mimics [40], without the physiologically high concentration of cholesterol. The artificial environment is found here to affect the population of ligand poses drastically: the pose found in the X-ray structure, at 3.6-Å of resolution, is the most populated one in none of our 0.8 μs-long MD simulations. This study suggests that computer-aided studies of hA2AR in nearly physiological conditions may give key contributions to the investigation of receptor's function as well as to the development of CFF derivatives that retain CFF's neuroprotective benefits with much higher affinity for the target than CFF [41,42].

Materials and Methods

Homology Modeling of hA2AR

hA2AR is a class A GPCR, composed of 7 transmembrane helices (H1–H7) and a helix lying at the membrane-cytoplasm interface (H8). The X-ray structure of the CFF/hA2AR complex has been solved at a 3.6-Å resolution (PDBid: 3RFM) [27]. Amino acid sequence after residue 317 was deleted to remove the highly mobile cytoplasmic C terminus. The truncated sequence was joined by a polyhistidine tag (residues 318–329). Residues 1MPIMGS6, 150KEGKNHSQ157, 306HVLRQQEPFKAAAAHHHHHHHHHH329 are not detected in the crystal structure. Moreover the crystalized receptor contains 8 mutations (A54L, T88A, R107A, K122A, L202A, L235A, V239A, S277A). The missing regions were complemented and the mutations were mutated back to wild type by multiple-template-based homology modeling (S1 Table).

12 X-ray structures of hA2AR are deposited in the Protein Data Bank [15,27,43-47]. Among these we selected 3RFM and 4 other templates with resolution below 3.0 Å (3VG9, 3EML, 4EIY, 2YDV). Notice that the 2YDV template where hA2AR is bound with an agonist N-ETHYL-5'-CARBOXAMIDO (NEC), is the only one used for modeling the residues 291–325 in the C terminus without including any other residues, as it has the longest resolved helix H8 among all the available hA2AR X-ray structures. 100 models using the 5 templates were generated in Modeller 9.11 [48]. The best model, in terms of both DOPE score [49] and stereochemistry PROCHECK analysis [50] underwent to loop refinement procedure [48]. 500 models were generated. The best model was selected as the optimal initial model for MD simulations.

The first six residues in the N terminus are predicted as a helical segment. The missing 8 residues (150KEGKNHSQ157) in the loop connecting helix 4 and helix 5 are also so as to form a short helical segment. As for the missing residues at the C terminus, residues 306HVLRQQEPFKAAAAHHHH323 are modeled as the H8 [46]. The last six residues, all histidines, are modeled as a loop. The backbone Root Mean Square Deviation (RMSD) between the model and 3RFM is 0.9 Å. The backbone RMSD of the residues located in hA2AR ligand binding site (within 7.0 Å of CFF) against their counterparts in 3RFM is 3.3 Å.

Simulation Details

Membrane models in systems I, II and III were generated using the MemBuilder tool [51]. The inflateGRO code [52] was used to pack lipids around the hA2AR constructed in the previous step. The systems were inserted in a simulation box of size 14.8 nm x 11.1 nm x 10.0 nm. With this choice, the minimum distance between periodic images of the protein was larger than 1.5 nm in all systems. Water, sodium and chloride ions were added in order to solvate and neutralize the systems at an ionic strength of 0.15 M. The final systems comprised ~150,000 atoms (Table 1).

The AMBER99SB-ILDN force fields [53], the Slipids [54,55], the TIP3P [56] force fields were used for the protein and ions, the lipids, and the water molecules respectively. The General Amber force field (GAFF) parameters [57] were used for CFF, along with the RESP atomic charge using Gaussian 09 [58] with the HF-6-31G* basis set [59,60]. MD simulations were performed using Gromacs v4.5.5 package [61] on JUROPA supercomputer. The Particle Mesh Ewald method [62] was used to treat the long-range electrostatic interaction with a real space cutoff of 1.2 nm. A 1.2 nm cutoff was used for the short-range non-bonded interaction. A timestep of 2 fs was set. The LINCS algorithm [63] was applied to constrain all bonds involving hydrogen atoms. Constant temperature and pressure conditions were achieved via independently coupling protein, lipids, solvent and ions to Nosè-Hoover thermostat [64] at 310 K and Andersen-Parrinello-Rahman Barostat [65] at 1 atm. For each system, the receptor in the free state underwent minimization, 1-ns simulated annealing and 10-ns equilibration with positional restraint using a force constant of 1000 kJ mol-1 nm-2 on the heavy atoms of the protein. 40-ns equilibration was further carried out with positional restraint on the side chains of the residues within the binding cavity (residues within 7.0 Å of CFF on backbone alignment to 3RFM [27]). This allowed water molecules to diffuse into the ligand-binding cavity. Next, CFF was inserted so as to fit the conformation it has in the X-ray structure [27] using backbone alignment in Pymol [66]. Energy minimization, annealing, 20-ns equilibration with positional restraint on the side chains of the residues belonging to the binding cavity and the CFF were performed before removing all the restraints. Then 0.8 μs MD at 310 K and 1 atm was performed for I, II, III, with one frame collected every 20 ps. The starting binding pose of CFF resembled fairly that in the X-ray counterpart (RMSD = 1.4 Å, 0.8 Å, 2.3 Å for I, II, III).

Trajectory Analysis

The RMSD, pairwise RMSD matrices [28] and secondary structure content are calculated over the entire trajectory with g_rms, do_dssp of the Gromacs v4.6.5 package [61]. The CFF orientational flipping angle is defined as arccos(μτ·μι), where μι is the vector in the plane of the bicyclic core of CFF, chosen so that it faces toward the extracellular side in the initial frame ι of the trajectory; μτ is the vector at each frame τ of the trajectory. Also this quantity is calculated over the entire trajectory.

The following properties are calculated over the last 400 ns of the three MD simulations: (i) The PAD flexibility index, using a in house code [67,68]. (ii) The density profile of lipid phosphate groups along the z axis, using g_density in Gromacs v4.6.5 package [61]. (iii) The CFF binding poses, identified using the Gromos cluster algorithm [69] with a 2-Å RMSD cutoff on alignment of protein backbone. The g_cluster module of Gromacs v4.6.5 [61] has been used. (iv) CFF center of mass is defined by the vector ri of the coordinates of CFF center of mass at a frame i, upon protein backbone alignment. (v) The average number of water molecules in the ligand binding cavity of hA2AR is calculated by following the similar procedure in [33,70]. Namely, for each system, the number of water molecules in a rectangular box, centered in the binding site and incorporating residues A63, T88, F168, M177, W246, L249, H250, N253, I274, H278 (i.e. residue within 4.5 Å from CFF in 3RFM [27]), is averaged over the trajectory’s frames. (vi) The average hydrogen bond occupancies are calculated as average number of hydrogen bonds formed between the receptor and lipids (), within lipids (), and between lipids and solvent molecules (), divided by the total number of acceptor oxygen atoms. (vii) The average area per lipid (APL) is calculated with the GridMat-MD program [71]. (viii) The PN vector of one lipid molecule is defined as the vector from the phosphorous atom to the nitrogen atom of its polar headgroup, as defined in [72]. The ΦPN of one lipid molecule is the angle between its PN vector and the z-axis. This is the axis perpendicular to the lipid bilayer surface. For each system, ΦPN of each lipid molecule is sampled over the entire MD simulation for the normalized ΦPN distributions of POPC and POPE lipids. (ix) the lateral diffusion coefficients of cholesterol molecules in III are calculated using the Einstein relation [73] for a 60 ns simulation of III in the NVT ensemble (see S2 Text for details). Cholesterol molecules are classified as molecules in close proximity of the protein if they have atoms within 0.35 nm of the receptor during the dynamics. The other molecules are classified as 'free'.

Molecular graphics are drawn using Pymol [66], VMD [74] and Ligplot+ [75].

Chemical structures of CFF, POPC, POPE and cholesterol molecules.

(PDF)

Click here for additional data file.

Pairwise RMSD matrix of Cα atoms.

(PDF)

Click here for additional data file.

Proteins' backbone.

(PDF)

Click here for additional data file.

Selected properties of systems I-III.

(PDF)

Click here for additional data file.

Flexibility of individual residues of hA2AR.

(PDF)

Click here for additional data file.

Cholesterol-induced conformational transitions of H2 residues in system III.

(PDF)

Click here for additional data file.

Lateral Mean-square displacements (MSDs) of two groups of cholesterol molecules in system III.

(PDF)

Click here for additional data file.

CFF orientational flipping angle.

(PDF)

Click here for additional data file.

The hydration of the ligand binding cavity of hA2AR.

(PDF)

Click here for additional data file.

Orientation of POPC and POPE headgroups.

(PDF)

Click here for additional data file.

Distribution of ΦPN of lipid headgroups for systems I-III.

(PDF)

Click here for additional data file.

Changing the ratio between the lipids in system II.

(PDF)

Click here for additional data file.

Available X-ray structures of hA2AR.

(PDF)

Click here for additional data file.

Selected membrane properties of systems I-III.

(PDF)

Click here for additional data file.

Description of most populated CFF binding poses A-D.

(PDF)

Click here for additional data file.

Lateral diffusion coefficient of cholesterol molecules.

(PDF)

Click here for additional data file.