A functional selectivity mechanism at the serotonin-2A GPCR involves ligand-dependent conformations of intracellular loop 2.

With recent progress in determination of G protein-coupled receptor (GPCR) structure with crystallography, a variety of other experimental approaches (e.g., NMR spectroscopy, fluorescent-based assays, mass spectrometry techniques) are also being used to characterize state-specific and ligand-specific conformational states. MD simulations offer a powerful complementary approach to elucidate the dynamic features associated with ligand-specific GPCR conformations. To shed light on the conformational elements and dynamics of the important aspect of GPCR functional selectivity, we carried out unbiased microsecond-length MD simulations of the human serotonin 2A receptor (5-HT(2A)R) in the absence of ligand and bound to four distinct serotonergic agonists. The 5-HT(2A)R is a suitable system to study the structural features involved in the ligand-dependent conformational heterogeneity of GPCRs because it is well-characterized experimentally and exhibits a strong agonist-specific phenotype in that some 5-HT(2A)R agonists induce LSD-like hallucinations, while others lack this psychoactive property entirely. Here we report evidence for structural and dynamic differences in 5-HT(2A)R interacting with such pharmacologically distinct ligands, hallucinogens, and nonhallucinogens obtained from all-atom MD simulations. Differential ligand binding contacts were identified for structurally similar hallucinogens and nonhallucinogens and found to correspond to different conformations in the intracellular loop 2 (ICL2). From the different ICL2 conformations, functional selective phenotypes are suggested through effects on dimerization and/or distinct direct interaction with effector proteins. The findings are presented in the context of currently proposed hallucinogenesis mechanisms, and ICL2 is proposed as a fine-tuning selective switch that can differentiates modes of 5-HT(2A)R activation.

Introduction

The G protein-coupled receptors (GPCRs) have key physiological roles in a large variety of signaling and cell–cell communication processes and hence are primary targets for a large percentage of therapeutic drugs.[1−3] Intense research focus on GPCRs, and especially on the rhodopsin-like class A, has produced significant insights into their structure and function, and various interpretations and models for ligand binding and receptor activation have been proposed.[3] Emerging GPCR activation models suggest a highly versatile receptor capable of signal transduction, with different efficacies, through various downstream signaling pathways in a ligand-specific manner.[4−6] The ligand-dependent differential efficacy for distinct downstream cellular responses initiated via the same receptor is well documented for many GPCRs[4,7] and has been named “functional selectivity” or “biased agonism”,[4] but for the most part the structural basis for this selectivity is not well understood. One way to address this gap in understanding is to identify the GPCR structural elements that adopt distinct ensemble conformations in response to ligands with different pharmacological properties. Notably, this would also have profound practical implications in the development of more effective therapeutic drugs that target specific signaling pathways with a minimum of detrimental side effects. We present here the results from such an investigation, targeting GPCRs in the family of serotonin (5-hydroxytryptamine; 5-HT) receptors.

Of the 15 different receptors activated by the neurotransmitter serotonin, the 5-HT2A subtype is of great interest not only because it plays a crucial role in cognitive processing but also because it is the target of a large number of medications including antidepressants and antipsychotics.[8−10] Remarkably, several 5-HT2A agonists, such as the classical psychotomimetic LSD,[11] are known to display hallucinogenic properties. Indeed, a large body of evidence indicates that the common target of all hallucinogens is the 5-HT2A receptor (5-HT2AR).[8,10,11] Functional selectivity by hallucinogenic (HL) and nonhallucinogenic (NHL) agonists has been measured for 5-HT2AR-mediated activation of Gαq/11 and pertussis toxin-sensitive Gi/o protein,[8] which signal through different pathways.[12,13] Furthermore, 5-HT2AR-mediated functional selectivity has been observed in β-arrestin-dependent signaling pathways.[14] Yet despite extensive research and significant progress in recent years, a mechanistic understanding of the distinct responses elicited by 5-HT2AR agonists that are known to produce LSD-like hallucinations (e.g., LSD, DOI, psilocybin, mescaline), and the differences in the manner in which they elicit the underlying activation of that receptor, remains elusive.

The remarkable functional selectivity of HL compounds on 5-HT2AR,[8,12−14] included in the ample literature on the experimentally determined properties of the receptor and of structure–activity relations for its ligands,[11] prompted us to investigate structural and dynamical elements associated with the functional selectivity of the HL and cognate NHL 5-HT2AR agonists[8] using computational modeling and simulation. To cover a chemically distinct ligand space, we selected 5-HT2AR complexes with the four agonists (i–iv) described below for extensive unbiased all-atom molecular dynamics (MD) simulations. Thus, we studied the receptor in complex with two HL compounds: (i) the hallucinogenic substituted amphetamine, 2,5-dimethoxy-4-iodoamphetamine (DOI), and (ii) the prototypical hallucinogen l-lysergic acid diethylamide (LSD). DOI has a relatively small and flexible chemical structure, whereas LSD is representative of the larger and more rigid chemical family of ergots (see Figure S1A in SI). We also selected two cognate NHL compounds: (iii) the endogenous 5-HT2AR ligand serotonin (5-HT), and (iv) the partial agonist R-lisuride (LIS) that belongs to the same chemical family as LSD but has a very different pharmacological fingerprint with regards to perceptional and cognitional phenotypes.[8] All four compounds have been extensively characterized with diverse biophysical and physiological techniques in vitro and in vivo with respect to serotonergic signaling efficacy across several downstream pathways and hallucinogenic phenotypes.[8]

