Molecular study of different effects of serotonin 2A receptor modulators

The general activation mechanism of GPCRs has been studied extensively. On a macroscopic scale, an outward movement of Transmembrane Helix 6 (TM6) occurs upon activation, accompanied with an inward movement of TM3 and TM7 as shown in Fig 1. Note that the Ballesteros-Weinstein nomenclature is employed in further sections. [32] At the level of individual interactions between side chains (Fig 1A), the activation of GPCRs involves breaking the ionic lock of the E/DRY motif between residues Glu^6x30 of TM6 and Arg^3x50 of TM3. Additionally, the hydrophobic interaction between Ile^3x46 and Leu^6x37 in these helices also breaks. This allows TM7 to rearrange itself into a kinked helix so its NPxxY motif (Asn^7x49, Pro^7x50 and Tyr^7x53) interacts with a hydrophobic residue of TM3 (Ile^3x46). In turn, the side chain of the tryptophan "toggle switch" (Trp^6x48) rotates and TM6 bends at the level of the Phe^6x44 residue, forming an interaction between the Ile^3x40 and Phe^6x44 residues of the P^5x50I^3x40F^6x44 motif. Furthermore, agonists can create a hydrogen bond with Ser^5x46 resulting in an inward bulge in TM5. [33] Interestingly, this serine is substituted by an alanine residue in both rodent 5-HT_2A receptors and human 5-HT_2B/C receptors. [34] Therefore, this residue is thought to be involved in differences in ligand affinity through different species, sub-type selectivity binding, binding kinetics [35] and biased signaling. [36] It should be mentioned that these microswitches can be toggled independently, i.e., that ligands can alter the energetics for each switch differently. [37]

Throughout the GPCR family, the highly conserved tyrosines Tyr^5x58 and Tyr^7x53 (the latter is part of the NPxxY motif) forms a mediated hydrogen network which is believed to stabilize the active state. [38] More specifically for 5-HT_2A/B/C receptors, residue F329^6x41 stabilizes an alternative rotamer of Y254^5x58 through π-π stacking whereas this residue in other GPCRs contains a hydrophobic non-aromatic side-chain.[36]

Previous studies [33,39 –42] have shown that multiple combinations of "on" and "off" microswitches are likely relevant in the activation pathway and/or the differential response to different types of GPCR ligands. To survey the resulting landscape of states in the specific case of the serotonin-2AR, we performed MD simulations with 7 different ligands; zotepine, an antipsychotic antagonist, risperidone, an antipsychotic inverse agonist, IHCH-7086, a β-arrestin biased partial agonist, (R)-69, a G-protein biased partial agonist (both with retained hallucinogenic properties), lisuride, a G-protein biased partial agonist, and lastly, 25CN-NBOH and LSD, two β-arrestin biased partial agonist with hallucinogenic properties. These structures are displayed in Fig 1C. As it became clear during the study that an intracellular binding partner is necessary for the receptor to reach its active state, we performed each MD simulation in the absence and presence of an intracellular G protein construct. This allowed us to correlate the various pharmacological profiles of the selected ligands with the dynamics of the 5-HT_2AR by mapping their different effects on key microswitches as well as on large-scale helix movements during MD simulations.

