Abstract From computational simulations of a serotonin 2A receptor (5-HT 2A R) model complexed with pharmacologically and structurally diverse ligands we identify different conformational states and dynamics adopted by the receptor bound to the full agonist 5-HT, the partial agonist LSD, and the inverse agonist Ketanserin. The results from the unbiased all-atom molecular dynamics (MD) simulations show that the three ligands affect differently the known GPCR activation elements including the toggle switch at W6.48, the changes in the ionic lock between E6.30 and R3.50 of the DRY motif in TM3, and the dynamics of the NPxxY motif in TM7. The computational results uncover a sequence of steps connecting these experimentally-identified elements of GPCR activation. The differences among the properties of the receptor molecule interacting with the ligands correlate with their distinct pharmacological properties. Combining these results with quantitative analysis of membrane deformation obtained with our new method (Mondal et al, Biophysical Journal 2011), we show that distinct conformational rearrangements produced by the three ligands also elicit different responses in the surrounding membrane. The differential reorganization of the receptor environment is reflected in (i)-the involvement of cholesterol in the activation of the 5-HT 2A R, and (ii)-different extents and patterns of membrane deformations. These findings are discussed in the context of their likely functional consequences and a predicted mechanism of ligand-specific GPCR oligomerization.

Author Summary The 5-HT 2A receptor for the neurotransmitter serotonin (5-HT) belongs to family A (rhodopsin-like) G-protein coupled receptors (GPCRs), one of the most important classes of membrane proteins that are targeted by an extensive and diverse collection of external stimuli. Recently we learned that different ligands targeting the same GPCR can elicit different biological responses, but the mechanisms remain unknown. We address this fundamental question for the serotonin 5-HT 2A receptor, because it is known to respond to the binding of structurally diverse ligands by producing similar stimuli in the cell, and to the binding of quite similar ligands with dramatically different responses. Molecular dynamics simulations of molecular models of the serotonin 5-HT 2A receptor in complex with pharmacologically distinct ligands show the dynamic rearrangements of the receptor molecule to be different for these ligands, and the nature and extents of the rearrangements reflect the known pharmacological properties of the ligands as full, partial or inverse activators of the receptor. The different rearrangements of the receptor molecule are shown to produce different rearrangements of the surrounding membrane, a remodeling of the environment that can have differential ligand-determined effects on receptor function and association in the cell's membrane.