From the studies described herein we identified a special role for the second intracellular loop (ICL2) in observed hallucinogen-specific conformations. ICL2 is a structural element known to be involved in activation of downstream signaling, and detailed analysis of the results leads to a proposed role for ICL2 in mediating ligand-specific receptor responses to HL versus NHL ligands. Thus, the functional selectivity evidenced by the different pharmacological outcomes of the binding of HL as compared to NHL ligands is proposed to emerge from molecular mechanisms affected by the divergent properties of the ICL2 in HL-bound versus NHL-bound 5-HT2AR. The affected molecular mechanisms relate to the propensity of the ligand-bound 5-HT2AR systems for selective dimerization and/or direct interaction with different effector proteins as a result of differences in ICL2 properties.

Results

An all-atom model of the 5-HT2AR generated with homology modeling and embedded in a physiologically relevant cholesterol-containing membrane model was used as the starting point of five different microsecond-long MD simulations for the following ligand-bound complexes of the 5-HT2AR: 5-HT2AR/5-HT, 5-HT2AR/LIS, 5-HT2AR/DOI, 5-HT2AR/LSD, and the ligand-free state 5-HT2AR (APO), (see Methods section and Figure S1 in SI). Both the 5-HT2AR molecular model and the ligand contacts were found to be in very good agreement with the published X-ray structures of the related 5-HT1B and 5-HT2B serotonin receptors (see text in SI). Conformational elements known to be associated with GPCR activation were analyzed[15−17] from the second half of each of the equilibrated trajectories, i.e., from 500 to 1000 ns, unless otherwise indicated (Figure S1B in SI).

The ionic lock between R3.50 and E6.30 (superscripts identify throughout the residues in the generic Ballesteros and Weinstein numbering system,[18] see Methods) is a known landmark in GPCR activation[19] and was found here to undergo open-to-close dynamics in all the agonist-bound systems but not in the APO construct where it remains closed for the duration of the simulation (Figure S2 in SI). Also, the χ1 angle of the aromatic residue W6.48 that triggers the toggle switch proposed as another key element of the activation[20,21] changed similarly in all the agonist-bound systems but not in the APO. These conformational changes elicited by GPCR agonists, but not by antagonists or inverse agonists, have been described and associated with active-like protein conformations.[3,22] Not surprisingly, all the ligand-bound systems behave similarly regarding the aforementioned activation conformational elements since all of them, HL and NHL, are 5-HT2AR agonists.

To identify structural and dynamic changes that would not be shared among the ligands with identified functional selectivity properties, such as the HL versus NHL, we broadened the analysis from changes generally associated with GPCR activation (including the ionic lock and arginine–cage interactions and also toggle switch at W6.48, NPxxY motif, and Y5.58–Y7.52 interactions[15,17,19]) to address as well rearrangements in portions of the GPCR that are closer to the intracellular region and are likely to modulate the interface with effectors and regulators (e.g., G proteins, arrestins, etc.). As described herein, we identified particular properties of interest in the dynamics of the ICL2 which, in the class A subfamily, has been found to play a major role in the interaction of the receptor with different intracellular effector proteins.[3,23−25]

ICL2 Adopts Distinct Conformations in 5-HT2AR Complexes with Different Ligands

From the analyses of the microsecond MD simulation trajectories of 5-HT2AR with different ligands, we found that ICL2 conformations favored in the HL-bound systems are different from those favored in the NHL-bound and in the unbound constructs. The distinct conformations were monitored in an internal coordinate system defined by taking advantage of the observation that the secondary structure content of ICL2 remains helical for most of the trajectories obtained for the five systems (see Figure S6 in SI). Using the internal coordinate system analysis (see Figure S7A and text in SI), we uncovered the differences in the conformations of ICL2 sampled by the receptor binding HL compounds (LSD and DOI) compared to those visited when NHL are bound (5-HT and LIS) or in the APO form (compare blue and yellow bars in Figure S7B in SI). The ICL2 conformations were also characterized by defining the center of mass of the helical segment on ICL2 as a collective variable and calculating the root-mean-square deviation (rmsd) of the center of mass of the ICL2 along the trajectories, relative to the center of mass of the ICL2 in the initial structure. The distributions of the rmsd values show two distinct conformations for the ICL2 (Figure 1A), consistent with the analysis in the internal coordinate framework (see Figure S7 in SI). The more “outward” and more “upward” oriented ICL2 conformations (colored blue in Figure 1B) are seen to be highly favored by HL (DOI and LSD, see middle panel in Figure 1B) in contrast to the more “inward” and more “downward” ICL2 conformations (colored yellow in Figure 1B) adopted when the NHL (5-HT and LIS) are bound or when the unbound (APO) receptor is simulated. The representative structures of the ICL2 segment conformations in each of the studied systems (Figure 1B) show that the more outward conformations (favored by HL) situate the ICL2 segment farther away from the axis of the TM helical bundle, whereas more upward conformations place the ICL2 segment closer to the center of the membrane bilayer (see also Figure 2). Representative structures in the 5-HT2AR/DOI complex are also depicted in Figure 2. In this particular complex, ICL2 selectively prefers more outward-upward conformations (colored blue) but explores as well the inward-downward ICL2 conformations preferred by the NHL (colored yellow), see Figure 2. All ligand-bound receptors exhibited dynamic transitions between states, but with notable preferences related to their pharmacological class (Figure 1B).

Figure 1

RMSD distribution and representative structures of ICL2. (A) The distributions of rmsd values, relative to the starting structure, are shown for the five simulated systems, 5-HT2AR/5-HT, 5-HT2AR/LIS, 5-HT2AR/LSD, 5-HT2AR/DOI, and 5-HT2AR (APO), respectively. The more outward-upward conformations (blue) are highly favored in just the hallucinogenic systems. (B) Representative ICL2 structures for the five simulated systems, as seen from the intracellular side, are shown. As a reference, the initial structure (gray) is also depicted in each case. The more outward-upward ICL2 conformations are colored blue, whereas the more inward-downward conformations are colored yellow (see also Figure 2 for a different perspective of the ICL2 conformations). In these views, the more outward ICL2 conformations correspond to larger values in the X-axis coordinate. Interestingly, the outward-upward conformations (blue) are preferentially stabilized in the hallucinogenic systems, LSD and DOI. The thickness of the ICL2 representation corresponds to the percentages of the distributions from (A). In the case of LIS, any of the conformations was sorted as part of the “blue” conformations. The helical axis of the TM bundle is represented by a magenta triangle in each case.

