The pentameric ligand-gated ion channel 5-HT 3A exhibits a transient preactive state that provides valuable insight into its binding and activation mechanisms. In this preactive state, the allosteric regulation of the channel structure, and thereby function, upon ligand binding occur at timescales that are too fast to be measured experimentally. This study demonstrates the power of microsecond-timescale MD simulations in detecting these transitional preactive states upon ligand binding and describing their effects on channel function at the molecular scale. Such a mechanistic understanding of the channel function is a critical element in the design of therapeutics for the regulation of 5-HT 3A , which are needed to reverse the effects of numerous pathological conditions.

Aided by efforts to improve their speed and efficiency, molecular dynamics (MD) simulations provide an increasingly powerful tool to study the structure–function relationship of pentameric ligand-gated ion channels (pLGICs). However, accurate reporting of the channel state and observation of allosteric regulation by agonist binding with MD remains difficult due to the timescales necessary to equilibrate pLGICs from their artificial and crystalized conformation to a more native, membrane-bound conformation in silico. Here, we perform multiple all-atom MD simulations of the homomeric 5-hydroxytryptamine 3A (5-HT 3A ) serotonin receptor for 15 to 20 μs to demonstrate that such timescales are critical to observe the equilibration of a pLGIC from its crystalized conformation to a membrane-bound conformation. These timescales, which are an order of magnitude longer than any previous simulation of 5-HT 3A , allow us to observe the dynamic binding and unbinding of 5-hydroxytryptamine (5-HT) (i.e., serotonin) to the binding pocket located on the extracellular domain (ECD) and allosteric regulation of the transmembrane domain (TMD) from synergistic 5-HT binding. While these timescales are not long enough to observe complete activation of 5-HT 3A , the allosteric regulation of ion gating elements by 5-HT binding is indicative of a preactive state, which provides insight into molecular mechanisms that regulate channel activation from a resting state. This mechanistic insight, enabled by microsecond-timescale MD simulations, will allow a careful examination of the regulation of pLGICs at a molecular level, expanding our understanding of their function and elucidating key structural motifs that can be targeted for therapeutic regulation.

The homomeric 5-hydroxytryptamine 3A (5-HT 3A ) serotonin receptor is a pentameric ligand-gated ion channel (pLGIC) located at the postsynaptic cleft that converts chemical signals to electrical responses in the central and peripheral nervous system (1, 2). The primary chemical signal responsible for 5-HT 3A activation is the neurotransmitter 5-HT (i.e., serotonin). The binding of 5-HT causes conformational changes in the structure of 5-HT 3A that permit the flow of ions through the channel formed between its 5 monomer subunits, generating an action potential at the postsynaptic cleft (3). Clinically, pLGICs including 5-HT 3A regulate physiological functions such as nausea and are implicated in numerous psychiatric disorders including major depressive disorder, posttraumatic stress disorder, and Parkinson’s disease (4, 5). However, the mechanism by which agonist binding activates pLGICs and the structural basis that governs the transition between functional states is not well understood and remains a critical element hindering the design of therapeutics for many psychiatric disorders (6, 7). This mechanism is postulated to be more complex than a simple binary binding by activating ligands to a pLGIC, instead requiring the “priming” of a pLGIC by dynamic binding through transitional preactive states (8, 9), a mechanism that is further supported by the conclusions of this work.

Molecular dynamics (MD) simulations are a useful technique for examining the basis of pLGIC function, including the mechanisms that govern rapid and dynamic transitions between states that cannot be observed through experimental techniques (10⇓–12). Numerous pLGICs have been investigated using all-atom MD simulations, including the nicotinic acetylcholine (nAChR) (13), glutamate-gated chloride channel (GluCl) (14), and glycine (GlyR) (15) receptors, which have provided valuable insight into the structural response of pLGICs to agonist binding. MD simulations of 5-HT 3A have only been performed more recently, since its structure (excluding the intrinsically disordered intracellular domain) was first reported using X-ray crystallography in 2014 [Protein Data Bank (PDB) ID 4PIR (16)] and later through cryo-electron microscopy (cryo-EM) in 2018 [PDB ID 6BE1 (17)]. These apo structures, i.e., structures without 5-HT bound, have been reported as nonconductive (17⇓–19), while more recently additional structures of 5-HT 3A have been resolved with various agonists and antagonists and have been reported as both conductive and nonconductive (18, 20). Such states are initially assigned based on pore radius through the channel, where a minimum pore radius greater than that of a given wetted ion, such as Na+ or K+, represents a pLGIC in a conductive state because it is sufficiently wide to permit the translocation of ions.