Citation: Shan J, Khelashvili G, Mondal S, Mehler EL, Weinstein H (2012) Ligand-Dependent Conformations and Dynamics of the Serotonin 5-HT 2A Receptor Determine Its Activation and Membrane-Driven Oligomerization Properties. PLoS Comput Biol 8(4): e1002473. https://doi.org/10.1371/journal.pcbi.1002473 Editor: Dennis R. Livesay, UNC Charlotte, United States of America Received: November 7, 2011; Accepted: February 26, 2012; Published: April 19, 2012 Copyright: © 2012 Shan et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: The work was supported by the National Institute on Drug Abuse (Grants P01 DA-012923 and R01 DA015170). Computational support was provided by the National Science Foundation Terascale Computing System at the Texas Advanced Computing Center (TG-MCB090132 and TG-MCB080079N), the computer facilities at the Institute of Computational Biomedicine of Weill Cornell Medical College, and New York BlueGene at Brookhaven National Laboratory. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Discussion The MD simulations of the 5-HT-, LSD- and KET-bound 5-HT 2A R reported here provide the first molecular representation of the different effects that pharmacologically distinct ligands can have on the 5-HT 2A R. The concepts of “functional selectivity” [49], [50] and “receptor bias” [51] are frequently being used to explain the increasingly common observation of differential responses elicited by different ligands from the same receptor (e.g., for 5-HT 2A R see [4], [52]). However, no structural context had been identified for the distinct effects on the dynamics produced in the same GPCR by the binding of pharmacologically different ligands. Here we simulated the dynamics of the 5-HT 2A R binding of such pharmacologically distinct ligands, and identified different effects on the SM/FMs of the receptor. These effects were shown to lead to different rearrangements that correspond to the different levels of activation known to be produced by these ligands. Notably, the differential effects were shown to be consonant with the pharmacological characterization of the three ligands as a full, partial and inverse agonist, respectively. To our knowledge, such inferences were obtained for the first time here from unbiased atomic MD simulations, but they are in line with the increasingly detailed experimental evidence about ligand-related functional selectivity [49], [50], [51], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], with the proposals of ligand-selective conformations in the 5-HT 2A R [67] and the D 2 R [68], and with structural data indicating that GPCRs such as β 2 AR are stabilized in distinct conformational states by inverse, partial, or full agonists - respectively [12], [13]. In the current simulations, structural changes associated with SM/FM characteristics of an “activated state” of the 5-HT 2A R appear in sub-microsecond trajectories. In contrast, experimentally determined GPCR activation timescales generally vary from microseconds (photoactivation of rhodopsin [69]) to seconds (β 2 AR in living cells [70]). We emphasize that the conclusions reached here do not require the simulations to have converged to an “active state structure” of the kind claimed for the constructs determined crystallographically. Indeed, a number of modes of activation proposed from experiment share this characteristic and can also be significantly faster [71], [72], [73]. But in general, there are many reasons for the time scale differences between our results and functional measurements. In particular, the simulated system is an idealized construct in that all interaction components are placed in optimal positions to be at or near their targets. Titratable groups are also assigned their final charge states, e.g., when the D3.49 and E6.30 are in the protonated form in some of the constructs. Interestingly, the specific protonation form does not determine whether the ionic lock is formed or not (see Figures 3–4, and Figure S5 in Text S1); rather, the determinant factor is seen from our results to be the nature of the dynamics induced by the binding of a specific ligand. But considering that inactive GPCR (β 2 AR) may be pre-coupled to G-protein Gs [31] and the protonation of E3.49 in rhodopsin (an activation step) depends on transducin [74], the degree of precoupling will likely play a role in the activation time. Moreover, the simulation conditions (such as pH, salt, lipid composition, and crowding) certainly do not mimic completely those surrounding the receptor in living cells (e.g., it is known that the highly flexible DHA chain of SDPC, included in the lipid mixture used here, facilitates GPCR activation [75]), and similar time-scale differences have been observed between computer simulations and experiments for other GPCRs [76], [77]. The response of the membrane environment to the different ligand-induced structural re-arrangements produces a reorganization of the membrane around the receptor. This is reflected in (i)-the involvement of Chol in direct interactions with the protein [43], [78], that was shown here to affect the dynamics of the SM/FMs, and (ii)-the membrane deformations around the TM bundle of a GPCR [48], [79], described here with the use of the CTMD method [32]. Because the different ligand-determined conformational changes in 5-HT 2A R establish different patterns of local perturbations in membrane structure around the receptor complex, they were suggested to promote different ligand-dependent receptor oligomerization patterns through the hydrophobic mismatch between the TMs and the surrounding membrane [32]. This is supported by observations in the literature that: (i)- oligomeric associations of the dopamine D 2 R [27], 5-HT 2C R [28], and the β 2 AR [31] is ligand-sensitive; and (ii)- GPCR self-assembly is regulated by the mismatch between the hydrophobic length of the TM segments of GPCRs and the hydrophobic thickness of the lipid bilayer, as suggested by both experimental results [80] and computational studies for rhodopsin [32], [48], [79]. Along these lines, the results presented here suggest that the dimerization interfaces of 5-HT 2A R oligomers will be different for inverse agonist-, partial agonist-, or agonist-bound complexes, and moreover that the inverse agonist KET would promote more extensive 5-HT 2A R oligomerization than the full agonist (5-HT). We note that these experimentally testable predictions regarding possible oligomerization interfaces were obtained by analyzing monomeric GPCRs in complex with different ligands, without the need to simulate the dimers or higher oligomers.