Figure 2

Conformations explored by the ICL2 in the 5-HT2AR/DOI complex. Lateral view of two representative structures of the 5-HT2AR/DOI system. The ICL2 structures are the same as those displayed in Figure 1B for this system (DOI). The more outward (relative to the helical axis of the TM bundle, shown here as a magenta line) and more upward (that is, closer to the center of the lipid bilayer) are preferred in the hallucinogenic systems (DOI and LSD), colored here in blue. In these views, the more outward ICL2 conformations correspond to larger values in the X-axis coordinate, while the more upward conformations correspond to larger values in the Z-axis coordinate. The more inward-downward conformations are preferentially sampled in the nonhallucinogenic systems (5-HT and LIS) and in the APO form, colored here in yellow. As a magnitude reference, the Cα atoms of residue H183 are depicted in both structures, and the distance for these particular structures is 5.3 Å (indicated as a red line). The intracellular boundary of the bilayer, as predicted by the orientations of proteins in membranes (OPM) database,[26] is depicted as a brown line.

The spatial similarity of the ICL2 conformations adopted preferentially in each of the different systems was investigated by generating the “clouds” of points representing locations sampled by the center of mass of the helical segment in ICL2 along the trajectories (see Figure 3A–H and Figure S8 in SI). To quantify such spatial similarities between the different conformational distributions, we calculated the overlap coefficient[27] of the center of mass of the helical part of the ICL2 (each of the points in the different clouds) by using eq 1 (see Methods). The overlap coefficient takes values from 0 (no overlap) to 1 (perfect overlap). As indicated by the overlap coefficient values (Figure 3I and Table S2 in SI), the conformations of ICL2 closely overlap for the HL systems (see LSD and DOI). On the other hand, the overlap between any of the HL systems with any of the NHL systems and the APO form is minimal (indicated by the low overlap coefficient values). Similarly, the NHL systems show significant overlap but have minimal overlap with any of the HL systems. This analysis confirms the similarity of ICL2 conformations for the receptor complexed with ligands with the same pharmacological phenotypes, i.e., HL or NHL. The distinct spatial overlap of the ICL2 conformations mediated by HL or NHL underscores the significance of the finding of distinct ICL2 conformations for the HL systems. The close spatial overlap between the ICL2 conformation from the 5-HT2AR/DOI and the 5-HT2AR/LSD complexes, quantified by the overlap coefficient, indicates that the ICL2 conformations of the systems containing the two HL are spatially more similar to each other than to any of the NHL complexes or the APO system.

Figure 3

Conformations explored by the ICL2. As a reference, the initial structure (gray) is also depicted in each case. (A–D) Intracellular views of the “clouds” of points that represent the conformation adopted by the ICL2 along the trajectories for the systems, 5-HT2AR/5-HT, 5-HT2AR/LIS, 5-HT2AR/DOI, and 5-HT2AR/LSD, respectively. The same region of space, in the XY-plane, is indicated for all the cases by blue ellipses. This region is explored preferentially in the hallucinogenic systems (C) 5-HT2AR/DOI and (D) 5-HT2AR/LSD. Here these regions represent ICL2 conformations that are positioned more outward, relative to the helical axis of the TM bundle (indicated by magenta triangles). This was not the case for the nonhallucinogenic counterparts or the APO system that explored more inward conformations. (E–H) Lateral views of the same systems. The same region of space, in the YZ-plane, is indicated by blue circles. The ICL2 conformations favored by the hallucinogens (DOI and LSD) are located more upward relative to the center of the membrane bilayer (larger positive values for the Z-coordinate in these orientation) when compared with the conformations favored by the nonhallucinogenic (5-HT and LIS) or the APO simulations. Here the intracellular boundary of the bilayer, as predicted by the OPM database,[26] is depicted as a brown line. The relative coordinate system is indicated in all the cases. Note that the cycles in the intracellular views and in the lateral views do not represent the same spatial region. Equivalent views for the APO system are shown in Figure S8 in SI. (I) Using the position of each of the points in the different clouds, we calculated the overlap coefficient for each of the five systems using eq 1.

To investigate differences in ICL2 dynamics in all five systems we performed principal component analysis (PCA) (see Methods). We first calculated the generalized correlation coefficient[28] between the center of mass and the first principal component (PC1) of the ICL2 motion, which indicated that the center of mass motion of ICL2 was strongly correlated with the PC1 in each system (see Table S1 in SI). This finding further supports the use of the center of mass as a collective variable, as described above. From the PCA we further found that the PC1 motion in the 5-HT2AR/DOI system accounts for a large fraction of the variation present in 5-HT2AR/DOI and 5-HT2AR/LSD systems, but not in the 5-HT2AR/5-HT, 5-HT2AR/LIS, or 5-HT2AR (APO) systems, indicating that this motion is HL-specific (see Table S1 in SI). These results, in combination with the previous analysis of the center of mass, bring into evidence the difference in the conformations and dynamics of ICL2 in our simulations of 5-HT2AR bound to the two groups of ligands exhibiting functional selectivity and identify ICL2 properties that are specific to the type of system involved, i.e., HL- or NHL-bound or APO.