However, assigning states to these static structures is not trivial due to limitations in structural resolution, symmetry assumptions, removal of highly flexible residues, and the addition of molecular components that may artificially constrain the protein in a given conformation, such as the nanobodies used in the original crystallization of 5-HT 3A (16). Therefore, MD simulations are performed to confirm the assignment of a state using careful analysis of structural dynamics including the pore radius profile, changes to the agonist-binding regions, and changes to secondary structure elements in the transmembrane domain (TMD) and extracellular domain (ECD). However, most simulations generally suffer from inadequately short timescales and usage of nonnative membrane lipids, and a lack of validation of critical simulation parameters including protein equilibration, membrane equilibration, and ligand representation in a force field that, if not properly validated, may bias a pLGIC in a nonnative conformation, subsequently yielding an inaccurate assignment of channel state.

In this work, we performed all-atom MD simulations of both apo and 5-HT–bound 5-HT 3A starting with PDB ID 4PIR (16) to determine the native protein state under conditions comparable with experiment (Fig. 1A). The study builds upon previous MD simulations of full-length 5-HT 3A (18⇓–20) in several notable ways. Firstly, each simulation was performed for 15 to 20 μs, which is an order of magnitude longer than any previous simulations of 5-HT 3A . We demonstrate such timescales are necessary to adequately equilibrate 5-HT 3A using the CHARMM36 (C36) force field (FF) for lipids (21, 22) and for proteins (23) in a lipid membrane and to observe the dynamic nature of 5-HT binding during channel priming, i.e., observation of a preactive state. Second, a ternary lipid membrane composed of 1-palmitoyl-2-oleoyl-SN-glycero-3-phosphocholine (POPC), polyunsaturated 1-stearoyl-2-docosahexaenoyl-sn-glyerco-3-phosphocholine (SDPC), and cholesterol was used to mimic the neuronal lipid environment in which 5-HT 3A resides (24), whereas past simulations used only POPC to represent the lipid environment (18⇓–20). The importance of lipid membrane type cannot be understated. It has been shown experimentally that cholesterol and polyunsaturated lipids interact directly with pLGICs such as nAChR (25⇓⇓⇓–29), which is largely homologous in structure to 5-HT 3A , to strongly regulate receptor structure and conductance. Lastly, we included several levels of simulation validation to ensure our model accurately represents experiment, including sufficient lipid membrane equilibration to ensure realistic lipid packing and careful parameterization of 5-HT for the CHARMM general force field (CGenFF) (30) using free energy perturbation (FEP) simulations to determine its solvation free energy, ensuring the accurate representation of 5-HT binding to 5-HT 3A . We specifically highlight these 2 parameters because insufficient lipid packing can result in a reduction of protein–lipid interactions around the TMD and because improper ligand parameterization can result in an unrealistic affinity of 5-HT to the binding pocket of 5-HT 3A . Together, these simulations provide a framework for how microsecond-timescale MD simulations must be used to examine the equilibration of pLGICs from a crystalized conformation and to examine allosteric regulation from ligand binding that can help reveal the nature and functional state of the ion channel.

Results

