Structural basis of ligand recognition at the human MT1 melatonin receptor.

Melatonin (N-acetyl-5-methoxytryptamine) is a neurohormone that maintains circadian rhythms1 by synchronization to environmental cues and is involved in diverse physiological processes2 such as the regulation of blood pressure and core body temperature, oncogenesis, and immune function3. Melatonin is formed in the pineal gland in a light-regulated manner4 by enzymatic conversion from 5-hydroxytryptamine (5-HT or serotonin), and modulates sleep and wakefulness5 by activating two high-affinity G-protein-coupled receptors, type 1A (MT1) and type 1B (MT2)3,6. Shift work, travel, and ubiquitous artificial lighting can disrupt natural circadian rhythms; as a result, sleep disorders affect a substantial population in modern society and pose a considerable economic burden7. Over-the-counter melatonin is widely used to alleviate jet lag and as a safer alternative to benzodiazepines and other sleeping aids8,9, and is one of the most popular supplements in the United States10. Here, we present high-resolution room-temperature X-ray free electron laser (XFEL) structures of MT1 in complex with four agonists: the insomnia drug ramelteon11, two melatonin analogues, and the mixed melatonin-serotonin antidepressant agomelatine12,13. The structure of MT2 is described in an accompanying paper14. Although the MT1 and 5-HT receptors have similar endogenous ligands, and agomelatine acts on both receptors, the receptors differ markedly in the structure and composition of their ligand pockets; in MT1, access to the ligand pocket is tightly sealed from solvent by extracellular loop 2, leaving only a narrow channel between transmembrane helices IV and V that connects it to the lipid bilayer. The binding site is extremely compact, and ligands interact with MT1 mainly by strong aromatic stacking with Phe179 and auxiliary hydrogen bonds with Asn162 and Gln181. Our structures provide an unexpected example of atypical ligand entry for a non-lipid receptor, lay the molecular foundation of ligand recognition by melatonin receptors, and will facilitate the design of future tool compounds and therapeutic agents, while their comparison to 5-HT receptors yields insights into the evolution and polypharmacology of G-protein-coupled receptors.

The pleiotropic functions of 5-HT are in part mediated through 12 different subtypes of G protein-coupled receptors (GPCRs), which are regulated by a variety of factors, including ligands and lipid molecules[5]. Phospholipids and cholesterol are known to bind GPCRs to modulate their activity[6,7]. Three major phosphoinositides are phosphatidylinositol (PI), phosphatidylinositol 4-phosphate (PI4P), and phosphatidylinositol 4,5-bisphosphate [PI(4,5)P2, or PIP2][8]. PI4P is the precursor of PIP2 and was reported as a major mediator of GPCR-stimulated production of the second messenger, diacylglycerol (DAG)[9]. In addition, various types of phosphoinositides have been reported to bind GPCRs and regulate their activities[10,11], but the lack of structural information limits the mechanistic understanding of the exact types of phosphoinositides and how they regulate GPCRs.

Many GPCRs also exhibit intrinsic basal activity in the apo state[12,13], including 5-HT receptors[14], whose high basal activities are required for their normal physiological functions and inhibition of these basal activities by inverse agonists causes adverse effects[15]. 5-HT1A, 5-HT1D, and 5-HT1e are among the 5-HT receptors that are primarily coupled to Gi/o. 5-HT1A is the prototype serotonin receptor that exhibits high basal activity and is subjected to lipid regulation[14,16]. 5-HT1A ligands, including aripiprazole, have been used to treat many CNS diseases such as schizophrenia and depression. 5-HT1D exhibits high sequence homology with 5-HT1B and they are the primary targets of the triptan class of antimigraine drugs[17]. 5-HT1e is expressed in brain regions important for memory[18]. No structure of 5-HT1A, 5-HT1D, and 5-HT1e has been reported. In this paper we report five cryo-EM structures of serotonin receptor-Gi complexes, including three 5-HT1A structures, in the apo state or bound to 5-HT and the aripiprazole; one 5-HT1D structure bound to 5-HT, and one 5-HT1e structure bound a 5-HT1e- and 5-HT1F- selective agonist BRL-54443. These results address a long-standing question of phospholipid regulation of GPCRs and reveal the basis for 5-HT pan-agonism and drug recognition at the serotonin receptor system.

Structures of 5-HT1A, 5-HT1D, and 5-HT1e

We assembled serotonin receptor-Gi complexes by co-expression of the receptors with a dominant-negative human Gαi1 [19] and human Gβγ (Extended Data Fig. 1a, b, 2a, b). The structures of the 5-HT1A-Gi complexes in the apo state, or bound to 5-HT and aripiprazole, were determined at the resolution ranges of 3.0-3.1 Å. The structures of the serotonin-bound 5-HT1D-Gi complex and the BRL-54443-bound 5-HT1e-Gi complex were both determined at 2.9 Å (Fig. 1, Extended Data Fig. 1-2). The EM maps are sufficiently clear to define the position of the receptor, the G protein trimer, and the bound ligands in the receptor-G protein complexes (Fig. 1, Extended Data Fig. 1-2). The overall structures of 5-HT1A, 5-HT1D, and 5-HT1e consist of a canonical transmembrane domain (TMD) of seven transmembrane helices (TM1-7), a short intracellular loop 2 (ICL2) helix, and an amphipathic helix H8 (Fig. 1b).

Extended Data Fig. 1

Sample preparation and cryo-EM of the 5-HT1A–Gi complexes.

a, Analytical size-exclusion chromatography of the purified complex. b, SDS-PAGE/Coomassie blue stain of the purified complex. Experiments were repeated three times with similar results. c, Representative cryo-EM image (scale bar: 30 nm) from 4,179 movies and 2D averages (scale bar: 5 nm) of 5-HT–5-HT1A–Gi complex. Experiments were repeated three times with similar results. d-f, Flowchart of cryo-EM data analysis and the local resolution shown for the density of apo- (d), serotonin bound- (e), and aripiprazole bound- (f) 5-HT1A–Gi complexes. g, ‘Gold-standard’ Fourier shell correlation (FSC) curves. h, the local resolution shown for the density of water molecules (W1-W4) in the ligand binding pocket of the apo 5-HT1A–Gi structure.

Extended Data Fig. 2

Sample preparation and cryo-EM of the 5-HT1D–Gi–scFv16 and the 5-HT1e–Gi–scFv16 complexes.