To identify specific molecular interactions involved in the observed differential conformations of the ICL2 segment, we analyzed comparatively the contacts involving residues in ICL2. The direct interaction between residue D1723.49 (from the conserved DRY motif) and H183 (located in the middle of the ICL2) was found to be more extensively maintained in the trajectories of HL systems compared to the NHL counterparts or the APO (see Figure 4A). Figure 4B shows that the minimal distance between any of the carboxylate oxygen atoms from the side chain of D1723.49 with any of the imidazole nitrogen atoms from the side chain of H183ICL2 in the HL systems fluctuates mainly to values ∼4 Å or shorter. In contrast in the NHL systems the values are mostly larger than 4 Å (the 4 Å is selected as reference distance to match the cutoff distance value used herein to define a molecular contact). This interaction is proposed to play a key role in determining the different conformational and dynamic properties of the ICL2 in the HL versus NHL systems.

Figure 4

Distances of residues D172 and H183. (A) Representative structure from the 5-HT2AR/DOI complex where the interaction of D1723.49 and ICL2 residue H183ICL2 is depicted. In the context of the DRY motif, this position in the ICL2 is a residue that can establish polar interactions by using its side chain and is located in the sequence position Z in the “DRY(X)5P(X)2Z” motif. (B) The minimal distances between any of the carboxylate oxygen atoms from the side chain of D1723.49 with any of the imidazole nitrogen atoms from the side chain of H183 are depicted. The distance (gray) and its moving average (black) are displayed. As a reference, a dashed line at 4 Å is also displayed (the same cutoff value used to define a receptor–ligand interaction contact).

Binding Site Interactions of 5-HT2AR Agonists

The receptor–ligand contacts were evaluated by considering all positions at which any heavy atom from the ligand comes within 4 Å of any heavy atom from the protein in the course of the trajectory (Figure S9 in SI). From this set we identified six contact loci (I2104.60, G2385.42, S2395.43, F2435.47, W3366.48, and N3436.55) that were found in all the simulated complexes but that exhibited differences in the frequency of contacts for HL versus NHL ligands (Figure 5A,B). An additional position (W1513.28) was also found to have differential contact frequencies between HL and NHL but only in the case of the larger ergoline ligands, LSD and LIS (Figure 5A,B). Figure 5A,B depicts the seven residues in the context of their positions inside the binding site, whereas Figure 5C displays their respective contact frequencies as the percentage of trajectory time in which each of the positions is in contact with the ligand. The location of this set of residues suggests that HL agonists preferentially interact with residues located in TM6, whereas their NHL counterparts preferentially establish contacts with residues in TM4 and TM3 (Figure 5A). Both classes of compounds interact with residues in TM5, but the HL preferentially contact residues that are located at the helical interface formed with TM6, whereas the NHL contact residues located at the helical interface formed with TM3 and TM4 (Figure 5A). Residues G2385.42 and S2395.43 present an interesting example of this selectivity because they occupy neighboring positions in the vicinity of the indole nitrogen of the 5-HT ligand (or equivalent atoms in the other ligands), see Figure 5A. Yet, position 5.42 is preferentially contacted by NHL (94%, 98% for 5-HT, and LIS versus 37% and 79% for DOI and LSD, respectively), whereas position 5.43 is contacted more extensively by HL compounds (24%, 31% for 5-HT, and LIS, versus 82% and 54%, by DOI and LSD, respectively). Interestingly, even though all the ligands contact the same residues in the orthosteric binding site (albeit with different frequencies; Figures S4 and S9 in SI), two of the residues preferentially contacted by the HL are large aromatic amino acids that are located deep in the orthosteric binding pocket, i.e., the highly conserved W3366.48, known to be implicated in signal transduction in different GPCRs,[22] and F2435.47, known to modulate DOI-dependent downstream signaling in 5-HT2AR (Figure 5B).[29]

Figure 5

Ligand binding contacts in the 5-HT2AR. (A) Extracellular and (B) lateral views that show the seven residues (W1513.28, I2104.60, G2385.42, S2395.43, F2435.47, W3366.48, and N3436.55) that display preferential frequency contacts between HL (blue) and NHL (yellow) ligands. (C) The percentage of time that each of the seven positions are in contact with the ligands along the trajectories are shown. Similar color code is used, 5-HT (orange), DOI (red), LIS (green) and LSD (cyan). The different agonist types are arranged: “small” agonists (first row), “large” agonists (second row), NHL (first column), and HL (second column). To discern contact frequency differences between NHL and HL compounds compare data in the different columns in each case. Similarly, by discern contact frequency differences between “small” and “large” agonists compare data in the different rows. The first three residues show a tendency to directly interact with NHL (W151 only interacts with the ergoline ligands, LSD and LIS), whereas the other four show a preference for the HL.

It is noteworthy that in spite of the minimal chemical and structural similarity of the HL ligands (Figure S1B in SI), they both have a positively charged nitrogen atom and an indole-like nitrogen atom (or equivalent) which have long been considered to be particularly important in interacting in the 5-HT2AR orthosteric binding site.[29,30] This is also the case for the NHL ligands (Figure S1B in SI). The lack of chemical and structural similarity within the groups, and the much greater similarity of compounds belonging to the different groups (cf. LSD and LIS), accentuates the significance of the identified common set of residues that establish different protein–ligand contacts in the HL versus the NHL systems.

Potential Implications of Distinct ICL2 Conformations

Dimerization

We considered the possible implications of the differences observed in the dynamic behaviors of ICL2 in the context of mechanisms that could be involved in differential signaling. One such mechanism involves the spatial organization of the 5-HT2AR-s in the cell membrane and in particular the homo- and heterodimerization potential. The physiological role of GPCR dimerization in signaling has been previously demonstrated.[31−33] In the case of 5-HT2AR, a heterodimer unit with the metabotropic glutamate 2 receptor (mGlu2R) has been proposed to be implicated in hallucinogenic responses.[34,35] The roles of specific dimer interfaces have also been suggested by the X-ray solution of different GPCR crystal structures.[33,36−40]