Methods Construction of the simulated systems Several model systems of the serotonin 5-HT 2A R receptor were studied with all-atom MD simulations in explicit models of the hydrated lipid membrane environment. The 5-HT 2A R was simulated in complex with three ligands known to exhibit different pharmacological efficacies: the full agonist 5-HT, the partial agonist LSD, and the inverse agonist KET (Figure 2A). 5-HT 2A R constructs. For the simulation of 5-HT bound 5-HT 2A R, the protein was modified twice, very slightly, in regions distal to the binding site and the SM/FMs. The original receptor construct had a specific truncation of IL3 so that it consisted of 296 residues, from H1.28 to D5.57 and from R6.21 to K7.73, with an Ala–Ala linker between them (H1.28–D5.75–AA–R6.21–K7.73, where “–” denotes truncation). To match observations in crystallographic structures of several GPCRs [20], [25], [81], we thus added, at 112.5 ns, four residues to the IL3 (H1.28–L5.79–AA–S6.17–K7.73) in order to extend helical parts of TM5 and TM6, respectively, by two turns. The extension was done as follows: an average structure of the protein (including the ligand and palmitoyl derivative) was obtained from the trajectory between from 83.5 to 112.5 ns. The averaged structure was minimized first with constrained protein backbone and ligand heteroatoms followed by minimization without constraints. To enhance the flexibility of the truncated IL3, we extended the intracellular ends of TM5 and TM6 by 2 turns of helix each using Modeller [82], and selected the representative model by clustering the 100 models using either the extended TM5 or TM6. The loop between the extended TM5 and TM6 was refined using Modeller. The protein with extended TM5 and TM6 together with the ligand and palmitoyl chain was minimized first with protein backbone and 5-HT heteroatoms constrained, followed by complete minimization. The minimized complex was inserted in the lipid/water/ion environment from the snapshot at 112.5 ns to conserve interactions, after which the entire system was minimized and equilibrated with constraints on the protein backbone (velocities were reassigned in a random distribution based on the temperature). For the second extension, at 174.2 ns (i.e., 61.7 ns after extending TM5 and TM6), we added three more residues at the N-terminus (S1.25–L5.79–AA–S6.17–K7.73) to allow TM1 to reach beyond the lipid phosphate group region of the model membrane so as to avoid artificial interactions between the positive N-terminus and negative phosphate groups in membrane lipids. In addition, the N-terminus was acetylated and the C-terminus was N-Methylamidated to further avoid charge-charge interactions between termini and lipids. The simulations were then continued and the results reported here are from the 350 ns trajectory. Note that the initial homology model of 5-HT 2A R includes an artificially open “ionic lock” between residues R3.50 and E6.30 due to the use of the β 2 adrenergic receptor (β 2 AR) template in the homology modeling [83]. In the β 2 AR X-ray structure [45] the ionic lock is broken due to the co-crystallized lysozyme, but has been shown to consistently reform in MD simulations of inactive β 2 AR without the lysozyme [40]. The simulations of LSD-bound and KET-bound 5-HT 2A R, started from the same conformation as for the 5-HT bound 5-HT 2A R except that they included the extensions from the very beginning. In addition, to test whether KET, as an inverse agonist, is capable of reversing the conformation induced by the bound agonist 5-HT, we substituted 5-HT with KET in the activated 5-HT 2A R structure obtained at the end of the 5-HT simulation, and restarted that simulation with KET for an additional 500 ns (termed “KET-substituted simulation”). The protocol for this switch of ligand was as follows: (i)-An average structure (protein+5HT+palmitoyl chain) was generated using the last 50 ns of 5-HT simulations, and then minimized; (ii) 5HT was substituted by KET so that the docking pose of KET (Figure 2D, left panel) is aligned with the minimized average structure using backbone atoms of binding site residues: D3.32, S3.36, S5.42 and S5.46. The complex (protein+KET+palmityol chain) was minimized by fixing the heteroatoms of KET and constraining backbone atoms of the protein; (iii)-The minimized complex was combined with the lipid/water/ion environment from a snapshot at 350 ns of the 5-HT simulation, to conserve the interaction between the protein and the environment. Lipid/water/ion was minimized and then equilibrated. Finally the whole system was equilibrated by reducing constraints on protein backbone atoms and KET heteroatoms. Velocities were reassigned based on the temperature. Residues D3.49 and E6.30 were protonated in the 5-HT and LSD simulations (see also Discussion section, above), and deprotonated in the KET simulations (including KET-substituted simulation). We note that the protonation state of the E6.30 residue does not affect the state of the ionic lock, as we show in the separate simulation of KET-bound 5-HT 2A R where E6.30 residue is protonated (see Figure S5 in Text S1). In all simulations, C7.70 was palmitoylated by moving the coordinates of the palmitoyl chain (PALM) from PDB 2RH1 [45] onto the C7.70 of 5-HT 2A R. Loop structures determined from ab inito loop prediction. To enable full-scale 5-HT 2A R simulations, we refined the loops in 5-HT 2A R homology model described recently [83] using the Monte Carlo-Scaled Collective Variables ab initio method [84], [85]. For details see Methods and Table S2 in Text S1. Initial placement of the ligands. Protonated 5-HT, LSD and KET were docked into 5-HT 2A R using several docking protocols, including Autodock 4 [86], Simulated Annealing-Docking [87], Glide (Schrödinger Inc.), and IFD (Schrödinger Inc.). In Autodock, the GA-LS algorithm and a maximum number of 2.5×107 evaluations were used. Simulated Annealing-Docking was carried out following a protocol previously established in our lab [87], [88] starting from a binding pose of 5-HT predicted by Autodock and consistent with experimental data. Glide [89] was carried out with and without H-bond constraints on D3.32 and/or S5.46. Applying H-bond constraints on S5.46 generated more docking poses that were consistent with the experimental data. IFD [90] was carried out starting either from scratch or from Glide docking poses that were consistent with experimental data. Other docking parameters were set to default values. These procedures generated docking poses consistent with experimental data in the literature [2], [6], [7], [91], [92] (Figure 1A–D). In particular, for KET, IFD produced a cluster of docking poses in which the amines of the ligands interacted with D3.32 and S5.46, and its quinazoline ring immersed deep into the binding pocket close to W6.48. The binding site remained almost unchanged except that F6.52 rotated to avoid steric clashes with KET (Figure 1D). In this docking pose, which was used in the simulations, KET was in direct contact with the aromatic cluster (F5.47, F6.44, W6.48, F6.51 and F6.52) by forming an edge-to-face interaction with W6.48. Internal waters. X-ray structures of several GPCRs show water networks around the toggle switch W6.48 and the NPxxY motif [45], [93], [94], [95], and these are hypothesized to be important for receptor activation [96]. Internal waters were therefore introduced by solvating the 5-HT 2A R with Grand-Canonical Ensemble simulations using the Monte Carlo program MMC [97]. The procedure placed waters around W6.48 and the NPxxY motif consistent with the X-ray structures of cognate GPCRs. Lipid membrane composition and protein-membrane complex preparation. The 5-HT 2A R-ligand complexes were embedded in identical mixed and hydrated 7∶7∶6 1-stearoyl-2-docosa-hexaenoyl-sn-Glycero-3-phosphocholine (SDPC)/phosphatidylcholine (POPC)/Chol membranes. The choice of the lipid mixture was dictated by several factors: (i)-Chol is known to be important for modulating ligand binding, G-protein binding and activation of serotonin receptors [44], and has even been found in complex with the X-ray structures of amine GPCRs elucidated recently; (ii)-POPC represents a typical phospholipid component of the bilayer membrane, with one saturated and one mono-unsaturated acyl chain; and (iii)-SDPC lipid has been implicated specifically in the function of various GPCRs [98] and is abundant in neuronal tissues. In addition, the use of this lipid composition enables a comparison of Chol dynamics around 5-HT 2A R with observations from earlier MD studies of rhodopsin in somewhat different Chol-containing mixed membranes [43]. The lipid bilayer model was generated using VMD [99] to construct first a 120 Å×120 Å (in the x-y plane) hydrated POPC membrane patch consisting of 406 lipids; then, half of the POPC lipids were transformed to SDPC by translating corresponding atoms, i.e., from the POPC headgroups to PCGL, from the 16∶0 sn-1 chain to STEA, and from the 18∶1 sn-2 tail to DHA (missing atoms were built using internal coordinates in the all-atom CHARMM27 force field [100] with CHARMM31b1 [101]). To reduce steric clashes between POPC and SDPC molecules, we made use of the relatively straight DHA chains from the equilibrated SDPC membrane bilayer (http://www.lipid.wabash.edu/), and replaced all the DHA chains in the current membrane patch with the straight DHAs. The 5-HT-, LSD-, or KET-bound 5-HT 2A R were inserted into the lipid matrix by aligning the backbone of its seven most conserved residues (one in each TM, see [18]) with those of rhodopsin immersed in a POPC membrane [22]. Lipids within 0.8 Å of the protein and PALM were then removed leaving 354 lipids in total. 26 SDPC and 26 POPC in each leaflet were randomly replaced with Chol (PDB 3D4S [46]), by fitting Chol's C4, C5 and C6 atoms to STEA's C5, C6 and C7 or POPC's C35, C36 and C37. Chol positions were then refined by lateral translation to avoid clashes with other Chol, SDPC or POPC lipids. Finally, the systems were solvated with TIP3 water and 0.15 M NaCl salt. The final simulated systems consisted of 125 SDPC, 125 POPC, 104 Chol and 20–22K water molecules resulting in a total of 106–114K atoms. Force-fields and MD simulations The parameters for 5-HT were taken from an earlier study [7]. For LSD and KET, the results of geometry optimization and electrostatic potentials obtained from quantum mechanical calculations with the Gaussian program (Gaussian, Inc., Wallingford, CT) were used as input to the Restrained-ElectroStatic-Potential fit method [102] implemented in Antechamber [103] to derive charges. CHARMM topology and parameter files were then prepared with Antechamber using Restrained-ElectroStatic-Potential charges and GAFF force field. Force field parameter files for 5-HT, LSD and KET are included in Text S1. For protein, PALM, lipids etc., the all-atom CHARMM27 force field with CMAP corrections [100] was utilized (this approach is similar to a procedure used successfully in previous studies [104], [105]). All MD simulations were performed with the NAnoscale Molecular Dynamics (NAMD) suite [106]. As established in similar studies in the lab (e.g., see [107]), the simulations were conducted under constant temperature and pressure conditions with anisotropic pressure coupling, and utilized PME for long-range electrostatics [108]. The Nose-Hoover Langevin piston method [106] was used to control the target pressure with the LangevinPistonPeriod set to 100 fs and LangevinPistonDecay set to 50 fs. All MD simulations were performed with rigidBonds allowing 2 fs time step. All the simulated systems were equilibrated following a procedure described recently [109]. According to this protocol, the 5-HT 2A R backbones and the heavy atoms of the ligands were initially fixed and then harmonically constrained, and water was prevented from penetrating the protein-lipid interface. Constraints were released gradually in four 300 ps-step MD simulations with decreasing force constants of 1, 0.5, 0.1 and 0.01 kcal/(mol·Å2), respectively. Following this equilibration phase, all three GPCR-membrane complexes were simulated for 350 ns. The stability of the simulated complexes was monitored from the backbone RMSDs of the TMs in 5-HT 2A R using the following definitions for TMs: L1.29–L1.59, A2.38–Y2.67, L3.24–N3.56, S4.38–V4.62, D5.35–K5.67, N6.29–I6.60, G7.32–F7.56 and K7.58–I7.68. As illustrated in Figure 1E, after initial equilibration, the RMSDs in all the three systems were stable and fluctuated around or below 2 Å. In all three simulations the ligands maintained key interactions with the receptor (Figure 1B–E), consistent with previous experimental data [2], [6], [7], [91], [92]. Analysis of MD trajectories To quantify the changes in protein structure produced by the simulations we used various analysis tools. Analysis of protein structural data was carried out with Ptraj in AMBER 9 [110] and other tools discussed below. To quantify helix distortion parameters in the simulations, we used the Prokink package [111] implemented in Simulaid [112]. The correlation analysis on the time-dependent data of different variables, such as helix bend angles, face-shifts, as well as Chol-protein distances, was conducted following the procedure described in [43]. Briefly, the correlation analysis was carried out on two separate sets of dynamic variables. In the first, we followed the time-sequence of m = 8 selected variables that included proline kink and face-shift angles in TM6 and TM7, the minimum distances between the Chol at the EC end of TM6 and the residues on TM6 (I6.53, M6.57, I6.60, C6.61). In the second set, m = 12 dynamic variables were selected that included proline kink and face-shift angles in TM6 and TM7, the minimum distances between the Chol at the IC end of TM6–7 and the residues on TM6 and TM7 (K6.35, I6.39, F6.42, V6.46, L7.44, V7.48, V7.52, L7.55, F7.56). For each set, we first studied pair-wise correlations between different variables, and constructed the matrix of coefficients of determination, R2 (Figure 7D of the main text) using Spearman's rank correlation test (see for instance Ref. [113]). In this method, given N frames pairs of observations, (x i , y i ), first the x i and y i values separately are assigned a rank, and then the corresponding difference, d i between the x i and y i ranks is found for each pair. The R2 is then defined as: (1)Because it uses rankings, Spearman's method eliminates the sensitivity of the correlation test to the function linking the pairs of values and thus is preferred over parametric tests when no a priori knowledge exists on the functional relationship between x i and y i pairs. Combined Essential Dynamics (Comb-ED) analysis To compare the conformational spaces of 5-HT 2A R stabilized by the different ligands (i.e., 5-HT, LSD and KET), a Combined Essential Dynamics analysis [42], [114] was performed on C α atoms of the protein using Gromacs 3.3 [115]. Essential dynamics analysis separates the configurational space into an essential subspace with a few degrees of freedom which describe overall motions of the protein that are likely to be relevant to its function, and a physically constrained subspace describing local fluctuations. The method is based on the diagonalization of the covariance matrix of atomic fluctuations defined by: (2)where x i are the three Cartesian coordinates of the carbon atoms C α of the molecule under study, and the angular brackets denote averages over an ensemble of configurations and over the simulation time. The diagonalization of Eq. (3) yields eigenvectors that describe the directions of correlated positional changes in the molecule, whereas the eigenvalues indicate the total mean square fluctuation along these directions. In the Comb-ED, the covariance matrix is calculated for two or more concatenated trajectories, which are fitted on the same reference structure. Given this construct, the eigenvectors resulting from Comb-ED do not represent the essential degrees of motion of the molecules, but rather reveal differences and/or similarities in the dynamical and structural characteristics of the compared simulated structures. To identify structural differences between 5-HT 2A R stabilized by the three ligands, Comb-ED analysis was performed on 3 concatenated trajectories obtained by combining the trajectories for the pairs 5-HT-LSD, 5-HT-KET, and LSD-KET, each for the last 100 ns, respectively. Analysis of membrane deformations and the residual mismatch. The properties of the membranes were analyzed from the simulation trajectories using the recently described CTMD method [32]. Briefly, to quantify membrane deformations in the simulations and the hydrophobic mismatch energies, we calculated the time-averaged hydrophobic thickness profile of the membrane surrounding 5-HT 2A R in all trajectories and used solvent accessible surface area calculations to calculate the energy of the residual mismatch which exposes TM residues participating in unfavorable interfacial hydrophobic/hydrophilic interactions. To identify these residues, we determined if the TM is thicker or thinner than the surrounding membrane by comparing the hydrophobic thicknesses of the TM domains (using the following TM definitions (given in the Ballesteros-Weinstein generic numbering [18]): 1.29–1.59 (TM1), 2.38–2.67 (TM2), 3.24–3.53 (TM3), 4.39–4.63 (TM4), 5.38–5.63 (TM5), 6.33–6.59 (TM6), 7.30–7.56 (TM7)) to the local membrane thickness d memb calculated from the membrane sectors corresponding to each TM, as described in [32].

Supporting Information Text S1. Including supplemental methods, Figures S1, S2, S3, S4, S5, S6, Tables S1, S2, S3, Topology and parameter files for 5HT, LSD and KET. https://doi.org/10.1371/journal.pcbi.1002473.s001 (DOC)

Acknowledgments We thank Dr. Junmei Wang for his assistance with the use of Antechamber.

Author Contributions Conceived and designed the experiments: JS GK SM HW. Performed the experiments: JS SM. Analyzed the data: JS GK SM ELM HW. Contributed reagents/materials/analysis tools: SM ELM. Wrote the paper: JS GK SM ELM HW.