Accurately modeling the native structure of 5-HT 3A is an essential step to unraveling the structural elements that govern its function on a molecular level. To this end, we describe the following: the equilibration of 5-HT 3A from a crystalized conformation to a more native, membrane-bound conformation; the allosteric regulation of its TMD from dynamic binding of 5-HT to its principal binding pockets located on its ECD; and differences in systems with varied 5-HT concentration and lipid membrane composition. Table 1 defines the key parameters for these systems: 3 are composed of 5-HT 3A embedded in a POPC/SDPC/cholesterol membrane without 5-HT (apo), with ∼5 mM 5-HT (5 docked 5-HT), and with ∼15 mM 5-HT (5 docked 5-HT plus 10 5-HT added to the aqueous phase). A fourth system is composed of 5-HT 3A embedded in a POPC membrane with ∼5 mM 5-HT (5 docked 5-HT) to directly compare differences that arise from lipid membrane composition (Table 1). An example docking pose of 5-HT can be found in SI Appendix, Fig. S1.

Table 1. Simulation parameters for the 4 5-HT 3A systems simulated with membranes composed of POPC, SDPC, and CHOL or POPC only

Structural Overview. The assignment of ion channel state between conductive and nonconductive is generally determined through the pore radius profile, which demonstrates the pathway available for ion translocation through the channel (SI Appendix, Analysis Methods). In Fig. 2A, pore radius profiles are shown for the initial crystal structure, PDB ID 4PIR (16) (light blue, dashed) and for 5-HT 3A-Apo (orange) averaged over 20 μs after the 250 ns of equilibration from the crystal structure, where the SE of the pore radius is within the thickness of the trend line. PDB ID 4PIR is nonconductive because while the entire ECD is hydrated and accessible to ion translocation, the minimum pore radius is <2 Å through the TMD (gray), below the threshold needed to permit the passage of hydrated K+ ions. SI Appendix, Fig. S6 clearly shows the contrast between PDB ID 4PIR and an open conformation of 5-HT 3A [PDB ID 6HIN (20)], which would permit the translocation of ions with a slight movement of pore-lining helices and the L260 side chain. Fig. 2. (A) Pore radius profiles for the starting structure (light blue, dashed) and 5-HT 3A-Apo (orange) averaged over 20 μs with the transmembrane domain (TMD) shaded gray and error bars smaller than the thickness of the profile trend. (B) TMD snapshot for 5-HT 3A-Apo (including L260) shown as secondary structure and lines (respectively), with the initial structure as transparent white and the final structure colored by monomer (A, green; B, purple; C, pink; D, yellow; E, orange), with lipids, water, and ions removed for clarity. Representative helix labels for M1–M4 shown for monomer C with red, dashed lines connecting the centers of pore-lining M2 helix to demonstrate symmetry, and red, solid arrows indicating the principal direction of M2 fluctuation. (C) Interior angle θ for the 5 monomers in A vs. time. (D) Representative trend of current (i, picoamperes) vs. time (t, seconds) for the life cycle of a 5-HT 3A receptor with instantaneous currents shown for the resting (1), preactive (2), activated/open (3), and desensitized (4) states. The arrows indicate the possible states suggested by 5-HT 3A-Apo . (E) Pore radius profiles for the starting structure and 5-HT 3A-5mM (blue) averaged over 15 μs with the TMD shaded in gray and error bars smaller than the thickness of the profile trend. (F) TMD snapshot of 5-HT 3A-5mM depicted the same as C with cartoon 5-HT indicating 5-HT binding between monomers for the entire 15 μs (solid) and transient binding (transparent with dashed arrows). (G) Same as for C, but for 5-HT 3A-5mM . (H) Same as for D, but for 5-HT 3A-5mM . The minimum pore radius along the ECD of 5-HT 3A-Apo deviates by less than 1 Å from the crystal structure. This allows the ECD to remain hydrated and thereby accessible to ion translocation. However, the pore closes across the TMD (−50 Å < z < 0 Å) as demonstrated by the pore radius approaching a minimum value of 0 Å for 20 μs of simulation time (Fig. 2A, gray region). Furthermore, the snapshot of the TMD shown in Fig. 2B reveals that symmetry is maintained between the M2 helices of each monomer in this closed conformation, as depicted by the red dashed line connecting the centers of the pore-lining M2 helices and shown in the associated plot of the interior angle θ for each of the 5 monomers (Fig. 2C). Any fluctuations in the orientation appear cooperative and mostly in parallel to these lines of symmetry, as opposed to orthogonal to the center of the pore. Fig. 2D shows the life cycle of activation for 5-HT 3A and indicates that this conformation is most consistent with a nonconductive resting state; however, a postdesensitized state cannot be ruled out, which is also functionally nonconductive, but exists after the activation and desensitization of 5-HT 3A and is unresponsive to 5-HT binding. In Fig. 2E, pore radius profiles are shown for the initial crystal structure, PDB ID 4PIR (light blue, dashed) and for 5-HT 3A-5mM (blue) averaged over 15 μs after the 250 ns of equilibration from the crystal structure, where the SE of the pore radius is within the thickness of the trend line. The average minimum pore radius along the ECD of 5-HT 3A-5mM deviates by less than 1 Å from the crystal structure, but remains hydrated and accessible to ion translocation (Fig. 2E). Similarly to 5TH 3A-Apo , the pore radius across the TMD (−50 Å < z < 0 Å) approaches a minimum value of 0 Å, and the pore remains closed for 15 μs (Fig. 2E, gray region). Alone, this result suggests that the presence of bound 5-HT only impacts the ECD and not the TMD of 5-HT 3A , suggesting a lack of allostery and a desensitized state of 5-HT 3A . However, in Fig. 2F, a snapshot of the TMD of 5-HT 3A-5mM demonstrates that there is a significant antisymmetric shift between monomers D (yellow) and E (orange). Such an antisymmetric shift is depicted by a similar M2 helix center connecting red-dashed line and red arrows to demonstrate the direction of monomer shift conversely to the symmetric closure observed in 5-HT 3A-Apo and by the associated plot of the interior angle θ for each of the 5 monomers (Fig. 2G) that demonstrates significant deviations in interior angle symmetry. The shift begins around 5 μs and appears correlated to the departures in RMSD of both the ECD and TMD at around 5 μs (SI Appendix, Fig. S3 B and C). The shift dynamics can be summarized as follows: M2–D shifts outward from the center of the pore in a cooperative fashion with M1–D, M3–D, and M4–D (M2–C also displays an outward shift, but to a lesser extent), which would result in incremental pore opening if not for M2–E subsequently shifting inward toward the center of the pore in a cooperative fashion with M1–E, M3–E, and M4–E and occupying the space voided by M2–D, resulting in the final conformation shown in Fig. 2F. M2–A and M2–B display fluctuations parallel to the pore center similar to the M2 helices of 5-HT 3A-Apo . Shown schematically in Fig. 2F, 5-HT binds within binding pockets (bp) formed between pairs of monomers in the ECD (removed for clarity), for example bpCD is formed between the complementary monomer C and primary monomer D, while bpDE is formed between complementary monomer D and primary monomer E (an example of the initial binding pose in shown in SI Appendix, Fig. S1). The asymmetric shifts observed in 5-HT 3A-5mM appear correlated to 5-HT binding in these binding pockets. Monomer D is bound on either side by 2 5-HT and subsequently undergoes a conformational change as described in the previous paragraph. On the other hand, monomer C is only bound by 5-HT on its complementary face and subsequently displays a minor M2 conformational change, while monomer E is only bound by 5-HT on its primary face and subsequently displays significant M2 collapse into the channel. Furthermore, Fig. 2F shows that monomers A and B do not undergo observable conformational changes and do not appear to be allosterically regulated by 5-HT. We find that 5-HT only rebinds transiently to bpAB and bpBC and does not rebind substantially to bpEA, in contrast to the binding observed in bpCD and bpDE. Because the conformational changes caused by the binding of only 2 5-HT does not result in the channel opening sufficiently to permit the translocation of ions (9, 34), we cannot assign an activated state to the simulated conformation 5-HT 3A . However, because this 5-HT binding causes conformational changes to a single monomer of 5-HT 3A that resemble the initialization of channel opening, we conclude that 5-HT appears to be priming the receptor for activation, i.e., the conformation resembles a preactive state, as indicated in Fig. 2H, and not a desensitized or postdesensitized state. Symmetric pore closure is also observed in 5-HT 3A-5mM-POPC and in 5-HT 3A-15mM (SI Appendix, Fig. S7), similar to 5-HT 3A-Apo . In 5-HT 3A-5mM-POPC (SI Appendix, Fig. S7A), only one 5-HT is bound over the duration of the simulation and no large conformational change is observed in the TMD, suggesting that the binding of a single 5-HT is also insufficient to cause channel opening or even priming of the channel. However, monomer D of 5-HT 3A-15mM (SI Appendix, Fig. S7B, yellow) is bound on either side by 5-HT, similar to 5-HT 3A-5mM , but does not demonstrate a large outward shift in M2–D from the center of the pore, in contrast to the shift observed for M2–D of 5-HT 3A-5mM . Subsequently, we explored the nature of the binding events at each 5-HT concentration to determine differences in 5-HT binding, which result in different apparent outcomes with respect to monomer D between these 2 systems. Additionally, a full comparison of the TMD conformations for all 4 systems compared to the open conformation (PDB ID 6HIN) is shown in SI Appendix, Fig. S8.