The X-ray and cryo-EM structures listed in Table 1 were downloaded from the RCSB PDB databank (http://www.rcsb.org/pdb)↗. [43] and prepared for MD using CHARMM-GUI. [44]

Preparation X-ray coordinates. PDB IDs 6A93 and 6A94 represent a crystallographic dimer, from which only chain A was retained. In PDB IDs 6A93, 6A94, 7WC7, and 7WC9, the highly flexible intracellular loop 3 (ICL3) was replaced with the rigid apocytochrome b562RIL for crystallographic purposes. This fusion protein may affect the receptor's conformational equilibria and was therefore removed. [45] This resulted in a chain break between residues 265 (in TM5, the 5th transmembrane helix) and 313 (TM6). The CHARMM-GUI PDB manipulator [46 –48] was used to cap the N- and C-termini of the protein - including the ICL3 chain break - with acetyl and methylamide groups, respectively. During this step, the mutations S372N, S162K, and M164W were reversed. Disulfide bonds were defined between Cys148-Cys227 and Cys349-Cys353. CHARMM-GUI's loop modeling functionality was employed to model the unresolved residues 181–187 (ICL2) in 7WC7 and 7WC9. In addition, the PDB of 7WC9 included 6 monoolein molecules. Accordingly, a topology file for monoolein was constructed manually from CHARMM lipid force field [49] atom types; see "available data".

Preparation cryo-EM coordinates. The unresolved residues 100–104 (ICL1), 143–144 (ECL1), 224 (ECL2), and 394–401 (helix 8) in PDB IDs 6WHA and 9AS4 were modeled using CHARMM-GUI's loop modeling functionality. Conversely, TM1 and extracellular loop 3 (ECL3) were not fully resolved in the cryo-EM structures, but the modeling for these regions resulted in unrealistically strained structures. Instead, for 6WHA, the coordinates of those amino acids (TM1; 71–83 and ECL3; 348–353) were copied from 7WC9 after alignment of the surrounding secondary structure elements. When the structures of 8UWL, 7RAN, and 9AS4 became available, a production run of more than 1,500 ns using 6WHA was already in progress. RMSD time series were calculated for the 6WHA trajectory, using the 8UWL, 7RAN, and 9AS4 cryo-EM structures as references. To improve model accuracy, the frame with the lowest RMSD value for each system was selected as the source for filling in its missing coordinates of TM1 and ECL3 (following structural alignment). Disulfide bonds between Cys148–Cys227 and Cys349–Cys353 were again defined.

The N- and C-termini, including those at the ICL3 chain break, were capped with acetyl and methylamide groups, respectively. The cryo-EM structures also contained a mini-Gq protein [50] complexed with a single-chain variable fragment (scFv16), both of which were initially excluded during system preparation.

Simulation box. The "preprocessed" proteins resulting from the above procedure were solvated in a rectangular periodic membrane system using CHARMM-GUI [49,51 –56] The box dimensions and water layer thickness were selected to ensure that the minimum distance between the receptor and its periodic images was at least 24 Å. The receptor's orientation within the membrane was determined based on the Orientations of Proteins in Membranes (OPM) database. [57] Ultimately, the X and Y dimensions of the box were both set to 75 Å and the water thickness to 13.5 Å. Next, a heterogeneous bilipid membrane layer was selected which represents the biological environment of the receptor. Specifically, based on the lipid composition of synaptic vesicles, [58] the membrane was built with 57 POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine), 18 POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine), 8 PSM (1-palmitoyl-2-stearoyl-sn-glycero-3-phosphoethanolamine), 12 CER160 (ceramide 160), 52 POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine), 4 POPI (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoinositol) and 42 cholesterol molecules. The replacement method was used to build the membrane. After performing a ring penetration check, the system was solvated in a 0.15 M KCl solution using the Distance Monte Carlo algorithm implemented in CHARMM-GUI.

G protein construct. For simulations 9–15, a small construct representing the portions of the G protein that influence the receptor was assembled by modifying the PDB files of 6WHA and 7RAN. Specifically, residues 25–39, 69–86, and 222–246 of the mini-Gq protein in chain B were retained at their original positions. All N- and C-termini resulting from breaking the chains were capped with acetylamide and methylamide groups, respectively; note that this does not include the actual C-terminus of the chain as used in the cryo-EM structure (residue 246). The retained residues are shown in S1 Fig. To preserve the geometry of the original mini-Gq protein, distance restraints were applied to a selected set of atoms (listed in S1 Table along with their corresponding force constants). After positioning the resulting construct next to the receptor based on structure alignment, the proteins were solvated in rectangular periodic membrane systems as described above (but resulting in a slightly larger box). Time series of the Z-distance between the helical segment of the G protein that interacts with the receptor and the receptor itself are presented in S2 Fig. The angle between this helix and the xy-plane is also shown. These plots indicate that no decoupling of the G protein construct ever occurred. In addition, ICL2, which adopts a helical conformation in the active state of 5-HT_2A, remained helical throughout these simulations. By contrast, in simulations where the G protein construct was removed, ICL2 transitioned to a disordered loop and TM6 shifted toward its inactive conformation within the first few hundred nanoseconds after the simulation began. These observations suggest that the G protein construct helps stabilize key structural elements of the receptor's active state.

Ligands. The mol2 files of zotepine, risperidone, lisuride, IHCH-7086, (R)-69, 25CN-NBOH, and LSD were generated by building the ligands in Avogadro, followed by protonation at their most basic nitrogen atom using OpenBabel [59]. For all ligands, this resulted in protonation of the basic amine, which is important because experimental structures suggest a polar hydrogen bond with D155^3x32 of the receptor. Topology and parameter files for each ligand were automatically generated with CGenFF, based on the mol2 files, during the PDB manipulation step [51] Simulations 2b and 10 were initiated from the 6WHA template, with 25CN-NBOH removed. In simulations 9 and 13, 25CN-NBOH was replaced with zotepine and IHCH-7086, respectively.

The system was heated and simulated with NAMD 2.14 [60] using the standard protocol and simulation parameters from CHARMM-GUI. [49,61] Briefly, a 10000-step minimization was followed by three 250 ps equilibration steps and three more 500 ps equilibration steps, each at 303.15 K but with gradually weakening restraints. In the subsequent production phase, a time step of 2 fs was used in conjunction with restraints on the hydrogen atoms using the SHAKE algorithm. The pressure of the system was kept at 1 atm using NAMD's Nose-Hoover Langevin-piston barostat with a piston period of 50 fs and a piston decay time of 25 fs. A constant temperature was maintained using Langevin dynamics with a damping coefficient of 1 ps^-1. Non-bonded interactions were smoothly switched off over 10–12 Å using a switching function for the electrostatics and a force switching function for the van der Waals interactions (as recommended for using the CHARMM force field in NAMD). The long-range electrostatic interactions were calculated using the particle-mesh Ewald method with a mesh size of 1 Å and a 6th order spline interpolation.

Input files and trajectories were deposited on the Zenodo repository (see "Data availability statement".

Trajectory visualization and alignment were performed using VMD 1.9.3. [62] Measurements of the degrees of freedom (as defined in S2 Table) and RMSD calculations (see timeseries in S3-S4 Figs) were carried out using MDAnalysis 2.4.2. [63] Principal component analysis (PCA) was applied to the time series of the degrees of freedom (DOF) using the Scikit-learn Python package. [64] For each system, 3,000 frames were evenly sampled from the trajectory before performing PCA. To improve visualization, the first part of the simulation—where the system had not yet reached a stable state—was removed based on the RMSD data. S3 Table shows which frames were discarded for each system. Histograms, 2D probability density plots, and contour plots were generated using SciPy's kernel density estimator.[65] Bin widths were set to 0.2 Å and 1º for distances and (dihedral) angles, respectively. The 2D probability density plots from the different MD simulations were Boltzmann inverted but these populations densities should not be interpreted as converged free energy profiles; the logarithm merely makes the wide range of population densities more manageable as well as expressing them in more intuitive units of energy. For these plots, the initial portion of the trajectory was excluded to ensure that only equilibrated regions of the simulation were visualized; the discarded frames are as well listed in S3 Table.

Given the wide range of microscopic modulations of the 5-HT_2A receptor in response to various ligands, we attempt to systematize the changes in the aforementioned macro- and microswitches. Known relevant DOF in the receptor are listed in S2 Table, along with their definitions and measurement methods. To find patterns in the high-dimensional space of these DOF, it was subjected to PCA, retaining the first two principal components (PCs) for closer interpretation. [64] Most important DOFs are displayed in Fig 1.

When considering the PCA loadings plot (Fig 2A-B), the hydrophobic interaction (I^3x46 L^6x37), ionic lock and the attributes related to the movement of TM6 are all situated along the positive axis of in PC1, indicating that they constitute the major source of correlated attributes. In addition, the helix distances related to TM3 and TM7 appear along the negative PC1 axis, indicating an anticorrelation. This is in excellent agreement with the consensus activation mechanism wherein they play important roles, which suggests that PC1 distinguishes between inactive and active states. This is confirmed by the RMSD values (including those of the NPxxY motif) calculated relative to the inactive and active conformations, which correlate and anticorrelate with the aforementioned degrees of freedom, respectively.

The trajectories of the MD simulations listed in Table 1 are represented by the two PCs in Fig 2C, where the circles and triangles represent the receptor's X-ray and CryoEM state, respectively, while the regions of dots and crosses represent frames of different trajectories. It is readily apparent from the figure that the simulations without G protein are clustered on the left while the simulations with G protein are clustered on the right of the plot. This is a first indication that the presence of the G protein may be necessary for the receptor to reach its fully active conformation (as has been shown for β₂ adrenoreceptor by Dror and colleagues [40]) More evidence and further discussion regarding this possibility are presented below.

Intriguingly, when examining the G protein bound simulations in Fig 2C, it would appear that positive values of PC2 are associated with the activating/hallucinogenic character of 25CN-NBOH (cyan, simulation 14) and LSD (yellow, simulation 15), whereas the apo structure (green) and the antagonist zotepine (red) have the most negative PC2 values; lisuride, (R)-69, IHCH-7086 are situated in-between (brown, purple and blue, respectively). The DOFs related to the NPxxY motif are a major contributor to the variance along PC2. RMSD time series, calculated relative to their initial geometries for all trajectories, are provided in S3-S4 Figs for both the receptor and the NPxxY motif specifically. In the G protein–bound simulations (S4 Fig), the NPxxY motif maintains a low RMSD during the 25CN-NBOH and LDS simulations. For the non-hallucinogenic ligand IHCH7086, (like the psychedelics, a β-arrestin-biased agonist), an equilibrium appears to exist between low and high-RMSD conformations of the NPxxY motif. This may explain the two population distributions (blue crosses) observed in Fig 2C. During the simulations of the G protein–bound apo form (10) and with the G protein-biased ligands (R)-69 (12), and lisuride (11a), the NPxxY motif gradually transitions to a high-RMSD state within the first 200–300 ns. These observations might suggest that NPxxY is indicative of β-arrestin bias, or alternatively, that both a strong activating ligand and an intracellular binding partner are necessary to maintain the NPxxY motif in an active-like conformation (the distinction between the two possibilities is somewhat ambiguous in the present data set).

Time scales in MD simulations are generally too short to reliably observe the transition of 5-HT_2AR to an active state. However, studies have shown that an intracellular binding partner is necessary for at least some GPCRs to maintain their active conformation. [38,66] By removing the G protein, we can study the deactivation mechanism instead. In this context, Dror and colleagues [40] observed three distinct states during the MD simulation of β₂ adrenoreceptor in which an intracellular nanobody stabilizing the receptor's active state was removed. Specifically, a plot of the RMSD of the NPxxY motif of an inactive conformation versus the TM3-TM6 distance showed three distinct clusters of conformations. The first cluster is positioned around the crystal structure. A second cluster is formed after the rearrangement of the NPxxY motif near the end of TM7. From this state, a last cluster is obtained when the distance between TM3 and 6 decreases. To replicate the observations leading to this non-canonical activation pathway for the 5-HT_2AR, similar plots were constructed in S3 Fig for simulations 2b, 5 and 7, in which the mini-G_q protein was removed from their corresponding PDB templates. In S5 Fig, the black star on the 2D plot represents the cryo-EM structure. In the first row are the TM3–6 and the RMSD of the NPxxY values plotted, revealing two clusters: one with a high RMSD value and with low RMSD values. Overall, these plots suggest that the NPxxY motif initially shifts toward a more inactive conformation, followed by a decrease in the TM3–6 distance before it subsequently increases again. An exception occurs in simulations 5a and b with ligand (R)-69 after the removal of the G protein. This simulation was run twice and transitioned to an inactive or partially active conformation along two distinct pathways. This would seem to imply that the rearrangements of the helices do not necessarily have to happen in a specific order.

Fig 3 shows how the TM2–6 distance and the RMSD of the NPxxY motif relative to the inactive state evolve during simulations 2b, 5 and 7 and their corresponding G protein-coupled counterparts (10, 12 and 14) as reference. As foreshadowed by the PCA, it becomes evident that the presence of the G protein construct prevents inward movement of TM6. However, it does not seems to prevent the NPxxY motif to assume an inactive-like conformation for the apo and (R)-69 simulation, indicating the ligand's ability to modulate this structural motif.

The same analysis was also performed for two other degrees of freedom: the intracellular TM2–6 and TM3–7 distances, as defined in S2 Table. This pair of distances showed a large anticorrelation in the PCA (Fig 2) and was chosen as representative for the shape of the intracellular binding surface. Indeed, as shown in Fig 1B, helices 2, 3, 6 and 7 form a square delineating this binding surface. Combining the information from the different experimental structures, a global anticorrelated motion must occur along the perpendicular TM6–2 and TM3–7 axes during the transition from an inactive apo state to a fully active state containing the mini-Gq protein and an activating ligand. However, our MD simulations indicate that two motions tend to happen sequentially rather than in a concerted fashion, as can be seen in S5H-J Fig. These motions are explained in more detail in § 3.4.

The outward movement of TM6 is often associated with a bend in the helix near residues W336^6x48 and F332^6x44. As suggested by Wallach and colleagues [30] and Liao and colleagues [41], different rotamers of W336^6x48 may play a role in biased signal transduction. Their studies indicate that when the switch undergoes a further downward movement in the binding pocket, corresponding to smaller χ₂ values, it promotes β-arrestin transduction. More recently, this was experimentally confirmed by Gumpper and colleagues [67] with the β-arrestin biased ligand RS130–180. Fig 4 displays the χ₁ and χ₂-angles of these residues. Consistent with their studies, the 25CN-NBOH simulations (7 and 14; cyan) display a transition in the χ₂ angle of the tryptophan toggle switch between ~120° ("off" state) and ~50° ("on" state) (Fig 4A). Additionally, the presence of the G protein promotes the "on" state in both the apo simulation (8) and the 25CN-NBOH–bound simulations (11).

However, Fig 4A reveals an unusual transition of W336^6x48 in the (R)-69 bound simulation (5b). Instead of moving further down the binding pocket, the indole moiety of W336^6x48 rotates away so that χ₁ assumes values between 30–60º (plotted in Fig 4A, and side chains illustrated in Fig 4C). This specific indole conformation is only observed in simulations with the G protein biased ligand (R)-69 and in the receptor-G protein complex in combination with a potent agonist like 25CN-NBOH. Therefore, this specific state depicted in Fig 4C could be an interesting target for future docking studies aimed at identifying additional G protein-biased ligands with low efficacies.

Additionally, Fig 4B displays χ₁-χ₂-plots for F332^6x44, which is part of the highly conserved PIF motif. The χ₁-value indicate that the benzyl moiety can adopt a gauche-like and antiperiplanar-like conformations (Fig 4D). Notably, the antiperiplanar conformation is observed only in the presence of strong activating ligands (i.e., simulations 7, 8, 14 and 15) or intracellular binding of mini-Gq (as in simulation 10). While in this anti conformation, the PIF-motif assumes an "active"-like conformation, this effect might initiate a cascade of conformational changes throughout the receptor's intracellular side, biasing it toward binding an intracellular partner. The agonist (R)-69-bound system (simulation 5; purple) jumps on a very few occasions also to this anti state while its Trp-switch assumes this "alternative state. Consequently, the χ₂ angle distribution of F332^6x44 increases, as observed in simulations 10, 14 and 11, suggesting that the phenyl-moiety of F332^6x44 obtains a higher mobility. In turn, this might as well lead to intracellular signaling. Conversely, the antagonist zotepine penetrates much deeper into the binding pocket, blocking the tryptophan toggle switch and preventing interactions between the side chains in the connector region (PIF motif) (Fig 4E).

These observations might help explain the pharmacological differences between the (R)-69 and 25CN-NBOH agonists. Although they occupy roughly the same binding pocket, 25CN-NBOH possesses an additional aromatic ring which seems to be incompatible with the inactive state of the tryptophan switch, thus forcing it into its "on" state, whereas (R)-69 biases the tryptophan residue to an alternative "on" position from where its inactive state remains accessible (Fig 5A-5B).

For the N-benzylated phenethylamine, 25CN-NBOH, it was hypothesized that pushing the W336^6x48 toggle switch was essential for its potency and efficacy. [30,36,67] Indeed, Fig 5B shows that as the indole ring is pushed downward, a hydrogen bond forms between W336^6x46 and N376^7x49 of the NPxxY motif. This interaction stabilizes the active, "kinked" helix of TM7, which may explain the consistently active NPxxY conformation observed during simulations 14. This hydrogen bond is not observed in the other simulations as can be seen from the distance distribution between W336^6x46 and N376^7x49 in Fig 5C. However, while this H-bond might further stabilize the active NPxxY state, we could observe for the RMSD timeseries of IHCH-7086 and LSD (S4 Fig) that this feature is not strictly necessary to achieve an active-like NPxxY state.

When TM7 adopts its kinked helical structure, the tyrosine residue of the NPxxY motif (Y380^7x53) shifts inward, corresponding to an active conformation (Fig 5B). As also indicated by the PCA, a shorter distance between Y380^7x53 and Y254^5x58 (TyrTyr) is associated with a more active conformation. For the β₂ adrenergic receptor it was proposed that both tyrosines face inward in the active state, stabilizing it by forming a water mediated H-bond network. [33,38,40] However, in the active states of 5-HT₂ subfamily, the phenolic group of Y254^5x58 is displaced outward from TM6 and it interacts with F329^6x41 by means of π-stacking. This has been suggested to contribute to the preference of the receptor to couple with proteins of the Gq family. [36] During the (R)-69 simulation, an inward shift of Y254^5x58 was observed while Y380^7x53 was directed outward (Fig 5A).

To further systematize ligand–receptor interactions, a second PCA analysis was performed, excluding the Apo trajectories from the data set and extending the list of variables with five DOFs (S6 Fig) associated with the ligands. Since the latter do not belong to a congeneric series, the additional DOFs were simply calculated as the shortest distance between a polar non-hydrogen atom of the residue (shown in Fig 6A and 6C) and the nearest polar non-hydrogen atom of the ligand. In this new PCA, the first PC essentially expressed the same degree of freedom as before, i.e., the large-scale displacements associated with the activation mechanism continue to dominate the PCA. In sharp contrast, PC2 is now entirely different, containing important contributions from the TM5 bulge (S242^5x46), the ligand, and TyrTyr. Notably, (R)-69, lisuride and LSD, the agonists with a pronounced ability for G protein signaling, are the only ligands that form a H-bond with S242^5x46 during our simulations (see S7 Fig). This interaction leads to a greater distortion of TM5 (the TM5 bulge; see Fig 5A), which explains the shorter S242^5x46 – G369^7x41 distances compared to the other ligands shown in Fig 6B. The presence of the intracellular binding partner further reduces these distance distributions, indicating an even greater distortion of TM5 near residue S242^5x46.

Like in the agonist-bound β₂-ARs, the TM5 "bulge" induces a conformational change in the PIF motif. [33] When S242^5x46 shifts inward, it moves P246^5x50 toward I163^3x40. As a result, the isoleucine residue adopts a different inward rotamer, which subtly influences the position of F332^6x44. This alteration in the PIF region is thought to modify the energetic landscape of the receptor, facilitating the transition of TMs 6 and 7 from the inactive to the active state.

In contrast to the above, the β-arrestin–biased ligands IHCH-7086 and 25CN-NBOH do not form a H-bond with S242^5x46 (Fig 6C). As a result, I163^3x40 adopts a different rotamer. As hypothesized by Cao and colleagues, this difference may underlie their distinct agonist behavior; indeed, following this hypothesis, IHCH-7086 was intentionally designed to form no H-bond with S242^5x46. The TM5 bulge in simulations 6 and 13 (IHCH-7086) is less pronounced compared to the other simulations, which might explain β-arrestin–mediated signaling without detectable G_q activity. More generally, GPCR signaling bias appears to involve shared allosteric mechanisms, as proposed by Weis and Kobilka. [33] We speculate that the TM5 bulge may influence the PIF motif, thereby communicating with the intracellular side of the receptor to stabilize the active-like NPxxY and Tyr–Tyr motifs. Interestingly, 25CN-NBOH seems to bypass this mechanism by directly pushing the toggle switch W336⁶ˣ⁴⁸, enabling direct interaction with the NPxxY motif.

Since the outward motion of TM6 does not happen simultaneously with the inward motions of TM3 and TM7 during activation, different partially-active conformations can be proposed. To further investigate the populations of these different states, 2D probability density plots of TM6–2 versus TM3–7 are displayed in Fig 7.

When comparing the non- and G-bound systems, a notable shift in the density to higher TM26 distances is observed. This shift indicates that the G protein is essential for TM6 to achieve its full outward position, consistent with the PCA results and the discussion in §3.2. In addition, when looking at the plots for the hallucinogen bound receptors (25CN-NBOH; 7 and 14 and LSD; 8 and 14), it can be observed that the minimum is shifted slightly downward, meaning that TM3 and TM7 are shifted inward in response to the presence of the G protein. Combing these observations, it would indeed appear that the hallucinogen-G protein-bound complexes are located the farthest along an activation pathway consisting of an increase in TM2–6 and decrease in TM3–7 distance. Conversely, if we compare the location of this minimum with the other simulations, we can carefully assume that the other G protein-bound complexes are in a partially activated state, since either TM6 is shifted more inward, TM3 and TM7 shifted outward or both. The only exception to this trend is simulation 7, which includes both the antagonist zotepine and the G protein. When comparing the minimum of simulation 7 with those of the non-psychedelic compounds (simulations 10–13), we observe that upon docking zotepine into the G-protein-bound active conformation, the receptor fails to transition to a less-active state over the course of the simulation. This contrasts with other G-protein-bound simulations containing a partial agonist or no small-molecule ligand, where such transitions do occur (indicated by the grey arrows). While the absence of a transition with zotepine is most likely a shortcoming of the computational time constraints, an alternative explanation would be that zotepine hinders the conformational freedom of the receptor. More specifically, it seems possible that zotepine would increase the free energy of one or more intermediary conformations in the activation mechanism, thus inhibiting both deactivation and activation. Indeed, in the absence of the G protein there is a trend away from the upper-left corner of the graph with increasingly activated ligands, while the simulation with zotepine displays a compact distribution corresponding to a fully deactivated state. We further speculate that zotepine could accomplish this by sterically hindering the unraveling of a network of interactions that holds TM6 in place in both the activated and the deactivated conformation. This is supported by Fig 4E, where zotepine hinders the tryptophan switch and the interaction between F^6x44 and I^3x40 of the PIF motif.

In apo simulation (10), upon removal of the hallucinogen 25CN-NBOH, the population density shifts toward higher intracellular TM3–7 distances, again indicating that both the G protein and the hallucinogen are necessary for the receptor to remain in the conformational basin corresponding to the active state. In addition, we found that upon G protein activation, the psychedelics 25CN-NBOH and LSD populates a structural different state than the other agonists (Fig 7B). To facilitate comparison across simulations, Fig 8 presents a contour plot of the 50% most frequently sampled TM26 and TM37 distances.

Comparing the conformational landscapes of the psychedelics (simulation 14 and 15) with the other agonists (11–13) leads us to speculate that the ideal β-arrestin biased receptor conformation would feature a TM37 distance close to 15 Å and a TM26 distance around 18 Å (see Fig 8; blue, cyan, yellow). In this hypothesis, TM26 should be somewhat higher (around 21 Å) for G protein signaling, which would be potentiated by a further decrease in TM37. Combined with 25CN-NBOH's and LSD's high absolute binding affinity for the receptor, this would make it plausible that it causes hallucinations by inducing excessive populations of activated 5-HT_2A while at the same time being β-arrestin biased. It would also explain why both IHCH-7086 and (R)-69 are weak agonists while being β-arrestin and G protein biased, respectively.

In summary, after some equilibration, the simulations without G protein largely assumed an inactive state, especially along the TM6−2 axis – however, they showed a clear tendency in the direction of the active state for increasingly agonistic ligands, which is most pronounced along the TM3−7 axis as can be observed from the distinguishable regions in Fig 8. Similarly, the simulations with G protein largely remained in an active state (again, especially along the TM6−2 axis) but showed a clear tendency in the direction of the inactive state for increasingly inactive or antagonistic ligands, at least along the TM7−3 axis. The only exception to the latter observation is the simulation with zotepine, which biased the conformational ensemble toward an inactive state without G protein, but toward an active conformational ensemble with G protein, despite it being known as an antagonist. This is most likely a consequence of insufficient sampling, as discussed above. Apart from this exception, the global trend is consistent with the hypothesis that an activating ligand by itself does not make 5-HT_2A transition to its fully activated conformation, but rather adds a bias to the conformational energetics such that the binding of the G protein concurrent with the transition to an active state (presumably near the lower right of Fig 8) becomes energetically more favorable. Indeed, this mode of action has been proposed on the basis of computational and experimental data for the β₂-adrenergic receptor [33,40] and might be more widespread in the GPCR family. Going one step further, our results regarding the G protein versus β-arrestin-biased ligands are in excellent agreement with the observation by Weis and Kobilka that despite their differential pharmacology "these [different] agonists likely share allosteric transmission mechanisms". [33] Taken together, this leads us to hypothesize that hallucinations result from excessive populations of activated 5-HT_2A and that a non-hallucinogenic antidepressant effect could be achieved by a sufficiently mild activating ligand. This may imply an energetic effect that weakly favors the "fully activated" conformation (i.e., a quantitative effect), either by binding relatively weakly to the receptor (like the endogenous ligand 5-HT, when compared to hallucinogens [68,69]) or by binding strongly but having a modest activating influence of the conformational energetics. Alternatively, a mild activating ligand could induce a strong bias in the receptor's conformational energetics, but toward a "non-ideal" conformation that only weakly enhances the intracellular binding of the G protein- (which we would call a qualitative effect). Such a conformation might lie between the inactive and fully active conformation or correspond to the β-arrestin biased conformation (which traditionally would lead to desensitization). Either way, the above distinction between a "qualitative" and "quantitative" effect is difficult to make when studying 5-HT_2A in isolation because the conformational difference between the tentative G protein biased and β-arrestin biased conformations is modest, so that Hammond's postulate predicts that both effects will often occur together. This equally implies that this distinction may be of lesser concern in practical efforts to design non-hallucinogenic 5-HT_2A antidepressants.