We therefore evaluated the possible role of the ICL2 conformational preferences in the stabilization of 5-HT2AR dimerization. To this end, we used the β1 adrenergic receptor homodimer as a structural template to construct pseudodimer structures of 5-HT2AR conformations resulting from our various MD simulations by aligning the simulated GPCR structures at each point in time, with each of the β1 protomers (PDB accession code 4GPO, TM4-TM5-ICL2 interface). For each of the five systems, we then constructed the histograms of the distance between the centers of mass of the ICL2 in the protomers of the pseudodimer (see Figure 6). Taking as reference the distance (11.2 Å) between the centers of mass of the ICL2 segments between the protomers in the dimer crystal structure of the β1 adrenergic receptor (see also Figure S10 in SI), we observed that the more outward-upward ICL2 conformations preferentially favored in the HL systems (DOI and LSD) would bring ICL2 residues in one protomer closer to their counterparts in the other protomer (gray bars in Figure 6). On the other hand, in the more inward-downward ICL2 conformations preferred by the NHL and APO systems, the separation between the ICL2 segments is larger compared to the crystallographic reference distance (white bars in Figure 6). Consequently, the outward-upward ICL2 conformations that are favored in the HL systems would favor dimer formation by exposing the ICL2 protein surface more when compared with the ICL2 exposure in the NHL-bound or APO receptor. This would increase the probability for the ICL2 to participate in this particular 5-HT2AR dimerization interface (TM4-TM5-ICL2) for the HL-bound states.

Figure 6

Dimer formation modulated by ICL2. After alignment to the crystal structure of the β1 adrenergic receptor in its ligand-free form (4GPO.pdb), the distance of the center of mass of ICL2 in each of the simulated systems was calculated. In the crystal structure of the β1 adrenergic receptor, the distance of each of the protomers center of mass of the ICL2 was 11.2 Å (indicated in red in the APO distribution). The distributions are shown for each case, where the gray area corresponds to distances that could favor dimer formation by this interface (base just on distance considerations and using a 12 Å cutoff). Interestingly, dimer formation by this particular interface will be favored in the case of the hallucinogenic systems, DOI and LSD. The values displayed correspond to the relative frequency in each case.

Interaction with Effectors

We reasoned that the distinct ICL2 conformations observed in the simulations could potentially result in differential accessibility to the intracellular side of 5-HT2AR. To address a possible preference for different intracellular signaling 5-HT2AR partners from this perspective, we calculated structural parameters that could relate to effector binding along the MD trajectories for each system. First, we investigated the structural properties of the intracellular-facing end of the receptor surface and used CASTp[41] to identify the residues forming a cavity observed in the crystal structure of the β2 adrenergic receptor in complex with the Gs protein (see Figure S11A in SI).[42] Then, by using the equivalent residues in the 5-HT2AR we characterized a similar cavity in the region (Figure 7A) and calculated the solvent accessible surface area (SASA) values of the component residues along the different trajectories. The average values for the distributions are 5-HT (1352.8 Å2), LIS (1484.6 Å2), DOI (1219.2 Å2), and LSD (1292.5 Å2). The NHL-bound systems display larger average values when compared with the HL-bound counterparts. To quantify the difference in the distribution of surface accessibility values between the HL- and NHL-bound receptor systems, we integrated the distributions into different aggregates and calculated the overlap between these combinations (see Methods section). The resulting overlap coefficients for the three pairs are 5-HT + LIS versus DOI + LSD (0.546), 5-HT + LSD versus DOI + LIS (0.932), and 5-HT + DOI versus LSD + LIS (0.842). Here, a clear distinction between the different distributions is only observed in the first case, that is, when the agonists displaying pharmacologically distinct phenotypes at 5-HT2AR (HL and NHL) are grouped together (Figure 7B). The distribution of the accessibility to the intracellular cavity for the HL-bound systems is notably shifted toward lower values than for the NHL counterparts (Figure 7B), suggesting that the NHL is more likely to support receptor interactions similar to those observed in the β2 adrenergic receptor/Gs protein complex. When the single residue accessibility of the serine residue S188, which has been implicated in the process of agonist-mediated desensitization of 5-HT2AR,[43] was subjected to similar analysis in the individual complexes, the average values for the distributions were found to be 5-HT (50.9 Å2), LIS (45.4 Å2), DOI (55.4 Å2), and LSD (60.6 Å2). The overlap coefficients for the three unique pairs are 5-HT + LIS versus DOI + LSD (0.734), 5-HT + LSD versus DOI + LIS (0.916), and 5-HT + DOI versus LSD + LIS (0.984). A clear distinction between the different aggregate distributions is again observed for the “HL versus NHL” case (Figure 7C), and the accessibility distributions for S188 is shifted toward higher values in the HL-bound systems relative to the NHL (Figure 7C).

Figure 7

Intracellular cavity in the 5-HT2AR. (A) The cavity in the 5-HT2AR, based on the cavity identified in the β2 adrenergic receptor/Gs-protein complex, is depicted as a gray color surface (see also Figure S11A in SI). The accessibility to this intracellular cavity along the MD trajectories was determined by calculating their SASA values. The position of residue S188 located at the end of ICL2 is also depicted. (B) The aggregate distributions of the cavity accessibilities calculated for the 5-HT2AR complexed with each class of agonists (HL: LSD + DOI) and (NHL: 5-HT + LIS) are shown in blue and yellow, respectively. The distribution of values for the HL systems peaks at lower accessibility values compared to the NHL distribution, which is shifted to the right. (C) The accessibility of the functionally relevant residue S188 in ICL2, as calculated by its SASA, also shows a difference in the aggregate distributions between HL and NHL. In this case, the distribution of values for the HL systems peaks at higher accessibility values relative to the NHL distribution.

Discussion