a, b, Analytical size-exclusion chromatography and SDS-PAGE/Coomassie blue stain of the purified 5-HT1D–Gi–scFv16 complex (a) and the 5-HT1e–Gi–scFv16 complex (b). Experiments were repeated three times with similar results. c, Representative cryo-EM image (scale bar: 30 nm) from 4,375 movies and 2D averages (scale bar: 5 nm) of 5-HT1D–Gi–scFv16 complex. d, Representative cryo-EM image (scale bar: 30 nm) from 5,249 movies and 2D averages (scale bar: 5 nm) of 5-HT1e–Gi–scFv16 complex. e, f, Flowchart of cryo-EM data analysis, the local resolution of the density, and the ‘Gold-standard’ Fourier shell correlation (FSC) curves of the 5-HT1D–Gi–scFv16 complex (e) and the 5-HT1e–Gi–scFv16 complex (f).

Fig. 1

Cryo-EM structures of the 5-HT1A–Gi complex, the 5-HT1D–Gi complex, and 5-HT1e–Gi complex.

a&b, Cryo-EM maps and structural models of the Gi complex of 5-HT1A, 5-HT1D, and 5-HT1e. Ligands are shown at right side of the complexes with surrounding density maps. c, Comparison of the ICL3 of 5-HT1B with 5-HT1A, 5-HT1D, and 5-HT1e.

All three 5-HT receptors share a similar conformation, including a pronounced outward movement of TM6 at the cytoplasmic end relative to the inactive 5-HT receptor structures[20], a hallmark of GPCR activation. Although the overall structures of the three 5-HT receptor-Gi complexes are similar to that of the 5-HT1B-Go complex[21], several distinct features were observed. Compared to 5-HT1B, 5-HT1A and 5-HT1D have helical extensions at the cytoplasmic end of TM6 (13 residues for 5-HT1A and 9 residues for 5-HT1D, Fig. 1c). 5-HT1D has an additional 7-residue extension at the cytoplasmic end of TM5 as part of intracellular loop 3 (ICL3, Fig. 1c). These additional structural elements in 5-HT1A, 5-HT1D, and 5-HT1e are close to the Gαi-Ras domain and form extra interactions in the receptor-G protein complex assembly (Fig. 1c).

PI4P regulates the activation of 5-HT1A

The most surprising observation in the three 5-HT1A-Gi structures is a bound phospholipid at the interface between the receptor and the Gαi (Fig. 2, Extended Data Fig. 3a-e). The structures of 5-HT bound, aripiprazole-bound, and the apo-state 5HT1A present a clear EM map that can be best fit with a PI4P molecule (Extended Data Fig. 3f). The two acyl chains with 14-16 visible carbons form extensive hydrophobic interactions with TM6 and TM7/H8 (Fig. 2a, b). The PI4P head group is inserted into a cavity formed between TM3/6/7 of the receptor and the α5 helix of Gi-protein (Fig. 2a, b, Extended Data Fig. 3d, e). The 4-phosphate group of PI4P forms a salt bridge with the conserved R1343x50 of the receptor and a hydrogen bond with C351 of the Gαi (Fig. 2b, Extended Data Fig. 3e). The myo-inositol group interacts with the receptor residues T3466x36, K3456x35, F4037x56, N4048x47, K4058x48, and with G352 of the Gαi protein (Extended Data Fig. 3e). We reasoned that the binding of PI4P at the interface between 5-HT1A and G protein might stabilize their complex formation. To test this hypothesis, we reconstitute 5-HT1A-Gi (WT) complex to peptidiscs[22] with different lipids, in which the composition of lipids was defined. Then the GTPase-Glo assay was adopted to determine the effect of phospholipids on GPCR-mediated G protein activation, by measuring the GTP hydrolysis activity of Gi when coupled to 5-HT1A in the presence or absence of PI4P[10]. We found that GTP hydrolysis was enhanced by 2.4-fold in the presence of PI4P (Fig. 2c), suggesting that PI4P improves G protein coupling and GTPase activity. PI and PIP2 also enhanced 5-HT1A-mediated GTP hydrolysis but to a lesser extent (Fig. 2c). Phosphatidylethanolamine (PE) and phosphatidylcholine (PC), which are the major membrane phospholipids, as well as phosphatidylglycerol (PG) and phosphatidylserine (PS), also had modest effects on GTPase activity but significantly lower than PI4P, demonstrating the specific role of PI4P in enhancing 5-HT1A activity (Fig. 2c). On the other hand, PI4P elevated the 5-HT1A basal activity and increased the efficacy of 5-HT1A-mediated recruitment of Gi-protein, suggesting that PI4P can act as a positive allosteric modulator of 5-HT1A (Fig. 2d). Mutations on the PI4P binding residues R1343x50A, K3456x35A, K4058x48A of 5-HT1A reduced the G-protein activation and abolished the regulatory function of PI4P (Fig. 2e).

Fig. 2

Regulation of 5-HT1A by PI4P and cholesterol.

a, The bound PI4P at the interface between 5-HT1A and Gi protein. PI4P is shown in red sticks with a surrounding density map at threshold 0.02. b, Interactions of PI4P at the 5-HT1A–Gi interface with interacting residues shown in sticks. Hydrogen bonds are shown with dashed lines. c, Receptor activation as measured by GTPase-Glo assay for different lipid-containing 5-HT1A–Gi peptidisc. Lower levels of residual GTP indicates higher levels of G protein activity upon receptor-mediated GDP/GTP exchange. The activation of 5-HT1A in different lipid-containing peptidisc was compared with that in LMNG micelle in the presence or absence of 5-HT. d, PI4P regulates 5-HT1A as a positive allosteric modulator. The G protein recruitment signal was detected by NanoBiT assay. Data are presented as mean ± SEM of three independent experiments performed in triplicate. e, Mutations on PI4P binding residues R134A, K345A, and K405A reduced the 5-HT1A-mediated G-protein activation and abolished the regulation of PI4P. Data in panels c-e are presented as mean ± SD of three independent experiments performed in technical triplicate. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, no significance; two-tailed paired t-tests. f, Cryo-EM map of the apo 5-HT1A and the surrounding lipids. The densities of the 5-HT1A (green), cholesterol or fatty acids (purple), and PI4P (red) are colored with a 0.027 threshold.

Extended Data Fig. 3

Lipids regulation in 5-HT1A receptor.