Domain Interactions. The mechanism that governs the transition between states of 5-HT 3A lies in the allostery between 5-HT binding and secondary structure elements in the TMD and ECD, namely the M2–M3 loop (TMD), the β1–β2 loop (ECD), the β8–β9 loop (ECD), and the pre-M1 region (TMD/ECD) (6, 18, 33), which are shown for 5-HT 3A-5mM monomers C, D, and E in Fig. 4A. The β1–β2 loop, the β8–β9 loop, and the pre-M1 region are all directly connected to residues involved in 5-HT binding through β-strands and all interact either directly or indirectly with the M2–M3 loop. The M2–M3 loop is critical to channel gating because it acts as a spring that regulates the orientation M2, which is critical because separation between neighboring M2 helices is needed to allow the turning of L260 away from the center of the pore [shown in an open conformation PDB ID 6HIN (20) as an indication of a conductive state; SI Appendix, Fig. S6], and outward tilt and/or translation of M2 allows this separation of helices to occur. The M2–M3 loop has been reported as natively restrained in the apo form (18), which we also observe in this work, meaning that pore opening is an active process that requires the release of the loop from a native restraint. Fig. 4. (A) Principal shifts in secondary structure motifs at the transmembrane domain (TMD)/extra cellular domain (ECD) interface for 5-HT 3A-5mM due to sustained binding events in bpCD and bpDE. Monomers C (pink), D (yellow), and E (orange) shown as secondary structure with initial structures colored as transparent white. Monomers A and B, and all M4 helices removed for clarity. Distances between the center of masses of the β1–β2 loops colored by monomers for 5-HT 3A-Apo (B) and 5-HT 3A-5mM (C). The disparity between the conformational change observed in monomer D and the collapse observed in monomer E initiates with the M2–M3 loop, which in monomer D coils into a helix that causes M2–D to tilt outward along with M3–D (Fig. 4A). The opposite occurs in the M2–M3 loop of monomer E, which becomes more elongated, allowing M2–E to tilt inward along with M3–E. Coiling of M2–M3–D is permitted because of a shift in β1–β2–D away from M2–M3–D (and toward β1–β2–E), while β1–β2–E and β1–β2–C do not demonstrate significant shifts (Fig. 4A). Distances between the center of masses of the 5 β1–β2 loops are shown in Fig. 4B for 5-HT 3A-5mM to demonstrate the deviation from the average distance between β1–β2 loops and the distances between β1–β2–D and its neighbors and in Fig. 4C for 5-HT 3A-Apo , which serves as a control case (no 5-HT is bound). We propose that, of the β1–β2 loops, only β1–β2–D shifts because bpDE residues D–Y46, D–W63, and D–R65 all lie on β-strands 1 and 2 of monomer D, which at the onset of binding are pulled toward the primary monomer E through 5-HT in bpDE as the intermediary, while the pentameter D residues of bpCD are pulled toward the complementary monomer C through 5-HT in bpCD as the intermediary, causing a twisting of monomer D, which includes twisting of β1–β2–D. Neither monomers C or E are bound on either side by 5-HT, so they do not display the same twisting as seen in monomer D, but rather are constrained in only one direction by 5-HT binding. Moreover, this twist of monomer D and coiled orientation of M2–M3–D is stabilized by a salt bridge that forms between D271 of M2 and K54 of β1–β2–D, which locks monomer M2–D in a conformation that reduces channel occlusion (SI Appendix, Fig. S12A). Alternatively, because β1–β2–E does not exhibit a twisting motion, K54 of monomer E is stabilized by forming a salt bridge with D52 and E53 of β1–β2–D when β1–β2–D twists toward β1–β2–E (Fig. 4A). The twisting of β1–β2–D further allows D271 of M2–E to also form a salt bridge with K54 of monomer D, but unlike D271 of M2–D, this salt bridge constrains M2–E in an orientation that creates channel occlusion instead of decreasing restriction through the channel (SI Appendix, Fig. S12B). Additionally, differences in shifts observed in the β8–β9 loops of monomers C and D appear to play a role in stabilizing the conformation of monomer D that reduces channel occlusions and the conformation of monomer E that increases it. The M2–M3 loop is natively constrained by hydrogen bonding to the β8–β9 loop and pre-M1 region of the complementary monomer, constraints that are broken, but then reconfigured between T276 of M2–M3–D and Q184 of β8–β9–C in the form of hydrogen bonding (SI Appendix, Fig. S12C). However, due to an upward shift in β8–β9–D as a result of twisting observed in monomer D, T276 of M2–M3–E and E186 of β8–β9–D also form a hydrogen bond (SI Appendix, Fig. S12D), but much like the D271–K54 salt bridge, this hydrogen bond contributes to lock M2–M3–D in a conformation that reduces channel occlusion and M2–M3–E in a conformation that increases it. In both cases, a salt bridge is formed between E186 and R218 (SI Appendix, Fig. S12 C and D), which strengthens the configurations governed by hydrogen bonding with T276. In summary, bound only to the primary face of 5-HT 3A-5mM monomer C, 5-HT appears to regulate the β8–β9 loop, which is connected through β-strand linkage to the primary face, but not to the β1–β2 loop, which is only connected to the complementary face (Fig. 4 and SI Appendix, Fig. S12). Conversely, bound only to the complementary face of 5-HT 3A-5mM monomer E, 5-HT appears to regulate the β1–β9 loop, which is connected through β-strand linkage to the complementary face, but not to the β8–β9 loop, which is only connected to the primary face, preventing the M2–M3–E loop from expanding outward, and allowing M2–E to collapse inward (Fig. 4 and SI Appendix, Fig. S12). Only in the case of 5-HT 3A-5mM monomer D, do we report full allosteric regulation and activation of a 5-HT 3A monomer in contrast to its neighboring monomers, which are only partially bound by 5-HT. Significant shifts of these loops in the other monomers of 5-HT 3A-5mM and in the 3 other systems were not observed, indicating that the binding of 5-HT that allosterically regulates monomer D of this system is indicative of preactivation, whereas at least 3 bound 5-HT are required to achieve complete activation (9, 34). Lastly, changes in the TMD are not solely related to agonist binding in pLGICs and can also be regulated by antagonists, glycosylation, phosphorylation, and the lipid membrane itself. This last effect of lipid interactions is relevant to the conformational transition of monomer D in 5-HT 3A-5mM . Therefore, we subsequently discuss the differences between the 3 systems with a membrane composed of a 7:7:6 mixture of POPC/SDPC/cholesterol and the final system with a membrane composed of pure POPC, as defined in Table 1.