With recent progress in determination of GPCR structure with crystallography, a variety of other experimental approaches (e.g., NMR spectroscopy, fluorescent based assays, mass spectrometry techniques) are also being used to characterize state-specific and ligand-specific conformational states.[44−47] MD simulations offer a powerful complementary approach to elucidate the dynamic features associated with ligand-specific GPCR conformations. To shed light on the conformational elements and dynamics of the important aspect of GPCR functional selectivity,[4] we carried out unbiased microsecond-length MD simulations of the human serotonin 2A receptor (5-HT2AR) in the absence of ligand and bound to four distinct serotonergic agonists. The 5-HT2AR is a suitable system to study the structural features involved in the ligand-dependent conformational heterogeneity of GPCRs because it is well-characterized experimentally and exhibits a strong agonist-specific phenotype in that some 5-HT2AR agonists induce LSD-like hallucinations, while others lack this psychoactive property entirely.[8,10,11]

That the signaling response of the 5-HT2AR, identified in vitro and in vivo, depends on the inherent nature of the agonist[8] makes this system suitable for the study of functional selectivity. Coupling to the heterotrimeric Gαq/11 protein results in activation of phospholipase C (PLC) causing an increase in the accumulation of inositol phosphates and calcium mobilization.[8] When bound to HL agonists such as LSD, 5-HT2AR can also signal via pertussis toxin-sensitive Gi/o protein, resulting in activation of phospholipase A2 (PLA2) and accumulation of arachidonic acid.[8] The variety of canonical signaling pathways mediated by the agonist-dependent activation of 5-HT2AR is also evident downstream, since differences in transcription factor expression in cells are produced by its interaction with either HL or with NHL agonists.[8] In addition to signaling via G proteins, 5-HT2AR signals via nonclassical pathways mediated by β-arrestin.[14] This interaction is also mediated differentially by the nature of the agonists that interact at 5-HT2AR.[14]

Our findings described here, identifying a role for the second intracellular loop of the 5-HT2AR in discriminant pathway activations, are consistent with previous observations about the signaling of class A GPCRs through various intracellular signaling partners.[3,23−25] Thus, the ICL2 of the 5-HT2AR has been shown to be involved in the interaction with G protein (including desensitization)[43] and with β-arrestin,[25] whereas for the related serotonin 1A receptor, ICL2 has been directly implicated in G protein coupling.[48] The more recent structural information, for the β2 adrenergic receptor complexed with the Gs protein, shows the ICL2 establishing extensive interactions with the β2/β3 loop in the N-terminus of the Gα subunit and with the C-terminus of helix α5.[6,42] In this context, the extensive unbiased MD simulations presented here provide evidence that different ligand classes bound to the 5-HT2AR can produce distinct conformations of the ICL2. Thus, ICL2 favors more outward-upward conformations in the HL-bound systems (i.e., the 5-HT2AR/DOI and 5-HT2AR/LSD complexes), while these conformations are not highly explored in the NHL systems or in the unbound receptor. The spatial distributions of the ICL2 conformations relative to the helical bundle are similar among the HL systems (DOI and LSD), as quantitatively depicted by the calculation of the overlap coefficient of the ICL2 center of mass and the projections of the principal components, and are different from those adopted by the NHL counterparts (5-HT and LIS) or the unbound receptor. This is consonant with previous results from Lefkowitz and co-workers who used quantitative mass spectrometry to identify ligand-specific conformations of the β2 adrenergic receptor and found that ICL2 adopts distinct conformations that differ between agonists.[44]

Our computational analysis shows that the ICL2 conformations are likely to be largely dependent on the extent of the interaction between D1723.49, from the conserved DRY motif in TM3, and H183 in the ICL2. Interestingly, interactions of D1723.49 with H183-equivalent residues in ICL2 has been observed in the crystal structures of several other GPCRs: in all the opioids and the aminergic muscarinic receptors (with an Arg in the corresponding position),[3] and in the serotonin 1B receptor (with a Tyr in that position).[49] Moreover, in another related GPCR, the aminergic β1 adrenergic receptor, a hydrogen bond is formed between D1723.49 and a tyrosine residue (Y149) in the equivalent ICL2 position, and introduction of the Y149A mutation, decreases receptor stability.[50] The relevance of this interaction is further emphasized by the fact the in the β2 adrenergic receptor, the phosphorylation of the equivalent tyrosine (Y141) shifts the conformational equilibrium so as to facilitate active state conformations.[51] In the context of the DRY motif in TM3, this particular ICL2 position is located in the sequence position Z in the “DRY(X)5P(X)2Z” motif, where in all the aforementioned examples position Z is a residue able to establish side chain polar contacts with D1723.49. In the context of the present analysis, it is tempting to suggest the probing of the D1723.49–H183ICL2 interaction by mutagenesis as a validation test of the predicted role of ICL2 conformations in the functional selectivity of the 5-HT2AR ligands we studied.

The distinct ICL2 conformations favored in the HL-bound systems are correlated with distinct protein–ligand interactions identified in the ligand binding site. The ligands considered in the simulations include two structurally and chemically related ergoline ligands, LSD (HL) and LIS (NHL), and two ligands with relatively smaller and more flexible chemical structures, the synthetic substituted amphetamine DOI (HL) and the simple tryptamine 5-HT (NHL). In spite of this chemical diversity within each set of pharmacologically diverse ligands (i.e., LSD and DOI, compared to LIS and 5-HT; see Figure S1B in SI) each set established common protein–ligand contacts which, however, differed for the HL versus the NHL compounds. This independence of structural and chemical features from the outcome of the protein–ligand interactions strengthens the validity of the effects on the GPCR structural dynamics as the significant determinant of the pharmacological phenotype. Moreover, the observation that the interaction differences relate to a small set of six residues in the binding site makes it possible to validate the computational results experimentally by addressing the residues that are involved in differential contacts for the two sets.