a-c, The EM map of apo–5-HT1A–Gi complex and the surrounding lipids are shown with different threshold of 0.025 (a), 0.03 (b), and 0.04 (c). d, Interactions of PI4P at the 5-HT1A–Gi interface. e, Interaction of the PI4P head group with the TM6/TM7/Gαi pocket. Hydrogen bonds are shown with dashed lines. f. Comparison of the density fitting for PI4P, PI, and PIP2. The place of density that is not fit well is circled by dash line. g. 5-HT1A-mediated Gi activity are regulated by PI, PI4P and PIP2 with the greatest degree of PI4P regulation. GTPase-Glo assay was perform in LNMG buffer. Lower levels of residual GTP indicates higher levels of G protein activity upon receptor-mediated GDP/GTP exchange. Data are presented as mean ± SD of three independent experiments performed in technical triplicate. **p < 0.01; ***p < 0.001; ****p < 0.0001; two-tailed paired t-tests.

In addition to the bound PI4P, 5-HT1A is surrounded by an extensive set of lipid molecules at the membrane interface, including 10 cholesterol molecules and three acyl tails from phospholipids (Fig. 2f). These lipid molecules are well defined in the EM density maps (Fig. 2f, Extended Data Fig. 3a-c). Interestingly, one of the two acyl chains of PI4P is sandwiched between two cholesterol molecules, suggesting that cholesterols play a role in stabilizing PI4P binding (Fig. 2f, Extended Data Fig. 4a-c). These observations illustrate that cholesterol directly participates in the binding of PI4P (CHL #2 and #3 in Fig. 2f, Extended Data Fig. 4a-c), therefore providing a structural basis for the long-standing observation that 5-HT1A signaling is modulated by cholesterol and phospholipids[23,24].

Extended Data Fig. 4

Cholesterol regulation in 5-HT1A receptor.

la, The model of the 5-HT1A–Gi complex shows multiple cholesterols bound to the surface of 5-HT1A. The 5-HT1A–Gi complex is shown as surface and lipids are shown as sticks. b, Interactions of CHL #1 with TM1/7 of 5-HT1A. c, Interactions of CHL #2 and #3 with 5-HT1A. d. The effect of cholesterol on the 5-HT potency to activate 5-HT1A. The effects of mutations at the binding residues of cholesterol #1 on 5-HT mediated activation of 5-HT1A (pEC50) were detected by NanoBiT recruitment assays. Data are presented as mean ± SD from at least three independent experiments performed in technical triplicate. *p < 0.05; **p < 0.01; two-tailed paired t-tests.

Ligand activation of 5-HT receptors

Our agonist- and G protein-bound structures are in the active state. Compared to the inactive 5-HT1B structure[20], the active 5-HT1A, 5-HT1D, and 5-HT1e structures undergo a common set of structural rearrangements[21,25], including agonist-induced outward movement at the cytoplasmic end of TM6. As illustrated by the 5-HT1A complexes, the binding of 5-HT or aripiprazole within the orthosteric pocket presses a downward swing of the toggle switch W6x48 (Extended Data Fig. 5c), which leads to conformational changes of F6x44 and I3x40 in the PIF motif (Extended Data Fig. 5d) and R3x50 in the DRY motif. These conformational changes lead to the breakage of the conserved ionic lock between TM3 and TM6, leading to a large outward movement of TM6 by as much as 9 Å at the Cα of E6x30 at the cytoplasmic end of TM6. In addition, the cytoplasmic end of TM7 has an inward movement of ~5.0 Å, which allows Y7x53 from the NPxxY motif, together with Y5x58, to form hydrogen bonds with R3x50 of the DRY motif. These conformational changes open the receptor cytoplasmic pocket for Gi binding.

Extended Data Fig. 5

Basal activity and ligand-induced activation of 5-HT1A.

a, Detection of the ligand-reduced activity and constitutive activity of the human 5-HT1A by NanoBiT G-protein recruitment assay. Three ligands, the full agonist 5-HT, the neutral antagonist WAY-100635, and the inverse agonist methiothepin were used. Data are presented as mean ± SEM of three independent experiments performed in technical triplicate. b, Water molecules are coordinated in the ligand binding pocket of the apo 5-HT1A structure. The density is shown at 3σ cutoff. c-f, Activation of the 5-HT1A by 5-HT and the binding to apo 5-HT1A by water molecules. Toggle switch in 5-HT bound 5-HT1A structure (c) and apo 5-HT1A structure (e). PIF motif of the 5-HT bound 5-HT1A structure (d) and the apo 5-HT1A structure (f). The 5-HT bound 5-HT1A structure is colored with turquoise; the apo 5-HT1A is colored with green; 5-HT is colored with orange. Aligned structures of the inactive state 5-HT1B (inverse agonist methiothepin bound) and the intermediate state 5-HT1B (agonist ERG bound) are colored with gray and light gray. Conformational changes of toggle switch residue W3586.48 and the residue F3546.44, which is part of a conserved PIF motif, are illustrated by arrows. g-i, The hydrogen-bonding network of the ligand-binding pocket observed in the MD simulations. Side view (g) and top view (h) of a hydrogen-bonding network linking the key residues of the active apo 5-HT1A receptor. i, Top view of water molecules accommodated in the inactive apo 5-HT1A receptor. The structure of the apo 5-HT1A is colored in light blue. A representative conformation from the active apo 5-HT1A simulations is colored in light green. A representative conformation from the inactive apo 5-HT1A simulations is colored in gray. The structured water W1 and W2 of the apo 5-HT1A-Gi complex structure are showed as sphere. Putative hydrogen bonds are showed as dash lines.

Basis for basal activation of 5-HT1A

We could generate a stable 5-HT1A-Gi complex in the apo state, consistent with the high basal activity of 5-HT1A, which was inhibited by methiothepin, an inverse agonist[20] and to a less extent, by WAY-100635, a neutral antagonist[26] (Extended Data Fig. 5a). The overall structure of the apo-complex is highly similar to the 5-HT bound complex (Extended Data Fig. 1). The EM map clearly reveals several structured water molecules that form hydrogen bonds with the receptor within the orthosteric binding pocket (Fig. 3a, b, Extended Data Fig. 5b), where the local resolution is in the range of 2.5-2.8 Å (Extended Data Fig. 1h). These water molecules are arranged in the same spatial plane as 5-HT (Fig. 3b). Three water molecules overlap with the positions of the serotonin hydroxyl (W1), indole nitrogen (W2), and primary amine (W3) functionalities, respectively (Fig. 3b). The latter water molecule (W3) forms a hydrogen bond with D1163x32, which is conserved in aminergic GPCRs where it binds the canonical primary amine of endogenous ligands (Fig. 3b). W2 is proximal to the toggle residue W3586x48, which is an important determinant of GPCR activation[27] (Extended Data Fig. 5e). The presence of water molecules resembles all three polar functionalities of serotonin in the active apo-5-HT1A-Gi complex, thus contributing to the stabilization of 5-HT1A in the active conformation.

Fig. 3

The structure of waters in the apo 5-HT1A pocket and comparison between 5-HT and BRL-54443 structures.

a, Waters in the apo 5-HT1A ligand binding pocket with EM maps of waters shown in mesh. Hydrogen bonds are shown with dashed lines. b, Superposition of 5-HT and waters in the structures of the 5-HT1A–Gi complexes. c, Comparison of 5-HT recognition between 5-HT1A and 5-HT1D. d, Comparison of ligand recognition between 5-HT bound 5-HT1D and BRL-54443 bound 5-HT1e. e-g, The difference in the orthosteric binding pockets among 5-HT1A (e), 5-HT1D (f), and 5-HT1e (g). h. Concentration-response curves of 5-HT1A, 5-HT1D, and 5-HT1e activation by 5-HT and BRL-54443. Data are presented as mean ± SEM of three independent experiments performed in triplicate.

To investigate the role of water molecules in the binding pocket, we performed additional MD simulations involving the active apo 5-HT1A receptor and the inactive apo 5-HT1A receptor. Two 100-ns simulations were produced for each system. In the active apo receptor simulations, the functionally important residues D1163x32, S1995x43 and T1213x37 could be stabilized by the hydrogen-bonding network of five water molecules (Extend data Figure. 5g, h). Such hydrogen-bonding network was maintained in 68.8% simulation trajectories of the active apo receptor system. However, in the inactive apo receptor system, the TM6 helix is positioned away from TM5 at the extracellular side, which yields to a larger pocket and accommodates more water molecules (Extend data Figure. 5i). Consequentially, the hydrogen-bonding network linking D1163x32, T1213x37 and S1995x43 seen in the active apo receptor system was only observed in 14.1% simulation trajectories of the inactive apo receptor system.

Basis for 5-HT pan agonism

To determine the basis of pan-agonism of 5-HT, we performed pairwise comparisons of the 5-HT1A, 5-HT1D, and 5-HT1e structures. 5-HT has very similar affinities for 5-HT1A and 5-HT1D (pKi of 8.4 and 8.5, respectively)[28]. A comparison with the binding mode of 5-HT in the 5-HT1A and the 5-HT1D structure reveals a nearly overlapped conformation of 5-HT be adopted in both 5-HT1A and 5-HT1D (Fig. 3c, Extended Data Fig. 6a-f), which were further supported by mutagenesis (Extended Data Fig. 7, Supplementary Fig.2). Residues that make up the 5-HT binding pockets have a moderate-to-high sequence conservation (73% identity and 95% similarity) among 5-HT receptors. Only eight out of 22 ligand pocket residues (within 5 Å around of serotonin) are identical among all 12 serotonin GPCRs (Extended Data Fig.8), whereas the other 14 residues can vary from receptor to receptor. Together, this shows that all parts of 5-HT can form varying interactions within the 12 serotonin GPCR subtypes, providing evidence of how GPCRs can interact with the same endogenous ligand with similar affinities and functional responses, a basis for the pan agonism of serotonin.

Extended Data Fig. 6

Ligand recognition of the 5-HT1A, the 5-HT1D and the 5-HT1e.

a, d, g, j, Conformation of the ligand binding pockets in the serotonin bound 5-HT1A (a), the serotonin bound 5-HT1D (d), the BRL-54443 bound 5-HT1e (g), and the aripiprazole bound 5-HT1A (j). b, e, h, k, Diagram of ligand recognition for serotonin in 5-HT1A (b), serotonin in 5-HT1D (e), BRL-54443 in 5-HT1e (h), and aripiprazole in 5-HT1A (k). c, f, i, l, Ligand binding pockets shown as surfaces. The orthosteric binding pocket is highlighted in orange.

Extended Data Fig. 7

Ligand-binding pocket mutagenesis data by NanoBiT Gi-protein recruitment assay.

Data are shown as mean ± SEM from at least three independent experiments performed in technical triplicate. The EC50 ratio, EC50(mutant)/EC50(WT), represents the shift between the WT and mutant curves, and characterizes the effect of the mutations on receptor activation.

Extended Data Fig. 8

Serotonin binding pocket alignment and ligands affinity among serotonin receptors.

a, Dendrogram and sequence alignment based on residues lining the serotonin binding pocket (5 Å cut-off). Identical residues are marked in white, whereas non-conserved are colored by their physicochemical properties. b, binding affinities (pKi values) for selected ligands of the 5-HT receptors. (https://pdsp.unc.edu/pdspweb/)

Drug recognition of 5-HT receptors

Selective and non-selective ligands of 5-HT receptors have been developed into treatments for many neuronal diseases and the five structures reported in this paper have provided a solid foundation to understand ligand selectivity and drug recognition by serotonin receptors. Specifically, BRL-54443 is a selective agonist for 5-HT1e and 5-HT1F with affinities 10- to 100-fold higher than for other 5-HT receptors and dopamine receptors[29]. G protein recruitment assays showed that BRL-54443 is a full agonist for 5-HT1e with an EC50 of ~7 nM, but has 15-fold and 100-fold weaker potency for 5-HT1D and 5-HT1A, respectively (Fig. 3h, Extended Data Fig.7, Supplementary Fig.2). BRL-54443 in 5-HT1e occupies the same orthosteric pocket as 5-HT in 5-HT1A and 5-HT1D (Fig. 3d). However, E3116x55 in the 5-HT1e pocket forms an extra hydrogen bond with the hydroxyl group of BRL-54443, while the corresponding interaction is missing for the homologous positions S3216x55 in 5-HT1D and A3656x55 in 5-HT1A (Fig. 3d, Extended Data Fig. 6). Mutations of E3116x55 in 5-HT1e to the corresponding residues in 5-HT1D and 5-HT1A reduced BRL-54443 induced 5-HT1e activation by 9-fold and 13-fold, respectively (Extended Data Fig. 7, Supplementary Fig.2). Conversely, 5-HT1A with an A3656x55E mutation enhanced the potency of BRL-54443 (Supplementary Fig.2a). These results indicate that the residue at position 6x55 is a key determinant for BRL-54443 selectivity of 5-HT receptors. Furthermore, mutations on K6x54, which interacts with E6x55 but not directly with BRL-54443, reduced the 5-HT1e mediated Gi recruitment efficacy by BRL-54443 (Supplementary Fig.2c). Mutation of E6x55 and K6x54 of 5-HT1e to the corresponding residues of 5-HT1A (A6x55 and V6x54) abolished the Gi recruitment efficacy (Supplementary Fig.2c). These results suggest the K6x54 stabilize the conformation of residue E6x55 allowing it to form a H-bond with BRL-54443 with a favored distance. Throughout the 5-HT receptor subfamily, only 5-HT1F share the same residues of E6x55 and K6x54 as 5-HT1e, consistent with the selectivity of BRL54443 to 5-HT1e and 5-HT1F over other 5-HT receptors.

In contrast to BRL-54443, 5-carboxamidotryptamine (5-CT), which is a 5-HT mimic with the replacement of the 5-hydroxyl by a bulkier carboxamide, exhibits high affinity (~1 nM) for 5-HT1A, 5-HT1B, and 5-HT1D receptors, which possess smaller Ala and Ser residues in position 6x55, but only modest to weak affinity (~1 μM) for 5-HT1e and 5-HT1F (Extended Data Fig. 8b, 9c). Replacement of E3116x55 for Ser or Ala in 5-HT1e increased the potency of 5-CT by 100 to 400 folds (Extended Data Fig. 9e), supporting that residue 6x55 is the determinant for this observed subtype selectivity[30]. To confirm this hypothesis, we correlated the residue types at 6x55 with the binding affinities of 5-HT, 5-MeOT, 5-CT, and Donitriptan (Extended Data Fig.8b). A regression of the affinity versus residue type as a qualitative variable gave R[2] values of 0.97, 0.79, and 0.71 for Donitriptan, 5-CT, and 5-MeOT, respectively (Extended Data Fig. 9a). A considerably weaker correlation was obtained for 5-HT (R[2] = 0.43), reflecting the promiscuity of 5-HT as a pan agonist of serotonin receptors.

Extended Data Fig. 9

Selectivity of 5-HT1 subfamily receptors.

a, Fitted regression model versus experimental binding affinities of 5-HT, 5-MeOT, 5-CT, and donitriptan for the GPCR 5-HT receptors. b-d, Serotonin- (b), 5-CT- (c), and donitriptan-(d) induced Gi activation assay using NanoBiT for wild type 5-HT1A, 5-HT1D, and 5-HT1e receptors. Data are shown as mean ± SEM from at least three independent experiments performed in technical triplicate. e, 5-CT-induced Gi activation assay using NanoBiT for 5-HT1e, and concentration-response curves for G-protein recruitment signals. Data are shown as mean ± SEM from at least three independent experiments performed in technical triplicate. f, g, The different side chains at TMD (f) and at ECL2 (g) that determine the recognition for donitriptan among 5-HT1A, 5-HT1B, 5-HT1D, and 5-HT1e. h, Docked pose of donitriptan in donitriptan-bound 5-HT1A (right), 5-HT1D (middle), and 5-HT1e (left).

Structural analyses of 5HT1A, 5-HT1D, 5-HT1e also revealed the basis for the recognition of aripiprazole, one of the best-selling antipsychotic drugs. Aripiprazole displays differential affinities to various members of the 5-HT1 subfamily, with a 5 nM affinity for 5-HT1A but with 10-1000 fold weaker affinities for 5-HT1B, 5-HT1D, and 5-HT1e [31]. The extracellular end of TM7, which directly participates in the formation of the extended ligand binding pocket, is shifted outward by approximately 3 Å in 5-HT1A relative to the TM7 position in 5-HT1B, 5-HT1D, and 5-HT1e (Fig. 4). Together with F1123x28 and Y962x63, TM7 stabilizes the quinolinone group of aripiprazole in 5-HT1A (Fig. 4). For 5-HT1B, 5-HT1D, and 5-HT1e, the inward movement of the extracellular end of TM7 and the larger side chain of W3x28 (corresponding to F1223x28 at 5HT1A) would overlap with the same binding space of the quinolinone group, resulting in a lower affinity of aripiprazole to these receptors (Fig. 4). Notably, one of the most visible cholesterols is inserted into a cleft between TM1 and TM7 and serves as a structural chaperone for these two helices (Fig. 2f, Fig. 4). This cholesterol is involved in the shaping of the ligand pocket and stabilizes the positions of TM1 and TM7 near the quinolinone group, resulting in the higher affinity of aripiprazole for 5-HT1A. This is consistent with the important roles of cholesterol in the functional regulation of 5-HT1A [32]. No cholesterol has been observed at the corresponding site in the structures of 5-HT1B [21], 5-HT1D, and 5-HT1e (Fig. 1). In addition, cholesterols are also directly involved in the binding of PI4P to enhance G protein coupling and signaling activity (Extended Data Fig. 4).

Fig. 4

The binding of aripiprazole is regulated by cholesterol.

a, A cholesterol at the surface between TM1 and TM7 of 5-HT1A. b, Top view of aripiprazole bound 5-HT1A structure. The 5-HT1A is shown in surface. c, d, Differences in the conformation of TM7 and its interaction with aripiprazole and cholesterol in 5-HT1A from 5-HT1B (6G79), 5-HT1D and 5-HT1e.

In summary, we have determined five structures of three different serotonin receptors. These structures reveal an unexpected role of phospholipid PI4P and cholesterol in G protein coupling and ligand recognition. The structures also reveal the basis of basal activation and 5-HT pan agonism as well as the basis for drug recognition at serotonin receptors. These observations have a wide implication in mechanistic understanding of serotonin signaling and drug discovery targeting this important family of receptors.

Methods

Constructs

The human 5-HT1A, 5-HT1D, and 5-HT1e were modified to contain the N-terminal thermally stabilized BRIL[33] to enhance receptor expression and the addition of affinity tags including N-terminal Flag tag and His tag. L1253x41W, L1273x41W, and L1113x41W mutations were separately introduced to 5-HT1A, 5-HT1D, and 5-HT1e to improve thermal stability[34,35]. The N-terminal 24-residue truncated 5-HT1A and the full-length receptors of 5-HT1D and 5-HT1e were cloned into the pFastBac (Thermo Fisher) vector using ClonExpress II One Step Cloning Kit (Vazyme Biotech Co.,Ltd). A dominant-negative human Gαi1 was generated by site-directed mutagenesis to incorporate mutations S47N, G203A, A326S, and E245A that improves the dominant-negative effect by weakening a salt bridge that helps to stabilize the interactions with the βγ subunits[19]. Human DNGαi1, human WTGαi1, human Gβ1, human Gγ2 and a single chain antibody scFv16[36] were cloned into pFastBac vectors. These constructs were generated in insect cell expression vectors.

Insect cell expression

Human 5-HT1A, 5-HT1D, and 5-HT1e, human DNGαi1, human Gβ1, and human Gγ2 were co-expressed in Spodoptera frugiperda Sf9 insect cells (Invitrogen) using the baculovirus method (Expression Systems). For the 5-HT1D–Gi and 5-HT1e–Gi complexes, scFv16 was co-expressed to stabilize the protein. As for the apo–5-HT1A–Gi (WT) complex used for GTPase assay, the DNGαi1 was replaced by WTGαi1. Cell cultures were grown in ESF 921 serum-free medium (Expression Systems) to a density of 2-3 million cells per ml and then infected with four separate baculoviruses at a suitable ratio. The culture was collected by centrifugation 48 h after infection and cell pellets were stored at -80°C.

Complex purification

Cell pellets were thawed in 20 mM HEPES pH 7.4, 50 mM NaCl, 10 mM MgCl2 supplemented with Protease Inhibitor Cocktail (Bimake). For the apo–5-HT1A–Gi complex, 25 mU/ml apyrase (Sigma) was added; for the 5-HT- or aripiprazole-bound 5-HT1A–Gi complexes, 20 μM 5-HT (TargetMol) or aripiprazole (TargetMol) and 25 mU/ml apyrase (Sigma) were added; for the 5-HT1D–Gi–scFv16 complex, 20 μM 5-HT (TargetMol) and 25 mU/ml apyrase (Sigma) were added; for the 5-HT1e–Gi–scFv16 complex, 20 μM BRL-54443 (TargetMol) and 25 mU/ml apyrase (Sigma) were added. The suspension was incubated for 1 h at room temperature and the complex was solubilized from the membrane using 0.5% (w/v) n-dodecyl-β-d-maltoside (DDM, Anatrace) and 0.1% (w/v) cholesteryl hemisuccinate (CHS, Anatrace) for 2 h at 4°C. Insoluble material was removed by centrifugation at 65,000 g for 30 min and the solubilized complex was immobilized by batch binding to Talon affinity resin. After that, the resin was packed and washed with 20 column volumes of 20 mM HEPES pH 7.4, 100 mM NaCl, 5 mM MgCl2, 0.01% (w/v) lauryl maltose neopentylglycol (LMNG, Anatrace), and 0.002% (w/v) CHS. The complex was then eluted in buffer containing 300 mM imidazole and concentrated using an Amicon Ultra Centrifugal Filter (MWCO 100 kDa). The complex was subjected to size-exclusion chromatography on a Superdex 200 Increase 10/300 column (GE Healthcare) pre-equilibrated with Size Buffer containing 20 mM HEPES pH 7.4, 100 mM NaCl, 0.00075% (w/v) LMNG, 0.00025% (w/v) GDN (Anatrace), and 0.00015% CHS to separate complex from contaminants. For the ligand-bound complexes, 20 μM corresponding ligand was contained in the Size Buffer. Eluted fractions were evaluated by SDS-PAGE (Supplementary Fig.1) and consisting of receptor–Gi-protein complex were pooled and concentrated for EM experiments.

Peptidisc reconstitution and GTPase assay

Peptide NSPr (FAEKFKEAVKDYFAKFWDPAAEKLKEAVKDYFAKLWD) was used for peptidisc reconstitutions with different lipid compositions. The reconstitution was based on methods previously reports[22,37]. In brief, the purified 5-HT1A–Gi (WT) complex was mixed with NSPr, lipids, LMNG, and cholesterol at optimized ratios. Reconstitution mixture was incubated overnight at 4℃. Then the samples were subjected to Bio-Beads SM2 (50 mg/ml; Bio-Rad) twice for total 8 h at 4℃ to remove detergent. The mixture was spun down, and the supernatant was loaded on to TALON beads to remove empty peptidiscs. After binding for 3 h at 4℃, the beads were washed with 20 mM HEPES pH 7.4, 100 mM NaCl and the lipid-containing 5-HT1A–Gi peptidiscs were eluted with 20 mM HEPES pH 7.4, 100 mM NaCl, 250 mM Imidazole. The samples were applied onto a Superdex 200 Increase 10/300 column (GE Healthcare) pre-equilibrated with 20 mM HEPES pH 7.4, 100 mM NaCl. The peak corresponding to the 5-HT1A–Gi peptidiscs was collected for the GTPase-Glo assay. A modified protocol of the GTPase-Glo™ Assay (Promega) was used to detect the GTPase activity of the G-protein[10]. In this assay, the amount of GTP depends on the degree of hydrolysis by G-protein or activated receptor-bound G-protein. The remaining GTP be converted to ATP by GTPase-Glo reagent, then converted to luminescent signal by detection reagent. To analyses the regulating effect of lipids, the GTPase reaction with 5-HT1A–Gi peptidisc started in an assay buffer containing 20 mM HEPES, pH 7.4, 100 mM NaCl, 1 mM DTT, 5 μM GDP and 5 μM GTP. After 60 min incubation, the detection reagent was added and then incubated for 30 min at room temperature. Luminescence was measured using EnVison® (PerkinElmer).

NanoBiT G-protein recruitment assay

Analysis of Gi-protein recruitment was performed by using a modified protocol of the NanoBiT system (Promega) described previously[30]. NanoBiT system is a two-subunit system based on NanoLuc luciferase that can be used for functional detection for binding and dissociation of GPCR and G-protein. LgBiT (17.6 kDa) of the NanoBiT luciferase was fused to the C-terminal of GPCR with 15-amino acid flexible linkers. SmBiT was C-terminally fused to Gβ subunit with a 15-amino acid flexible linker. Human 5-HT1A-LgBiT, 5-HT1D-LgBiT, and 5-HT1e-LgBiT, WTGαi1, SmBiT-fused Gβ1, and Gγ2 were co-expressed in Spodoptera frugiperda Sf9 insect cells using the baculovirus method (Expression Systems). Cell cultures were grown to a density of 2-3 million cells per ml and then infected with four separate baculoviruses at a suitable ratio. The culture was collected by centrifugation 48 h after infection and cell pellets were collected with PBS. The cell suspension was dispensed in a white 384-well plate at a volume of 40 μl per well and loaded with 5 μl of 90 μM coelenterazine (Yeasen) diluted in the assay buffer. Test compounds (5 μl) were added and incubated for 3-5 min at room temperature before measurement. Luminescence counts were normalized to the initial count and fold-change signals over vehicle treatment were used to show G-protein binding response.

Cryo-EM grid preparation and data acquisition

Three microliters of the purified agonist-bound 5-HT1A, 5-HT1D, and 5-HT1e complexes at ~5 mg/ml, 30 mg/ml, 45 mg/ml, respectively, were applied onto a glow-discharged Quantifoil R1.2/1.3 200-mesh gold holey carbon grid. The grids were blotted for 3 s under 100% humidity at 4 °C and then vitrified by plunging into liquid ethane using a Vitrobot Mark IV (Thermo Fischer Scientific). For the 5-HT1A complexes, cryo-EM data collection was performed on a Titan Krios at 300 kV accelerating voltage in the Center of Cryo-Electron Microscopy, Zhejiang University (Hangzhou, China) and the micrographs were recorded using a K2 Summit direct electron detector (Gatan) in counting mode at a calibrated magnification of 1.014 Å per pixel. Micrographs were obtained at a dose rate of about 8.0 e/Å[2]/s with a defocus ranging from -1.0 to -3.0 μm. Each micrograph was dose-fractionated to 40 frames with a total exposure time of 8 s. For the 5-HT1D and 5-HT1e complexes, cryo-EM data collection was performed on a Titan Krios at 300 kV accelerating voltage in the Center of Cryo-Electron Microscopy, Shanghai Institute of Materia Medica, Chinese Academy of Sciences (Shanghai, China). The micrographs were recorded using a K3 Summit direct electron detector (Gatan) with a Gatan energy filter (operated with a slit width of 20 eV) (GIF). The microscope was operated at a magnification of 47,847 in counting mode, corresponding to pixel size of micrograph at 1.045 Å. The images were recorded at a dose rate of about 26.7 e/Å[2]/s with a defocus ranging from -1.0 to -3.0 μm. The total exposure time was 3 s and intermediate frames were recorded in 0.083 s intervals, resulting in a total of 36 frames per micrograph. A total of 4179, 5269, 2005, 4375, and 5249 movies were collected for the apo-, 5-HT bound-, aripiprazole bound- 5-HT1A–Gi, 5-HT bound- 5-HT1D–Gi, and BRL-54443 bound- 5-HT1e–Gi complexes, respectively.

Image processing and 3D reconstruction

Image stacks were subjected to beam-induced motion correction using MotionCor2.1[38]. Contrast transfer function (CTF) parameters for each non-dose weighted micrograph were determined by Gctf[39]. Automated particle selection and data processing were performed using Relion-3.0-beta2[40]. For the dataset of 5-HT bound 5-HT1A–Gi complex, automated particle selection yielded 5,279,538 particles. The particles were extracted on a binned dataset with a pixel size of 2.028 Å and were subjected to reference-free 2D classification, producing 898,450 particles with well-defined averages. The map of NTSR1–Gi1 complex (EMDB-20180)[41] low-pass filtered to 40 Å was used as an initial reference model for 3D classification, which produces two good subsets showing clear structural features accounting for 559,323 particles. Further 3D classifications focusing the alignment on the complex except AHD of the G⍺, produced three high-quality subsets accounting for 472,338 particles. These particles were subsequently subjected to Bayesian polishing, CTF refinement and 3D refinement, which generated a map with an indicated global resolution of 3.0 Å at a Fourier shell correlation of 0.143. For the dataset of apo–5-HT1A–Gi complex, particle selection yielded 2,719,825 particles, which were subjected to reference-free 2D classification, producing 1,922,844 particles for further processing. The map of 5-HT–5-HT1A–Gi complex low-pass filtered to 20 Å was used as an initial reference model for 3D classification, which produces one good subset accounting for 680,449 particles. Further 2 rounds of 3D classifications focusing the alignment on the complex except AHD of the G⍺, produced the high-quality subsets accounting for 245,886 particles. These particles were subsequently subjected to Bayesian polishing, CTF refinement and 3D refinement, which generated a map with an indicated global resolution of 3.0 Å at a Fourier shell correlation of 0.143. For the dataset of aripiprazole–5-HT1A–Gi complex, a total of 1,486,169 particles were automatically picked and extracted on a binned 2 dataset. These particles were subjected to reference-free 2D classification, producing 748,746 particles with well-defined averages for further processing. The map of apo–5-HT1A–Gi complex low-pass filtered to 20 Å was used as an initial reference model for 3D classification, which produces one good subset. Further 3D classifications focusing the alignment on the complex, produced one high-quality subset accounting for 154,241 particles. These particles were subsequently subjected to 3D refinement and Bayesian polishing, which generated a map with an indicated global resolution of 3.1 Å according to the 0.143 criteria of the FSC. For the dataset of 5-HT–5-HT1D–Gi complex, particle selection yielded 4,212,737 particles. The particles were subjected to reference-free 2D classification, producing 2,361,598 particles with good features. The map of 5-HT–5-HT1A–Gi complex low-pass filtered to 60 Å was used as an initial model for 3D classification, which produces one good subset accounting for 896,076 particles. Further 3D classifications focusing the alignment on the receptor and the receptor–Gi complex, produced one high-quality subset accounting for 141,501 particles. These particles were subsequently subjected to CTF refinement and Bayesian polishing, which generated a map with an indicated global resolution of 2.9 Å. For the dataset of BRL54443–5-HT1e–Gi complex, automated particle selection yielded 4,977,538 particles, which were subjected to 2D classification. 2,175,234 particles with well-defined averages were selected for further 3D classification processing. The map of 5-HT–5-HT1A–Gi complex low-pass filtered to 60 Å was used as an initial reference model for 3D classification, which produces one good subset accounting for 1,156,885 particles. Further 3D classifications focusing the alignment on the receptor and the entire complex, produced the two subsets with good receptor density. The two subsets accounting for 163,354 particles were subsequently subjected to CTF refinement, Bayesian polishing and 3D refinement, which produced a map with an indicated global resolution of 2.9 Å according to 0.143 criterion of the FSC. Local resolution was determined using the Bsoft package[42] with half maps as input maps.

Structure determination and refinement

The cryo-EM structures of the apo–5-HT1A–Gi, 5-HT–5-HT1A–Gi, aripiprazole–5-HT1A–Gi, 5-HT–5-HT1D–Gi–scFv16, and BRL-54443–5-HT1e–Gi–scFv16 complexes were solved using 5-HT1B (PDB code: 6G79)[21] and rhodopsin–Gi complex (PDB code: 6CMO)[43] as initial models. The models were docked into the electron microscopy density maps using Chimera[44], followed by iterative manual adjustment and rebuilding in COOT[45] and ISOLDE[46], against the cryo-EM electron density maps. Real space and reciprocal refinements were performed using Phenix[47]. The model statistics was validated using MolProbity[48]. Structural figures were prepared in Chimera[44], ChimeraX[49] and PyMOL (https://pymol.org/2/). The final refinement statistics are provided in Extended Data Table 1.

Extended Data Table 1

Cryo-EM data collection, refinement and validation statistics

	Apo–5-HT1A–Gi(EMD-30971)(PDB 7E2X)	5-HT–5-HT1A–Gi(EMD-30972)(PDB 7E2Y)	Aripiprazole–5-HT1A–Gi(EMD-30973)(PDB 7E2Z)	5-HT–5-HT1D–Gi–scFv16(EMD-30974)(PDB 7E32)	BRL-54443–5-HT1e–Gi–scFv16(EMD-30975)(PDB 7E33)
Data Collection and Processing
Magnification	49,310	49,310	49,310	47,847	47,847
Voltage (kV)	300	300	300	300	300
Electron exposure (e-/Å2)	64	64	64	70	70
Defocus range (μm)	-1.0 ~ -3.0	-1.0 ~ -3.0	-1.0 ~ -3.0	-1.0 ~ -3.0	-1.0 ~ -3.0
Pixel size (Å)	1.014	1.014	1.014	1.045	1.045
Symmetry imposed	C1	C1	C1	C1	C1
Initial particle projections (no.)	2,719,825	5,279,538	1,486,169	4,212,737	4,977,538
Final particle projections (no.)	245,886	472,338	154,241	141,501	163,354
Map resolution (Å)	3.0	3.0	3.1	2.9	2.9
FSC threshold	0.143	0.143	0.143	0.143	0.143
Map resolution range (Å)	2.2-5	2.2-5	2.2-5	2.2-5	2.2-5
Refinement
Initial model used (PDB code)	6G79, 6CMO	6G79, 6CMO	6G79, 6CMO	6G79, 6CMO	6G79, 6CMO
Model resolution (Å)	3.0	3.2	3.2	3.0	3.1
FSC threshold	0.5	0.5	0.5	0.5	0.5
Model resolution range (Å)	50-3.0	48-3.2	48-3.2	48-3.0	48-3.1
Map sharpening B factor (Å2)	-115.387	-108.16	-125.267	-88.65	-88.26
Model Composition   Non-hydrogen atoms	7287	7035	7073	8839	8610
Protein residues	890	873	880	1116	1091
Lipids	15	5	6	1	1
B factors (Å2)   Protein	160.4	100.4	55.5	118.5	115.6
Ligand	161.6	109.7	60.9	55.5	99.8
R.m.s. deviations   Bond lengths (Å)	0.005	0.005	0.005	0.017	0.006
Bond angles (°)	1.13	1.16	1.11	1.10	0.96
Validation
MolProbity score	1.42	1.51	1.35	1.72	1.59
Clashscore	4.48	4.19	4.39	9.16	10.02
Poor rotamer (%)	0.00	0.00	0.00	0.00	0.11
Ramachandran Plot   Favored (%)	96.9	95.6	97.3	96.4	97.7
Allowed (%)	3.1	4.4	2.7	3.6	2.3
Disallowed (%)	0.0	0.0	0.0	0.0	0.0

Molecular docking

The structures of 5-HT1A, 5-HT1B (PDB code: 6G79), 5-HT1D, and 5-HT1e were employed to build the molecular docking models. All ligands and water molecules in these structures were removed to generate the models. The structures of the small-molecule ligands of interest were download from the PubChem database (www. pubchem.ncbi.nlm.nih.gov). A ligand was docked to the proposed ligand-binding pocket of its receptor using Schrodinger Glide software in SP mode with default parameters[50]. The ligand was initially placed in the center of the pocket and was constrained to move within a 1 nm diameter sphere, where it was allowed to move freely during the docking process. The extended conformation searches were performed using the Lamarckian Genetic Algorithm. Docked poses with most negative docking scores were selected for analysis.

Regression models for ligand affinity

Partial least squares regression models were made for the binding affinity of 5-HT, 5-CT, 5MeOT, and Donitriptan to the 5-HT receptors (Extended Data Figure 8a) and the amino acid residue at position 6x55 as a qualitative variable. Separate models were made for each ligand and the applied number of principal components was determined by cross-validation. The calculations were made using Simpca-P version 11.0 (https://umetrics.com).

Molecular dynamics simulation.

The active apo 5-HT1A simulation model was built based on the apo–5-HT1A–Gi complex structure. All structured water molecules were included in the simulation system. Using homology modeling method, the inactive 5-HT1A simulation model was built based on the structure of inactive 5-HT1B (PDB code: 5V54) with Modeller[51]. The default parameters were employed to construct the models. The missing backbone and sidechains were added. The models with the lowest root mean square deviations from their template structures were selected. To build a simulation system, we placed the complex model into a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine lipid bilayer. The lipid embedded complex model was solvated in a periodic boundary condition box (160 Å × 160 Å × 160 Å) filed with TI3P water molecules and 0.15 M KCl using CHARMM-GUI[52]. Each system was replicated to perform two independent simulations. On the basis of the CHARMM all-atom force field[53,54], molecular dynamics simulations were conducted using GROMACS 5.1.4[55,56]. After 100 ns equilibration, a 100 ns production run was carried out for each simulation. All productions were carried out in the NPT ensemble at temperature of 303.15 K and a pressure of 1 atm. Temperature and pressure were controlled using the velocity-rescale thermostat[57] and the Parrinello-Rahman barostat with isotropic coupling[58], respectively. Equations of motion were integrated with a 2 fs time step, the LINCS algorithm was used to constrain bond length[59]. Nonbonded pair lists were generated every 10 steps using distance cutoff of 1.4 nm. A cutoff of 1.2 nm was used for Lennard-Jones (excluding scales 1-4) interactions, which were smoothly switched off between 1 and 1.2 nm. Electrostatic interactions were computed using particle-mesh-Ewald algorithm with a real-space cutoff of 1.2 nm[59].