The functional implications of the salient difference we observed in the structural dynamics of ICL2 of 5-HT2AR bound to members of the two different sets can be understood in the context of available information regarding structure and function of the biogenic amine receptors. Thus, we proposed that the outward-upward conformations of ICL2 may mediate homo- or heterodimerization of 5-HT2AR[34,52] as the ICL2 conformations favored in the HL-bound systems would be more likely to mediate protein–protein interactions involving residues in the ICL2 and support dimerization at the TM4-TM5-ICL2 interface (Figure S10 in SI).[33,53] In this context, the involvement of ICL2 (but not ICL3) in functionally active homodimer constructs has been demonstrated for the dopamine D2 receptor.[31] Also, the role of ICL2 has an important implication in view of the proposed role of heterodimer formation between the 5-HT2AR and the metabotropic glutamate 2 receptor (mGlu2R) as the physiologically functional unit in hallucinogenesis[34,35] (indeed Sealfon, Gonzalez-Maeso, and co-workers identified the TM4-TM5 interface of mGlu2R as responsible for the formation of the heterodimer complex, and a simplified model was suggested which also involves the TM4-TM5 interface of 5-HT2AR in the formation of the 5-HT2AR/mGlu2R complex[34]). Interestingly, mGlu2R residues located at the intracellular end of TM4 were identified as necessary for the formation of a heterodimer with 5-HT2AR under physiological conditions, suggesting a possible role for the TM4-TM5-ICL2 interface in the mGlu2 receptor in a forming heterodimer complex with 5-HT2AR.[54] Similarly, regarding the effector coupling mechanism, we evaluated the implications of the differential ICL2 structural dynamics with respect to the site of interaction with G proteins. Thus, structural information from the crystal structure of the β2 adrenergic receptor/Gs protein complex,[42] the receptor’s intracellular cavity that directly interacts with the Gs protein, and the respective residues that constitute such cavity guided the identification of the equivalent cavity in 5-HT2AR. The differences in ICL2 orientation for HL and NHL complexes with the receptor are reflected in this region (e.g., the accessibility of residue S188ICL2 is larger in the HL systems when compared with the NHL counterparts, and this residue, together with S241C-term, is required in the process of agonist-mediated desensitization of 5-HT2AR[43]). These observations suggest a differential propensity for interaction with different intracellular effector proteins that ultimately will modulate the efficacy to activate different signaling pathways. These intrinsic differences could be further enhanced in the presence of different effector proteins that establish specific interactions with the intracellular part of the receptor. Interestingly, examples that demonstrate the importance of small conformational changes on protein function are known.[46,55] By using resonance energy transfer approaches, small values (∼2.5 Å) have been observed in the difference of distances at the intracellular segment of TM6 and the C-terminus for the arginine-vasopressin type 2 GPCR between conformations favored when Gs- and β-arrestin-biased agonists activate the receptor.[46] Moreover, based on enzyme activities and in receptor activation studies, it has been demonstrated that small conformational variations (as small as 1 Å) can cause significant changes in downstream function if amplified by the cell machinery.[55]

Concluding Remarks

Together, these findings provide what, to our knowledge, is the first indication of particular agonist-mediated conformational substates implicating specific conformations of ICL2 in systematically different responses of the 5-HT2AR to ligands with known differences in receptor-dependent phenotypes. The ligand-dependent behavior of the ICL2 may be characterized as a fine-tuning selective switch that depends on the intrinsic chemical and structural features of the GPCR ligands. As a specific example of the mode in which functional selectivity can be achieved, this study paves the way for further characterization of the heterogeneity of GPCR conformational states in the context of receptor functional selectivity/biased agonism.[4]

Methods

Microsecond unbiased all-atom MD simulations were carried out in the membrane-embedded 5-HT2AR in the unbound form (APO) and in complex with four different agonists: 5-HT, LSD, DOI, and LIS (Figure S1 in SI). For two of the systems described here (5-HT2AR/5-HT) and (5-HT2AR/LSD), shorter segments of the simulations (relative to the extent of the trajectories presented in this work) were part of a previous study from our group.[16] The respective overlap with the trajectories from our previous study is ∼ 0–175 ns for 5-HT2AR/5-HT and ∼0–250 ns for 5-HT2AR/LSD complexes (Figure S1B in SI). As a control, MD simulations of two closely related 1B (5-HT1BR) and 2B (5-HT2BR) human serotonin receptors in complex with the 5-HT ligand were carried out (100 ns each). All analyses were performed on the second half of the trajectories.

5-HT2AR Structure Complexes

The different systems were constructed as described previously.[16] Briefly, the 5-HT2AR model was created with homology modeling using as templates, the high-resolution X-ray crystal structures of the β2 adrenergic receptor (PDB accession code, 2RH1) and bovine rhodopsin (PDB accession code, 1U19).[56] As discussed in the Supporting Information, the crystal structures of two closely related human serotonin receptors, the 1B (5-HT1BR) and 2B (5-HT2BR) receptors, were solved after the MD simulations presented here were collected,[49,57] and thus, they were not considered as template for the 5-HT2AR structure but were used for validation and controls (see text in SI for detailed validation of the structure). The resulting 5-HT2AR structure is comprised of the segment S67 to K400 (a 28-residue segment in the long ICL3, the first 66 N-terminal residues, and the last 70 C-terminal residues were not included, see Figure S3A in SI) and was capped at its N- and C-termini by the acetyl and N-methylamide groups, respectively. A palmitoyl moiety was attached at position C397 based on the structural information on the β2 adrenergic receptor (PDB accession code, 2RH1). The protonation states of the amino acid residues were those more likely to be present at pH = 7.0, with the exception of positions: D1723.49 and E3186.30 which were protonated in the ligand-bound systems based on experimental considerations about GPCR activation but left charged in the APO form of the receptor.[16] All MD simulations started from the same 5-HT2AR structure, and the initial positioning of the agonists in the ligand binding pocket of 5-HT2AR was carried out by using several docking protocols (i.e., Autodock 4,[58] Simulated Annealing Docking,[59] Glide, and IFD (Schrödinger Inc.)), and were consistent with experimental information.[16] The 5-HT2AR systems were embedded in a physiologically relevant lipid membrane composed of a symmetric 7:7:6 mixture of SDPC (1-stearoyl-2-docosa-hexaenoyl-sn-glycero-3-phosphocholine):POPC (phosphatidylcholine):cholesterol, respectively. The GPCR-membrane systems were then hydrated by using the TIP3P[60] water model, followed by neutralization of the entire system by introducing ions to generate a NaCl salt concentration of 0.15 M.[16]

The parameters for the different ligands were obtained as described previously.[16] Briefly, quantum mechanical calculations were used to obtain an optimized structure and the electrostatic potential. These were subsequently used in the restrained electrostatic potential method (RESP) to generate the corresponding partial charges.[61] The topology and additional parameter files were prepared with Antechamber using the RESP charges and GAFF force field.[62]

All-Atom MD Simulations

Details of the 5-HT2AR simulations are as described previously.[16] Briefly, unbiased all-atom MD simulations were performed using NAMD[63] with the all-atom CHARMM27 force field with CMAP corrections for proteins and lipids[64] for trajectories of at least 1000 ns. Langevin dynamics and the hybrid Nosé–Hoover Langevin piston were used to maintain constant temperature (310 K) and constant pressure (1 atm), respectively.[65] Full electrostatics were evaluated using PME techniques with grid spacing <1.0 Å in each dimension and a fourth-order interpolation.[66] Bond lengths involving hydrogen atoms were constrained to their equilibrium values[64] by the SHAKE algorithm.[67] All MD simulations were performed with a 2.0 fs time step.

Control All-Atom MD Simulations in the 5-HT1BR/5-HT and 5-HT2BR/5-HT Complexes

The same protocol used in the simulations involving the 5-HT2AR systems with regards to setup conditions was used. The simulated systems consisted of the crystal structures of 5-HT1BR and 5-HT2BR (4IAR.pdb and 4IB4.pdb, respectively) complexed with the endogenous ligand 5-HT. To position the 5-HT ligand inside the binding pockets of the crystal structure of 5-HT1BR (PDB accession code, 4IAR) and 5-HT2BR (PDB accession code, 4IB4), we considered the equivalent atoms between ergotamine/dihydroergotamine and 5-HT. Since the ICL3 in both crystal structures are not resolved, the ICL3 loop was modeled so as to be consistent in length with the homology model of the 5-HT2AR. Other missing residues in the extracellular loops were modeled using the loop refinement algorithms in MODELLER. For 5-HT1BR, two additional residues were added at the C-terminus (C388 and T389) to achieve similar length as the 5-HT2AR model to introduce the conserved post-translational modification (palmitoylation). The coordinates for the palmitoyl chain were obtained from the crystal structure of the β2 adrenergic receptor (PDB accession code, 2RH1). In the case of the 5-HT2BR, one additional residue was added at the N-terminus (V47) for length consistency.

GPCR Identification Indexing

The Ballesteros and Weinstein indexing for GPCR residues[18] is adopted throughout the entire document as superscript accordingly. This generic numbering identifies particular residue positions in different GPCRs by assigning a pair of numbers A.B. In this notation, the first number corresponds to the transmembrane helix and the second number corresponds to the residue number relative to most conserved residue in transmembrane helix, which is assigned to 50.

Structural Alignment

For the structural analyses, all the structures were aligned to the structure of the β2 adrenergic receptor (PDB accession code, 2RH1) oriented with respect to the lipid bilayer according to the OPM database[26] by using the Cα atoms of the TM helices (for details about the segments considered in the structural alignment see Figure S3A in SI). Such alignment ensured that the Z-coordinate axis coincided with the helical axis of the TM bundle.

Overlap Coefficient

To quantify the similarity in the ICL2 conformations sampled from the different systems, the overlap coefficient was calculated.[27] After structural alignment, the center of mass of the ICL2 for the different conformations was calculated. Next, a grid of 1.0 × 1.0 × 1.0 Å was created; the different points, which represent the ICL2 center of masses, were sorted; and the frequency of each 1.0 Å3 cubes was obtained (ρ). The overlap coefficient (RAB) between the two sets of points, A and B, was calculated from eq 1. ρA and ρB are the frequencies of each of the sets for the “i” number of “cubes” that are compared. The RAB values range from 0 to 1, where 0 indicates no overlap between the two distributions and a value of 1 indicates a perfect overlap (identical distributions).

Principal Component Analysis

To quantify the major motions in ICL2, we used principal component analysis (PCA).[68] Using the Cα and heavy atom covariance matrices, we first found the first principal component (PC1) of the ICL2 movement in each system, which represented a large portion of the variance in all systems except for APO (see Table S1A in SI). To investigate differences in ICL2 dynamics in all five systems, we calculated PC1 for each simulation and then calculated the variance across that principal component for each other simulation. Atomic fluctuation correlations were calculated using carma,[69] and PCA was performed with in-house programs.[70]

Solvent Accessible Surface Area (SASA) Analysis

The SASA calculations were carried out with VMD.[71] The distributions for the four systems were integrated into six different aggregate distributions by averaging two component distributions: 5-HT + LIS, 5-HT + LSD, 5-HT + DOI, DOI + LSD, DOI + LIS, and LSD + LIS. From these, the three unique combinations were generated (5-HT + LIS versus DOI + LSD, 5-HT + LSD versus DOI + LIS (0.932), and 5-HT + DOI versus LSD + LIS), and the overlap between the different distributions was calculated by using eq 1.