Exploring the Conformational Impact of Glycine Receptor TM1-2 Mutations Through Coarse-Grained Analysis and Atomistic Simulations

Introduction

Pentameric ligand-gated ion channels (PLGICs) form a large family of integral membrane proteins with a central role in signal transduction from prokaryotes to eukaryotes (Dacosta and Baenziger, 2013; Changeux, 2014; Taly et al., 2014). Their ring-like pentamer architecture, with a fivefold symmetry axis centered on the ion-conducting pore, is conserved from bacteria to humans and mediates an incredibly sophisticated mechanism to allosterically propagate signals from the extracellular binding site to an ionic gate situated up to 50 Å away. In animals, PLGICs share a conserved extracellular cysteine bridge, which gives its name to the so-called Cys-loop family of ionotropic receptors. Given their key role in chemical synapses, Cys-loop receptors are major drug targets in neurological conditions from Alzheimer to rare genetic diseases. Mostly expressed at post-synaptic neurons, upon pre-synaptic neurotransmitter release, they mediate passive ion fluxes that shift the membrane potential. Depending on pore-lining residues, PLGICs are selective for cations like sodium ions (Na⁺) and result in excitatory effects, e.g. the nicotinic acetylcholine receptor (nAch-R) or the serotonin type-3 receptor (5-HT₃-R), or are selective for anions like chloride ions (Cl⁻) resulting in inhibitory effects, e.g., the α-aminobutyric acid receptor (GABA_A-Rs) or the glycine receptor (Gly-Rs) (Figure 1).

FIGURE 1

FIGURE 1. Structure and mutations of the glycine receptor (GlyR). (A,B) Human GlyR closed (5cfb) and desensitized (5vdh) structures used as reference in this study (left), zoom onto TM2 pore helices with key pore residues highlighted with licorice representation (center) and corresponding profile of pore radii calculated by HOLE (right). (C) COSMIC GLRA3 mutations colored by frequency; note the top mutation in blue (S241L) at the end of pore lining TM2 (orange). (D) Chord plot of the network of conserved contacts from 5CFB, with doubly HPX-COSMIC mutated residues in bold (top) and mutational clustering across GLRA1-3 genes weighted by conservation score (bottom). Note how multiple peaks overlap with positions mutated in startle disease (red) and the highest one corresponds to the TM1-2 region where HPX changes concentrate. See related Supplementary Figure S1, Supplementary Table S1.

Whether prokaryotic or eukaryotic, PLGICs share a topology characterized by a large N-terminal ligand-binding extracellular domain (ECD), followed by four transmembrane domains (TM1–4). Monomers assemble into a pentameric cylinder, with the five orthosteric ligand binding sites (LBS) located at ECD subunit interfaces, TM2 helices shaping the ion channel pore across the central symmetry axis, and TM4 helices facing the plasma membrane (Figures 1A,B, left). This universal topology is linked to incredibly conserved functional mechanisms to regulate ion gating (Gielen and Corringer, 2018). Classically, PLGIC activation was interpreted with simple two-state models like the MWC (Monod et al., 1965), in which receptors spontaneously sample resting/shut and active/open states until agonist binding shifts the equilibrium. Nevertheless, single-channel studies and growing structural data show a far more complex conformational cycle, with striking similarities across species. In the initial pre-activation or “priming” step, agonist binding stabilizes the ECD in a contracted higher-affinity conformation (“un-blooming”), triggering a key revolving motion of TM2-3 loop at the ECD–TMD interface and subsequent ECD–TMD rotation in opposite directions (“quaternary twist”) (Nemecz et al., 2016). This sequence of events is captured by the so-called “locally closed” conformations, where the ECD has undergone the transition toward the active state-like conformation, but the TMD still remains in a resting conformation (Prevost et al., 2012). In the activation step, global twisting couples to cooperative tilting of pore-lining helices (“iris-like gating”) (Martin et al., 2017), which widens the upper part of the channel (“activation gate”, TM2 9′ and 13′; see Figures 1A,B right), formed by two or three rings of hydrophobic residues that create a barrier to ion permeation (Althoff et al., 2014; Du et al., 2015). Once the gate is open, ions flow according to their electrochemical gradient and channel selectivity, determined by the “selectivity filter” at the cytoplasmic end of the pore (TM2 -1′ or -2′). Remarkably, this “unbloom-and-twist” allosteric mechanism that propagates a signal (ligand) from the ECD to a remote TMD pore gate is encoded in the PLGIC fold, as predicted by elastic network models (ENMs) (Bahar et al., 2010), and it has been further confirmed by molecular dynamics (MD) (Calimet et al., 2013) and coarse-grained eBDIMS simulations (Orellana et al., 2016). Apart from agonist-elicited activation, PLGICs also share another key physiological property: desensitization, in which the sustained agonist presence causes the channels to transit from the active to an agonist-bound inactive state (Katz and Thesleff, 1957) to prevent over-activation. Once the agonist disappears, receptors slowly recover, although the detailed transitions from agonist-unbound desensitized states to unbound resting or open states are still unclear. Desensitization mainly involves pore closure at its intracellular end (TM1-2 loop), which therefore also acts as the main “desensitization gate” (Hibbs and Gouaux, 2011).

Among Cys-loop receptors, strychnine-sensitive GlyR, the major inhibitory ionotropic receptor in the brainstem and spinal cord, has become by far the better characterized (Howard, 2021): there are over 40 structures deposited in the Protein Data Bank, mostly in open, closed and desensitized-like states (Du et al., 2015; Huang et al., 2015, 2017b, 2017a; Kumar et al., 2020; Yu J. et al., 2021, Yu et al., 2021 H.); a few are trapped in a super-open states, whose physiological significance has been questioned, specially by MD (Cerdan et al., 2018; Cerdan and Cecchini, 2020; Dämgen and Biggin, 2020; Dämgen et al., 2020). Apart from these rich structural data, GlyR stands out among the eukaryotic PLGICs due to its role in neurological diseases and particularly in a rare Mendelian condition known as hyperekplexia (HPX) or “startle disease” (Lynch et al., 2017). Similar to the way the alkaloid strychnine antagonizes glycine binding, hyperekplexia disrupts glycinergic neurotransmission, resulting in exaggerated “startle” responses and muscle stiffness. Analysis of hyperekplexia patients has resulted in an exceptional amount of information on GlyR mutations and their phenotypic impact being gathered during the past two decades (Lewis et al., 1998; Chung et al., 2010; Bode et al., 2013; Bode and Lynch, 2014). As the GABA receptor, GlyR usually functions as an heteropentamer of alpha and beta subunits, but only alpha subunits form functional homopentamers. Over 50 mutations have been linked to hyperekplexia, mostly targeting the alpha1 subunit (GLRA1), thus providing an exhaustive mutational scanning map in terms of structural areas where missense changes result in GlyR disruption (Chung et al., 2010; Bode et al., 2013; Bode and Lynch, 2013, 2014). In general, HPX mutations (Supplementary Table S1, Figure 1D) are either recessive and typically associated with low surface expression (loss-of-function), or dominant, mostly located around TM2 and causing prolonged desensitization and/or spontaneous activation (gain-of-function), which leads to reduced maximal currents.

Despite the fact that its central function is synaptic signaling, PLGICS are also expressed in non-neural cells where they play a diversity of roles, including stem cell and cancer proliferation (Young and Bordey, 2009; Zhang et al., 2013; Bhattacharya et al., 2021). GlyR is no exception and is known to be expressed in cells as diverse as hepatocytes, spermatozoa, pancreatic, endothelial, or renal cells (Lynch, 2004; Van den Eynden et al., 2009). As tumor genomic screenings advance, an increasing number of GlyR mutations are being reported in a surprising variety of tumors. For the three GlyR genes, GLRA1-3, cancer-reported mutations display an intriguing clustering partially overlapping with HPX positions (Figures 1C,D, Supplementary Table S1). Here, we perform a preliminary exploration of selected GlyR mutations using as model the GLRA3 homomer, which apart of being structurally well characterized in humans also carries the highest frequency of mutations, specially focused on the lower TM2 section (Figure 1C). Based on mutational clustering, conservation, and ENM analysis, we select for MD study three TM1 and 2 mutations found in tumors, which are close or overlap with hyperekplexia-mutated positions (S241L, R252S, and V260M) of uncertain functional effect. Our results suggest that these residue substitutions are far from neutral but profoundly impair channel permeability, selectivity, and desensitization.

Methods

Structural Data, Sequence Alignment, and Conserved Network Analysis

The reference wildtype structures for closed (PDB ID: 5CFB, Figure 1A) and desensitized state (PDB ID: 5VDH; Figure 1B) human GlyR GLRA3 homopentamer were obtained from the Protein Data Bank and the series of mutants, S241L, R252S, and V260M were constructed from both states, resulting in eight systems. To select candidate mutations for further study, we focused on the resting 5CFB structure to perform a simple weighting analysis based on the 3D space distribution and degree of conservation as in Orellana et al. (2019b). First, we fetched the information for GLRA1-3 missense mutations reported in the COSMIC database (Tate et al., 2019) (Supplementary Figures S1, S2). Conservation scores were retrieved from the ConSurf database (Ashkenazy et al., 2016); to get a simpler overview of residue conservation, sequences for GLRA1-3, GBRA1, and the intensely studied PLGICS G.violaceus GLIC and C. elegans GluCL were retrieved from the UniProt database (The UniProt Consortium, 2021) and aligned with ClustalW (Supplementary Figure S1). Then, in order to evaluate spatial 3D clustering, we applied a simple counting algorithm: each amino acid is represented by its C-alpha carbon, and the number of mutations reported for each position and its neighbors within a 3D-sphere of cutoff radius 9 Å (typical to evaluate residue pairwise interactions) is added to obtain a raw number of hits (i.e., reported mutations within the 3D-sphere) (Supplementary Figure S2A), which are then weighted according to their ConSurf scores (Figure 1D, bottom). Hence, random isolated mutations or changes in non-conserved areas are filtered out to obtain a final estimate of the conserved spatial mutation concentration around each amino acid. Finally, the contact network between highly conserved residues (i.e., those with ConSurf scores 8–9) was plotted with a chord diagram to visualize interactions across different regions and mutated areas (Figure 1D, top).

Principal Component Analysis

Principal component analysis (PCA) (Jolliffe, 2002) is a statistical technique to reveal dominant patterns in noisy data. The diagonalization of the covariance matrix of the system allows obtaining the major axis for statistical variance or principal components (PCs). In this way, complex multidimensional data are mapped to a reduced set of coordinates, which contain the dominant trends explaining data variation. PCA has been widely applied in structural biology to analyze ensembles, usually coming from MD simulations (Amadei et al., 1993, 1996). Protein structures are aligned to a reference in order to compute a covariance matrix, which describes the mean-square deviations in atomic coordinates from their mean position (diagonal elements) and the correlations between their pairwise fluctuations (off-diagonal elements). Diagonalization then yields a set of eigenvectors (principal components, PCs) and eigenvalues representing the motions that explain the variation in the atomic coordinates. In the structurally rich ensembles here analyzed, the first two PCs contain, on average, around 60–80% of the ensemble structural variation (Orellana et al., 2016), and provide excellent coordinates to assess mutation effects on MD sampling (Chen et al., 2021a; 2021b). On this framework, a structure i containing N residues is thus accurately characterized by its projections onto the conformational space defined by the major components, PC_k (k = 1, 2 .… 3N-6) (see Figure 2)

$P{C}_k= \left[T_{i-0}\right]\cos \alpha .(1)$

where T _i-0 is the difference between the coordinates of i-structure and the apo reference, PC is one of the major axes, and α is the angle formed by PC_k and T _i-0. Here, we retrieve all available GlyR structures in the Protein Data Bank, corresponding to H. sapiens, D. rerio, and S. domesticus, and after elimination of structures with missing gaps and alignment to the common conserved core, we obtain an ensemble of 33 structures (3jad, 3jae, 3jaf, 5cfb, 5tin, 5vdh, 5vdi, 6plo. 6plp, 6plq, 6plr, 6pls, 6plt, 6plu, 6plv, 6plw, 6plx, 6ply, 6plz, 6pm0, 6pm2-6, 6ubs, 6ubt, 6ud3, 6vm0, 6vm2-3, 7mlu) that aligned to 5cfb with low RMSD. As previously shown by us (Orellana et al., 2016), major PCs are captured, i.e., correlate with the heuristic variables that typically characterize PLGIC conformations (Supplementary Figure S2B). Projections were used to also track the time evolution of trajectories (Figure 4).

FIGURE 2

FIGURE 2. Principal component analysis and elastic network modeling of GlyR. (A) Simplified scheme of the conformational cycle depicting the three main meta-stable states: closed, open, and desensitized, with ECD in blue and TMD in pink depicting blooming and open and closed pore TM2 orientations (see the main text). (B) Experimental GlyR ensemble containing human and zebrafish structures (top) and principal components from the GlyR ensemble (bottom): the major component tracks the anticorrelated ECD-TMD twist, while the second captures iris-like pore gating. (C) Projections of the GlyR ensemble onto PC1-2, with structures colored according to their radius. Note how structures separate onto three clusters corresponding to the assigned functional states, with the only exception of open blocked 6ud3. The human closed structure 5cfb is used as a reference and thus projects at the origin. See related Supplementary Figure S1, Supplementary Table S1.

Molecular Dynamics Simulations and Essential Dynamics

We used the CHARMM-GUI web interface to build systems for MD simulation (Jo et al., 2008, 2017), which allowed us to repair missing residues in the crystal structures and build a membrane bilayer containing ∼300 phospholipid POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) molecules. A hexagonal water box was used to reduce the cell volume and overall system size. Membrane-embedded proteins were then solvated using the TIP3P water model. The CHARMM36m force field was used to describe the system (Huang et al., 2016). Potassium (K⁺) and chloride (Cl⁻) ions were added to maintain the physiological salt (150 mM KCl) concentration to mimic intracellular conditions. Energy minimization, equilibration, and production runs were carried out with GROMACS (Pronk et al., 2013; Abraham et al., 2015), following the CHARMM-GUI Membrane Builder standard protocols (Wu et al., 2014). The temperature was maintained at 303.0 K using the Nose–Hoover thermostat (Nosé, 1984; Hoover, 1985) and pressure was set to 1.0 bar using the Parrinello–Rahman barostat (Parrinello and Rahman, 1981) with semi-isotropic pressure coupling. Hydrogen bonds were constrained using the LINCS algorithm (Hess et al., 1997), short-range van der Waals (vdW), and electrostatic interactions cutoffs were set to 12 Å, and long-range electrostatic interactions were described using the particle mesh Ewald approach (Ewald, 1921; Essmann et al., 1995) with periodic boundary conditions. Production runs were carried out using a 2 fs time step and writing at every 1 ps interval. Each system was simulated for 300 ns, with four replicas starting from different random seeds. Therefore, in total, we simulated 1.2 μs for each one of the eight GlyR systems, i.e., from 5CFB/5VDH, with WT/S241L/R252S/V260M sequences. To filter out noise and extract the main collective motions, we performed essential dynamics (ED) (Amadei et al., 1993; Daidone and Amadei, 2012) for each system’s 1.2 μs meta-trajectory (Figure 3A) using in-house scripts. The relative free energy landscape (FEL) at 300 K was obtained from the probability distribution of the reaction coordinate, R (PC1, PC2) (Figure 3B). GROMACS tools with defaults were used to perform the RMSD cluster analysis of the TMD and calculate the average TM1-2 helicity. For cluster analysis (Figure 8), all simulated systems were combined at a 1 ns interval to yield a single and long C_α atom 9.6 μs trajectory of 9600 frames, and an RMSD cut-off of 0.20 nm was selected to obtain lesser and larger clusters. See the summary of all simulations in Supplementary Table S2.

FIGURE 3

FIGURE 3. Conformational space sampling by elastic network models and essential dynamics across the three major axis for transition defined by closed (5CFB), desensitized (5VDH), and super-open (6PM0) structures from each PC cluster. (A) Overlap between the first 10 ENM modes computed from the three state representatives versus the three transitions (left) and between the first 10 ED modes from 5CFB and 5VDH simulations versus the three transitions (right). (B) Relative free energy landscapes for 1.2 us meta-trajectories upon projection onto PC1-2. Note how simulations from 5CFB (top row) extend in the direction of super-open structures, and while for WT GlyR populate two minima, for mutations sample an elongated one. On the contrary, simulations from 5VDH evolve toward the closed cluster, while for WT GlyR they point to 5cfb, mutants tend to evolve to collapsed closed states that project at negative PC2 values. All simulations from 5VDH are indicated with *.

Elastic Network Normal Mode Analysis

Elastic network models (ENM) (Atilgan, 2018) are minimal coarse-grained representations of protein structures as beads-and-springs. In the anisotropic network model (ANM), the potential energy of a protein structure—assumed to be at an energy minimum—is defined by a network of interactions with nodes at the C^α atoms coupled by uniform springs (Tirion, 1996; Atilgan et al., 2001). Once the pair force constants K_ij are defined, the network potential energy is approached by the parabola:

$E= {\displaystyle \sum _{i=j}{K}_{\text{ij}}{\left(R_{\text{ij}}-R_{ij}^0\right)}^2}(2)$

where R_ij and R_ij⁰ are the instantaneous and equilibrium distances between all nodes. Within normal mode analysis (Case, 1999), the Hessian matrix of the energy second derivatives is diagonalized, yielding a set of orthonormal 3N-6 eigenvectors representing the normal modes, which have been shown to accurately predict intrinsic collective motions. Here, we use the nearest-neighbors MD-derived ED–ENM algorithm, which predicts experimental conformational changes (Orellana et al., 2010) and, implemented onto ED–ENM Brownian dynamics (eBDIMS), entire sequences of intermediate states along transition pathways (Orellana et al., 2016; 2019a). Average thermal fluctuations for each residue pair i, j and each residue i are evaluated as in Atilgan et al. (2001) (Supplementary Figure S2C).

Essential Dynamics and Normal Mode Analysis Similarity Against Experimental Transition Vectors

An experimental conformational change between two conformations R₁ and R₂ can be described by the normalized difference 3N-dimensional vector ∆r = (R₂ - R₁)/║R₂ - R₁║ between the two sets of alpha-carbon coordinates, after optimal superimposition of the structures. Therefore, given a motion space from ED or NMA, the degree of similarity or overlap between the directions of the experimental ∆r vector and a given kth mode υ_k is measured by their angle cosinus (Yang et al., 2009):

${\alpha }_k= \frac{\mathrm{\Delta }r.v_k}{\Vert \mathrm{\Delta }r\Vert \Vert v_k\Vert }(3)$

where ∙ designates the dot product and the bars denote the vectors’ modulus or magnitude; a cosinus close to 1 means that the directions are parallel (Figure 3A). Therefore, the similarity between ANM or ED modes and the transition is evaluated by the cumulative contribution of the first M modes

${\delta }_k= {\left({\displaystyle \sum _1^M{\alpha }_k^2}\right)}^{1/2}(4)$

Here, we consider M = 10, which typically cover >90% of the variance (Orellana et al., 2010) (Figure 2A).

Heuristic Channel-Defined Structural Variables

To evaluate channel descriptors as in Orellana et al. (2016) (Figures 5–7, Supplementary Table S3, S4), we used VMD along with in-house tools and scripts written in python, C++, and FORTRAN. The quaternary twist motion is the anti-correlated rotational movement of the ECD versus TMD around the channel axis, which decreases as the channel transitions to the open state. Channel closing is also coupled to ECD collapse (un-bloom) to a narrower diameter. The iris-like gating motion can be broken down into two components in the form of tilt and twist motions of the M2 helices that alter the pore radius. Hence, blooming was evaluated as the maximal radius of the extracellular domain defined by maximally separated residues at the tip of the five subunits. Quaternary twist was calculated as the average rotation angle of each subunit with the vector from extracellular domain and transmembrane domain CM to overall CM on the XY plane. The tilt and the twist angle of M2 helices were calculated as an average over five subunits. In order to get comparable angles between subunits, the reference structure (5CFB) was aligned to the center of mass of each M2 helix. X-axis was characterized between two centers of mass: M2 helix and the protein center of mass. The Z-axis was chosen as the principal of inertia parallel to the symmetry axis. Finally, the Y-axis was defined as the vector normal to the XZ plane. With these axes, the tilt angle is calculated between the projected helical axis onto the XZ plane and the Z-axis and twist angle between the projected helical axis onto the YZ plane and Z-axis. Pore calculations were performed with HOLE (Smart et al., 1996) using 0.5 Å step from C-alpha atoms only. The pore radius at 9′ was averaged over a 2.5 Å window in both directions from the center of mass of 9’ residues. The average hole profiles were obtained for the structures during the production run using only M2 helices and are reported with their standard deviation. Hydration in the GlyR channel was quantified by calculating the number of water molecules inside the channel pore along the axis normal to the bilayer. To study the conductivity of the GlyR channel in wild type and mutant proteins, we identified the water/ions in the pore and their z axis coordinates (bilayer normal) at every snapshot and plotted them versus simulation time.

Results and Discussion

GlyR Cancer Mutations and Dominant Hyperekplexia Changes Partially Overlap on TM1-2

As of March 2022, there are around 1000 cases of tumors reported to carry GlyR mutations: 952 for GLRA1, 1196 for GLRA2, and up to 1546 for GLRA3. After removing nonsense and deletion changes, we retrieved 245, 255, and 299 mutations, for GLRA1-3, respectively. Mutations for GLRA1 concentrate in 175 positions (70%), for GLRA2 in 174 (68%), and for GLRA3 in 208 (69%), which represent approximately 50–60% of the GlyR chain (347 positions); a majority of these changes are considered to be passenger mutations due to defective repair mechanisms in neoplastic cells. Nevertheless, as the mutations are not evenly distributed along GLRA sequences but preferentially focused on key conserved regions (Figures 1C,D, Supplementary Figure S1, Supplementary Table S1), we hypothesized they could impact on the receptor function. Particularly, for GLRA3, mutation S241L has high frequency, with over 20 cases reported, in contrast with top GLRA1 and GLRA2 mutations. Moreover, nearly half of the cases (12) of S241L mutation are linked to melanoma, followed by skin carcinoma (8), while in GLRA1-2, they tend to spread across multiple tumor types. A closer look also reveals specific changes recurrently appearing across GLRA1-3 genes and/or often overlapping with HPX mutations sometimes identical (see Figures 1C,D, Supplementary Table S1). Overall, these mutational patterns strongly suggest that at least a fraction of these changes could affect the channel function, which we decided to investigate in more detail for GLRA3, due to the availability of solved human X-ray structures and its significant mutational frequency. We focused our analysis on the resting state, represented by structure 5CFB. To filter out as many random changes as possible, we quantified mutation clustering (Figure 1D bottom, Supplementary Figure S2A) based on the 3D space distribution and degree of conservation as in Orellana et al. (2019b); we also mapped GLRA1-2 and hyperekplexia mutations to better detect overlapping patterns. We considered not only the raw number of mutations reported but also the recurrence of multiple allele changes onto the same positions (Supplementary Figure S2A). Upon weighting the spatial concentration of mutations by their ConSurf scores, we were able to identify four broad mutational clusters, which display overlap between hyperekplexia and cancer-reported GLRA1-3 changes: three minor ones at the ECD, located around loop D (residues 65–75), loop A (residues 95–110), and the Cys-loop ECD–TMD interface (residues 140–150), and a major one across TM1-2 (residues 222–244) to TM2 (residues 250–272); note that almost all reported HPX mutations overlap with one of these peaks with a few exceptions (Figure 1D bottom), with the notable exception of the most frequent mutation, R271L/Q/P, located at TM2-3 loop. Highest clustering scores were achieved at the area surrounding the TM1-2 loop region, where the selectivity filter and desensitization gate locate, and which contains HPX mutations such as P250T known to prolong desensitization (Saul et al., 1999). This area contains the highest frequency mutation, S241L/P (ConSurf score 8), located at the end of TM1, and its contacting neighbor R252S (ConSurf score 9), at the end of TM1-2 loop/beginning of TM2. While S241 is close to dominant hyperekplexia change W239C, of unknown functional effect, R252 is mutated in HPX as R252H/C, which results in low expression and activity. Another high-scoring neighbor area upstream TM2, close to 9′ gate, contains mutation V260M (ConSurf score 8), which corresponds exactly with a dominant HPX change (del Giudice et al., 2001; Castaldo et al., 2004) known to disrupt gating. Importantly, both R252 and V260 are located in the pore-opposing TM2 face, at the level of lower 2′ and upper 9′ gates (see Figures 1A,B, right). To further evaluate the significance of mutations, we also compared mutational patterns with the conserved contact networks of 5CFB (Figure 1D, top) and 5VDH (not shown). In contrast with S241 and R252, mostly involved in the local TM1-2 loop connectivity, V260 participates in longer-range intrachain contacts connecting TM2 with TM3; moreover, it also contacts TM2 (T262) across adjacent subunits in the closed but not the desensitized state (not shown). Overall, the concentration of these dramatic size, charge, and polarity changes overlapping HPX mutants near the conserved activation and desensitization gates suggests that they are not neutral but can impact receptor stability and conformation.

GlyR Ensemble Encodes Opening and Desensitization Transitions

In spite of the diversified roles of PLGICs, structural studies have revealed an astonishing degree of fold and conformational conservation across species. Our previous studies of the channel GLIC and other model proteins show that the PCA of structural ensembles containing multiple conformations can reveal the pathways for inter-connecting transitions providing an accurate framework to monitor MD sampling (Orellana et al., 2016; Orellana, 2019). Currently, there are dozens of zebrafish and human GlyR structures solved in closed, desensitized, and open and super-open states (Figure 2A). We aligned n = 33 nearly-intact homo- and heteromeric structures (Figure 2B; RMSD = 2.2 ± 0.7 Å) and performed PCA to extract the dominant ensemble motions and then investigate how they correlate with channel descriptors and annotated functional status. Similar to what is observed for prokaryotic GLIC (Orellana et al., 2016), the first mode (PC1, 63% of the variance) tracks the global quaternary twist and blooming motion (R = 0.9 and 0.6), as well as TM2 tilt (R = −0.8), while the second mode (PC2, 10% of variance) captures most of TM2 twist and pore gating (R = 0.7 and 0.6) (Figure 2B, Supplementary Figure S2B), separating super-open from desensitized structures (see Methods for definitions). Together, PC1-2 (73% variance) split the structural ensemble onto three to four major clusters (Figure 2C): to the right, closed/antagonist-bound structures (pore radius ≈1.5 Å, 1.2 Å for 5CFB), and to the left, those with an un-bloomed ECD and a wide-open (top left corner, pore radius > 3 Å e.g., 6PMO) or desensitized central channel (lower left corner, ≈ 1.5 Å, 1.4 Å for 5VDH); the only exception is the open-blocked structure 6UD3. Significantly, structures annotated as open states (6PLO, 6PM2, 6PLY, etc.) appear correctly located as a sub-cluster along the path from fully desensitized to super-open structures (pore radius ≈2.5 Å), which supports that PC1-2 space is a suitable framework to annotate channel status based on global correlated features of multiple descriptors. Finally, we also investigated to which extent these large-scale gating movements coupled to opening/desensitization are intrinsic to the different GlyR states, as previously suggested (Bertaccini et al., 2010; Zhu and Hummer, 2010; Zheng and Auerbach, 2011). We took as representatives of the three main PC1-2 clusters our reference human closed (5CFB) and desensitized (5VDH), together with one of the super-open zebrafish structures (6PMO); these three structures altogether broadly define three main transition directions along which the conformational cycle could potentially proceed. ENM from these three distinct states (Figure 3A) indicates that while the opening transition could be fairly spontaneous in the absence of a ligand (overlap ≈60%), the same was less likely for re/desensitization. In contrast, normal modes computed from 5VDH displayed poor overlap with the associated transitions, suggesting a rigid structure with more uncorrelated local motions. Residue fluctuations from 5CFB–ENM (Supplementary Figure S2C) also indicated that these mutations are located in a transition region from high [TM1-2 loop (S241L)] to low flexibility (TM2, from R252S to V260M) and thus could have diverse effects despite their proximity, which we then explored with fully atomistic MD.

Molecular Dynamics Reveals Mutations Perturb GlyR Conformational Dynamics

We performed MD simulations from the closed inhibited (PDB: 5CFB) and desensitized (PDB: 5VDH) human GLRA3 homopentamers after removing all ligands (the antagonist strychnine, in the first, and agonists glycine, AM-3607 and ivermectin in the latter), for WT and mutant sequences. Root-mean-square fluctuations (RMSDs) of the C_α atoms versus the initial experimental structures were calculated for all systems to confirm simulation convergence (Supplementary Figure S3A) and estimate the overall stability and flexibility of unbound closed and desensitized states upon ligand removal (Supplementary Figure S3B). Despite the fact that differences across replicates and GlyR variants were small, a closer look reveals a slightly more rigid desensitized state, with distributions shifted to lower values (Supplementary Figure S3B). This is often observed for bound/unbound systems and in agreement with our ENM preliminary analysis. Nevertheless, for both closed and desensitized state simulations, TM1-2 mutations in general shifted RMSD distributions to the right, suggesting increased conformational flexibility. While for R252S, the peak height is clearly shifted toward higher RMSDs suggesting greater thermal fluctuations, in V260M, the distribution spreads over multiple peaks suggestive of different conformational clusters. Interestingly, S241L displays virtually no difference versus WT GlyR for 5CFB simulations, in contrast with 5VDH simulations. Locally, the introduced mutations disrupt highly stable hydrogen bond and salt bridge WT interactions in both the closed and the desensitized receptors, both at the local TM1-2 level and globally (Supplementary Figures S4A,B). Significantly, a WT salt bridge between Asp247 and Arg252 present in 5CFB but absent in 5VDH is broken by all three mutations in the closed state. Among mutants, R252S has a major impact on the closed state introducing multiple non-native salt bridges not seen in any WT simulations far away from the mutation site, around TM3-4. Native contacts in this area, which are maintained in both 5CFB and 5VDH simulations, are the most perturbed by mutations; interestingly, V260M closed state simulations display interactions seen in the open-like/desensitized state like Glu300-Lys320. Similar changes are seen in hydrogen bond patterns, with mutations mostly affecting longer-range interactions with TM3 and TM4. Despite these profound reshaping of interactions, TM1-2 local helicity is mostly unaffected (Supplementary Figure S4C), only displaying transient and very limited unfolding in desensitized state simulations, focused at helix termini surrounding the TM1-2 loop junction.

Given the long-range impact of mutations in H-bonding and salt bridge connectivity, we then explored their effect on global dynamics by extracting the MD essential modes (ED; see Methods) explored by each system (Figure 3A) and computing their alignment with the three major transition directions defined from PC1-2 clusters. Consistent with ENM predictions, we observed that simulations from the closed state tend to spontaneously sample toward the direction of the super-open cluster (51% overlap with 5cfb ↔ 6pm0) more than toward the open-like desensitized area (overlap with 5cfb->5vdh, 0.45%), while simulations from 5VDH tend to relax toward the closed state and barely sample along the open-desensitized axis (76% overlap with 5cfb ↔ 5vdh versus 20%). By contrast, all mutants and specially V260M enhance the opening transition displaying even better alignment with 5cfb ↔ 6pm0 (55–60% versus 50%) while they are slightly less efficient sampling in the recovery direction toward the closed state (overlaps 5cfb ↔ 5vdh 70%), especially in the case of TM1-2 loop mutations. These trends are also visible upon MD projection onto PC1-2 to build the corresponding free energy landscapes (Figure 3B) and to examine trajectory time evolution (Figures 4A,B). While WT closed simulations sample one major minima around the starting structure and a minor one skewed toward the direction of super-open structures (Figure 3B top), mutants show one larger elongated minimum shifted again toward the same super-open direction. Inspection of the simulation PC1-2 time evolution reveals that indeed a fraction of 5CFB trajectories proceed along PC1(blooming) and PC2 (gating) toward un-blooming/pore opening directions characteristic of the open/desensitized state, especially in the case of V260M (Figure 4A). By contrast, WT desensitized simulations show two small minima in the direction of 5CFB closed state (at 0.0), while mutant ones, also sampling across the same PC1 direction, point toward lower PC2 values, which are only explored in closed mutant simulations (see below). Accordingly, all 5VDH trajectories uniformly proceed along PC1 toward 5cfb minima (Figure 4B), although in terms of PC2 pore gating, the majority evolve toward the super-open cluster. Overall, this suggests a general trend for mutants to favor channel opening while slowing down recovery from desensitization, as is often reported for HPX mutations. Interestingly, although the biological relevance of super-open structures has been questioned, projections onto PC1-2 space clearly show that all simulations from the closed state are sampling PC2 toward solved super-open structures, while simulations from the desensitized 6PMO sample mostly PC1 toward the closed state cluster.

FIGURE 4

FIGURE 4. Time evolution of MD replicates along PC1 and PC2. (A) Time evolution of MD trajectories along PC1, which separates closed from open/desensitized structures, for each of the four replicates from each system. Note how for simulations from the closed state (top), trajectories evolve in both closed and open/desensitized directions or stay around the starting structure, with the exception of V260M, in which all runs proceed toward the open region to the left. On the contrary, all desensitized state simulations (bottom) progress in the closing direction. (B) Time evolution of MD trajectories along PC2, which separates closed from open/desensitized structures, for each of the four replicates from each system. Trajectories stay or evolve toward pore closing, with the exception of S241L and V260M, which have replicates progressing toward opening. On contrast, most simulations from the desensitized state (bottom) evolve toward pore opening, except for R252S*.

Finally, to better characterize the conformational effect of mutations, we also examined the classical heuristic variables describing PLGIC geometry: quaternary twist, blooming, and TM2 tilt and twist angles (Supplementary Table S3, Figures 5, 6). As could be expected, the ECD quaternary twist (Figure 5A) was higher not only for closed versus open state simulations but also for mutants versus WT trajectories, displaying in all cases a neat gradient across PC1. However, less marked gradients were observed for the rest of channel parameters, as could be expected from their poorer correlations with PCs (Supplementary Figure S2B). Apart from the quaternary twist, the ECD configuration is also defined by its maximum radius or blooming, which was shifted to the right by all mutations in both closed and desensitized simulations (Figure 5B). Interestingly, ECD blooming appeared more constrained, i.e., un-bloomed (higher peaks) in TM2 R252S* and V260M*, in comparison with S241L*. The impact of mutations was also apparent in local TM2 variables like helix twist and tilt (Figures 6A,B), especially upon projection, which revealed heterogeneous distributions across the PC1-2 space. Globally speaking, TM2 twist and tilt increased along PC1 toward the open/desensitized state directions. Nevertheless, while R252S increased tilt values dramatically in desensitized state simulations but also in the closed state, reaching in both cases the highest values (10.5 and 9.9 Å, respectively), V260M displayed relatively low tilt angles, especially in 5VDH simulations. Notably, S241L displayed only a mild increase, only noticeable in 5VDH simulations. Overall, although TM1-2 mutations have a similar mild allosteric effect on ECD features, they profoundly and differently disrupt the local configuration of the pore, sampling extreme values for TM2 twist and tilt angles, which can appear uncorrelated to PC1 and associated blooming and quaternary twist.

FIGURE 5

FIGURE 5. WT and mutant GlyR global ECD dynamics from 1.2 us meta-trajectories. (A) Global ECD-TMD quaternary twist angle and (B) blooming. Note that the most global descriptor of channel status, quaternary twist, preserves a clear gradient-like pattern across the PC1-2 space both in WT and mutant simulations. See average values in Supplementary Table S3.

FIGURE 6

FIGURE 6. WT and mutant GlyR local TM2 features extracted from 1.2 us meta-trajectories. (A) TM2 helix twist and (B) TM2 helix tilt. Mutations perturb not only TM2 features but also display long-range allosteric effects on ECD blooming versus WT. See average values in Supplementary Table S3.

Pore Analysis Shows That TM2 Mutations Perturb Pore Gates and Water-Ion Permeation

As these heterogeneous TM2 mutant features suggested a major impact on the pore, we investigated in more detail their effect on pore radius and its permeability for water (Figure 7, Supplementary Table S3). As could be expected from trajectory time evolution (Figure 4B, bottom), WT desensitized simulations from 5VDH, which start with a closed radius of 1.48 Å, quickly evolve opening the pore as seen by broader radius distributions up to a maximum of 3.26 Å (1.7 ± 0.6 Å), in comparison with 5CFB simulations, which start at 1.28 Å and reach up to a maximum of 2.85 Å (1.3 ± 0.5 Å) (see Figure 7A, Supplementary Table S3). Significantly, mutations tend to display broader distributions, generally shifted to the right and with multiple peaks that suggests a diversity of open and closed pore configurations. Differences versus WT simulations are particularly significant for 5CFB simulations, with mutants displaying a clear tail (S241L, V260M) or secondary peak (R252S) at higher nearly open radius close to 2 Å (Figure 7A). This shift to the right is also seen for 5VDH mutant simulations with the exception of V260M; notably, secondary peaks at collapsed pore states barely sampled by the WT are also seen for mutants. These pore radius changes translate to variations in water permeation versus the WT protein. The number of water molecules in the channel ranges from 40 to 120 for the closed state GlyR protein, with an average of ≈80 water molecules in the channel. This seems to shift toward higher values and with a broader distribution for mutant proteins, which have an average of ≈90 water molecules in the channel and maximal values for all mutations around 130–140 (Supplementary Table S3). By contrast, simulations from 5VDH have on average less water molecules in the channel, although in the case of TM2 mutants can reach higher hydration than the WT (up to 170 water molecules for TM1-2 loop mutants versus 160 for WT-GlyR) (Figure 7B).

FIGURE 7

FIGURE 7. Minimum pore radius (A) and hydration (B) calculated for WT and mutant GlyR from MD simulations. While WT simulations display Gaussian-like distributions of both pore radius and hydration, mutations shift them toward the right and/or split them in multiple sub-peaks suggestive of enlarged and more hydrated pore states. See average values in Supplementary Table S3.

To gain more insight into these differently hydrated pore states explored by MD, we performed TMD cluster analysis on the complete trajectory set with GROMACS (Figure 8, Supplementary Table S4). Simulations initiated from the desensitized state sampled four clusters: cluster 1, which accounts for nearly half of the total conformations sampled by GlyR from closed and desensitized states, and clusters 3, 9, and 10, which in total, account for 7% of the total population and sample expanded pores up to 2Å. Interestingly, we observed that this major semi-open cluster 1 conformation (1.5 ± 0.4 Å) was also sampled by both V260M (13% of cluster 1) and R252S (4% of cluster 1) simulations initiated from the closed state, which allowed the observed increased water permeation (Figure 8A right). Cluster 2 was mostly sampled by WT closed state simulations (51%) as well as the three mutants studied. In contrast, cluster 3, which features an expanded pore around 1.7 Å, is distinctively assigned for R252S* (95% of the cluster). Clusters 5 and 7 are also unique of R252S closed state simulations and feature a slightly collapsed pore (0.8Å), also found in cluster 6, which has a mixed population of conformers from neighbor S241L and R252S mutants (see Supplementary Table S4). Notably, similar collapsed states, previously observed in simulations (Dämgen and Biggin, 2020), are characterized by loss of the fivefold symmetry of the TMD pentamer (cluster 5 in Figures 8B,C). In our simulations, these states accounted, however, for a minor fraction of total sampling and were mostly populated by R252S, which, on the other hand, sampled also expanded pore states like the other mutants as well as WT* trajectories.

FIGURE 8

FIGURE 8. C-alpha RMSD cluster analysis of the GlyR TMD. (A) Projections of major clusters onto PC1-2 space (left) and cluster representatives showing water permeation (right): top row, cluster 2 (WT) versus cluster 1 (WT* and V260M); bottom row, cluster 3 (R252S*), cluster 4 (S241L), and cluster 5 (R252S). (B) Top and bottom views of representative TMD clusters and M2 helix residues shown with surface representation. Note the similarity of the open state from 5VDH in WT simulations (WT*) versus the V260M initiated from 5CFB (V260M). (C) Top and bottom views of representative TMD clusters, with side chains for pore gates (Leu261, Pro250) shown as licorice. See Supplementary Table S4.

From PLGIC structures and mutagenesis work, it is known that there are two major constriction points in the GlyR pore: an upper ring of hydrophobic residues at the 9′ position of TM2 (the L9′ activation gate, formed by Leu261) and a lower ring of Pro residues at −2′ (the P-2′ desensitization gate, corresponding to Pro250). While the central L9′ gate is closed in antagonist-bound structures (5CFB, 3JAD), in desensitized structures (5VDH, 3JAF) it is P-2′ that closes the intracellular “mouth” of the channel. The P-2′ is wide-open in the so-called super-open structures (3JAE). In general, despite the observed minoritary collapsed TMD clusters, all our simulations resulted in pore relaxation versus the initial crystallographic structures, but with substantial differences upon mutation. A closer inspection of pore profiles (Figure 9A) revealed that TM1-2 mutants dramatically reshape the central channel in a similar way. In the closed state, mutations expand the lower (P250, 9′) but decrease the upper (A272, 20’) pore sections, while in the desensitized simulations, they tend to destabilize and enlarge both. Mutations also displayed wider pore fluctuations, remarkable for S241L and V260M. Although the short time scale of our simulations did not allow us to observe complete opening of the channel, some of the sampled clusters, especially for the mutants, reached pore radius around 1.7–1.8 Å, which can allow the passage of a partially hydrated chloride ion (1.8Å); in such wider pore conformations, water was observed forming continuous molecule chains across the channels. A closer examination of water and ion penetration into the channel revealed that chloride ions penetrate indeed deep into the channel from the enlarged intracellular “mouth”, although only for WT* and V260M* simulations, this resulted in complete permeation events (Figures 9B,C, Supplementary Movies S1, S2). Notably, the enlarged open mouth at the level of the selectivity filter also allowed the entry of potassium ions to the vestibule in 5CFB simulations.

FIGURE 9

FIGURE 9. Pore radius profile for WT and mutant GlyR from MD simulations. (A) Pore profiles for closed/resting state WT (5CFB) and mutants constructed from it are shown to the left; WT* (5VDH) and mutants constructed from it are shown to the right side. Note how in both closed and desensitized/open state simulations all mutants tend to widen the P250 gate. (B) Water and ion permeation across the channel. The closed activation gate in 5CFB simulations (top) appears as a white area around Leu261, missing in desensitized state simulations (bottom). Chloride ions cross the area in WT* and V260M* in relation with permeation events (see Supplementary Movies S1, S2); there is also an entry of potassium ions through the selectivity filter (P250) in closed state simulations (yellow), but they stay in the vestibule. (C) Chloride ion distribution inside the channels; note the shift to the right for mutants due to the wider channels accommodating a larger number of ions.

Summary and Conclusion

Here, we aim to perform a preliminary study of uncharacterized GlyR mutations reported in human tumors. Based on 3D-clustering, conservation, and vicinity or overlap with hyperekplexia changes, we focused on three mutations potentially disrupting gating for further study: one located in the TM1-2 loop, at the “mouth” of the channel (S241L), and two at the back of pore TM2 (R252S and V260M). Our goal was to explore whether they were neutral or perturbed channel dynamics and potentially function. To evaluate functional/conformational status and monitor our simulations, we performed PCA of the current ensemble of GlyR structures, in order to obtain suitable reaction coordinates to build an “experimental” conformational landscape (Orellana et al., 2019a; Orellana, 2019). We confirmed the resulting PCs tracked multiple key channel descriptors providing an automatic and reasonable classification of the solved structures. MD projections on PC1-2 space allowed us to quantify how trajectories sample large-scale motions correlated with multiple channel features. This provided a framework to analyze sampling by MD simulations of human WT and mutant GlyR. Surprisingly, we found evidence that not only the closed-super-open transition is encoded in the experimental structures but also that they are indeed sampled in unbiased WT simulations in the absence of ligand. While WT simulations from the closed state clearly evolve in the direction of super-open solved structures, simulations from the desensitized state tend to sample the recovery transition back to the closed state and not toward the super-open, as suggested by ENM. Our integrated structural analysis also revealed that mutations perturb sampling in a similar way, enhancing the exploration of the opening direction and disrupting relaxation to the native closed state. These global features relate to the local disruption of TM1-2 contacts, which results in an altered TM2 orientation that favors a wider pore and increased water permeation. Nevertheless, the impact of mutations was clearly distinct, with the mutation R252S, which targets the most conserved residue, resulting in non-native interactions and sampling of collapsed pore states. On the contrary, mutation V260M was the closest in behavior to the closed state and displayed enhanced sampling of the opening transition. Importantly, chloride permeation was only observed for WT* and V260M* simulations from the desensitized state. Functional evidence suggests that the mutation of Arg252 results in non-conducting channels, which would agree with our observations of a profoundly impaired channel. On contrast, the V260M HPX mutant has been shown to result in spontaneous activity and prolonged desensitization, also in agreement with our simulations. Functional data for mutations surrounding S241L are still lacking, but our results suggest a behavior in between the other two, probably closer to WT and V260M in terms of channel permeability and conduction despite its proximity to R252S.

In summary, our analysis indicates, on one side, that the concentration of GLRA1-3 mutations across TM1 and TM2 is very likely not random, and, on the other, that these changes certainly disrupt the GlyR function in the tumor cells carrying them, as seen for HPX. Whether and how disrupting GlyR function can represent and advantage for cancer cells is beyond the scope of this work but certainly deserves further computational as well as functional investigation.

Data Availability Statement

The data and code used for this study are available from the corresponding author LO upon request.

Author Contributions

AM performed MD simulations and analysis, and prepared figures under LO supervision. OY provided scripts for channel analysis and contributed to discussions. LO conceived the original idea, performed mutation clustering, ENM, PCA, and ED analysis, designed MD simulations, interpreted results, prepared figures and wrote the manuscript. All authors participated in data interpretation and manuscript revision.

Funding

This work was supported by Cancerfonden Junior Investigator Award (21 0305 JIA 01 H) and Project Grant (CF 21 1471 Pj 2022-24), Jeanssons Foundation 2021, Vetenskapsrådet Starting Grant VR2021-02248 2022-25, and Karolinska Institute KIS, KIF Grants to LO. Simulations were run using the Swedish National Infrastructure for Computing (allocations Q13 SNIC 2021/5-87 and 2021/5-409 to LO).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

LO thanks the Swedish Cancer Foundation (Cancerfonden), the Swedish Research Council (Vetenskapsrådet), the Jeanssons Foundation (Harald och Greta Jeanssons Stiftelse) and Karolinska Institute for generous support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.890851/full#supplementary-material

References

Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 1–2, 19–25. doi:10.1016/j.softx.2015.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Amadei, A., Linssen, A. B. M., and Berendsen, H. J. C. (1993). Essential Dynamics of Proteins. Proteins 17, 412–425. doi:10.1002/prot.340170408

PubMed Abstract | CrossRef Full Text | Google Scholar

Amadei, A., Linssen, A. B. M., de Groot, B. L., Berendsen, H. J. C., and Berendsen, H. J. (1996). An Efficient Method for Sampling the Essential Subspace of Proteins. J. Biomol. Struct. Dyn. 13, 615–625. doi:10.1080/07391102.1996.10508874

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashkenazy, H., Abadi, S., Martz, E., Chay, O., Mayrose, I., Pupko, T., et al. (2016). ConSurf 2016: An Improved Methodology to Estimate and Visualize Evolutionary Conservation in Macromolecules. Nucleic Acids Res. 44 (W1), W344–W350. doi:10.1093/nar/gkw408

PubMed Abstract | CrossRef Full Text | Google Scholar

Atilgan, A. R., Durell, S. R., Jernigan, R. L., Demirel, M. C., Keskin, O., and Bahar, I. (2001). Anisotropy of Fluctuation Dynamics of Proteins with an Elastic Network Model. Biophysical J. 80, 505–515. doi:10.1016/S0006-3495(01)76033-X

CrossRef Full Text | Google Scholar

Atilgan, C. (2018). Computational Methods for Efficient Sampling of Protein Landscapes and Disclosing Allosteric Regions. 1st ed. Elsevier. doi:10.1016/bs.apcsb.2018.06.001

CrossRef Full Text | Google Scholar

Bahar, I., Lezon, T. R., Bakan, A., and Shrivastava, I. H. (2010). Normal Mode Analysis of Biomolecular Structures: Functional Mechanisms of Membrane Proteins. Chem. Rev. 110, 1463–1497. doi:10.1021/cr900095e

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertaccini, E. J., Trudell, J. R., and Lindahl, E. (2010). Normal Mode Gating Motions of a Ligand-Gated Ion Channel Persist in a Fully Hydrated Lipid Bilayer Model. ACS Chem. Neurosci. 1, 552. doi:10.1021/cn100026t

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, D., Gawali, V. S., Kallay, L., Toukam, D. K., Koehler, A., Stambrook, P., et al. (2021). Therapeutically Leveraging GABAA Receptors in Cancer. Exp. Biol. Med. 246, 2128. doi:10.1177/15353702211032549

CrossRef Full Text | Google Scholar

Bode, A., Wood, S.-E., Mullins, J. G. L., Keramidas, A. T. D., Cushion, R. H., Pickrell, W. O., et al. (2013). New Hyperekplexia Mutations Provide Insight into Glycine Receptor Assembly, Trafficking, and Activation Mechanisms. J. Biol. Chem. 288, 33745–33759. doi:10.1074/jbc.M113.509240

PubMed Abstract | CrossRef Full Text | Google Scholar

Bode, A., and Lynch, J. W. (2013). Analysis of Hyperekplexia Mutations Identifies Transmembrane Domain Rearrangements that Mediate glycine Receptor Activation. J. Biol. Chem. 288, 33760–33771. doi:10.1074/jbc.M113.513804

PubMed Abstract | CrossRef Full Text | Google Scholar

Bode, A., and Lynch, J. W. (2014). The Impact of Human Hyperekplexia Mutations on glycine Receptor Structure and Function. Mol. Brain 7, 2. doi:10.1186/1756-6606-7-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Calimet, N., Simoes, M., Changeux, J.-P., Karplus, M., Taly, A., and Cecchini, M. (2013). A Gating Mechanism of Pentameric Ligand-Gated Ion Channels. Proc. Natl. Acad. Sci. 110, E3987–E3996. doi:10.1073/pnas.1313785110

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A. (1999). Molecular Dynamics and Normal Mode Analysis of Biomolecular Rigidity. Mol. Biol. 18, 329–344. doi:10.1007/0-306-47089-6

CrossRef Full Text | Google Scholar

Castaldo, P., Stefanoni, P., Miceli, F., Coppola, G., del Giudice, E. M., Bellini, G., et al. (2004). A Novel Hyperekplexia-Causing Mutation in the Pre-transmembrane Segment 1 of the Human glycine Receptor α1 Subunit Reduces Membrane Expression and Impairs Gating by Agonists. J. Biol. Chem. 279, 1200. doi:10.1074/jbc.M311021200

CrossRef Full Text | Google Scholar

Cerdan, A. H., and Cecchini, M. (2020). On the Functional Annotation of Open-Channel Structures in the Glycine Receptor. Structure 28, 690. doi:10.1016/j.str.2020.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerdan, A. H., Martin, N. É., and Cecchini, M. (2018). An Ion-Permeable State of the Glycine Receptor Captured by Molecular Dynamics. Structure 26, 1555. doi:10.1016/j.str.2018.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Changeux, J.-P. (2014). Protein Dynamics and the Allosteric Transitions of Pentameric Receptor Channels. Biophys. Rev. 6, 311–321. doi:10.1007/s12551-014-0149-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Wang, L., Wang, W., Sun, H., Pang, L., and Bao, H. (2021a). Conformational Transformation of Switch Domains in GDP/K-Ras Induced by G13 Mutants: An Investigation through Gaussian Accelerated Molecular Dynamics Simulations and Principal Component Analysis. Comput. Biol. Med. 135, 4639. doi:10.1016/j.compbiomed.2021.104639

CrossRef Full Text | Google Scholar

Chen, J., Zhang, S., Wang, W., Pang, L., Zhang, Q., and Liu, X. (2021b). Mutation-Induced Impacts on the Switch Transformations of the GDP-And GTP-Bound K-Ras: Insights from Multiple Replica Gaussian Accelerated Molecular Dynamics and Free Energy Analysis. J. Chem. Inf. Model. 61, 1470. doi:10.1021/acs.jcim.0c01470

CrossRef Full Text | Google Scholar

Chung, S.-K., Vanbellinghen, J.-F., Mullins, J. G. L., Robinson, A., Hantke, J., Hammond, C. L., et al. (2010). Pathophysiological Mechanisms of Dominant and Recessive GLRA1 Mutations in Hyperekplexia. J. Neurosci. 30, 9612–9620. doi:10.1523/jneurosci.1763-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Dämgen, M. A., and Biggin, P. C. (2020). A Refined Open State of the Glycine Receptor Obtained via Molecular Dynamics Simulations. Structure 28(1). 130. doi:10.1016/j.str.2019.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Dacosta, C. J. B., and Baenziger, J. E. (2013). Gating of Pentameric Ligand-Gated Ion Channels: Structural Insights and Ambiguities. Structure 21, 1271–1283. doi:10.1016/j.str.2013.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Daidone, I., and Amadei, A. (2012). Essential Dynamics: Foundation and Applications. WIREs Comput. Mol. Sci. 2, 762–770. doi:10.1002/wcms.1099

CrossRef Full Text | Google Scholar

Dämgen, M. A., Zaki, A. M., and Biggin, P. C. (2020). Comment on “On the Functional Annotation of Open-Channel Structures in the Glycine Receptor. Structure 28 (6), 601–603. doi:10.1016/j.str.2020.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Giudice, E. M., Coppola, G., Bellini, G., Cirillo, G., Scuccimarra, G., and Pascotto, A. (2001). A Mutation (V260M) in the Middle of the M2 Pore-Lining Domain of the glycine Receptor Causes Hereditary Hyperekplexia. Eur. J. Hum. Genet. 9 (11), 873–876. doi:10.1038/sj.ejhg.5200729

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, J., Lü, W., Wu, S., Cheng, Y., and Gouaux, E. (2015). Glycine Receptor Mechanism Elucidated by Electron Cryo-Microscopy. Nature 526 (7572), 224–229. doi:10.1038/nature14853

PubMed Abstract | CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 103, 8577–8593. doi:10.1063/1.470117

CrossRef Full Text | Google Scholar

Ewald, P. P. (1921). Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. 369 (3), 253–287. doi:10.1002/andp.19213690304

CrossRef Full Text | Google Scholar

Gielen, M., and Corringer, P.-J. (2018). The Dual-Gate Model for Pentameric Ligand-Gated Ion Channels Activation and Desensitization. J. Physiol. 596, 1873–1902. doi:10.1113/JP275100

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 18 (12), 1463–1472. doi:10.1002/(SICI)1096-987X(199709)18

CrossRef Full Text | Google Scholar

Hibbs, R. E., and Gouaux, E. (2011). Principles of Activation and Permeation in an Anion-Selective Cys-Loop Receptor. Nature 474, 54–60. doi:10.1038/nature10139

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoover, W. G. (1985). Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 31 (3), 1695–1697. doi:10.1103/PhysRevA.31.1695

PubMed Abstract | CrossRef Full Text | Google Scholar

Howard, R. J. (2021). Elephants in the Dark: Insights and Incongruities in Pentameric Ligand-Gated Ion Channel Models. J. Mol. Biol. 433 (17), 167128. doi:10.1016/j.jmb.2021.167128

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., De Groot, B. L., et al. (2016). CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 14, 4067. doi:10.1038/nmeth.4067

CrossRef Full Text | Google Scholar

Huang, X., Chen, H., Michelsen, K., Schneider, S., and Shaffer, P. L. (2015). Crystal Structure of Human glycine Receptor-Α3 Bound to Antagonist Strychnine. Nature 526 (7572), 277–280. doi:10.1038/nature14972

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Chen, H., and Shaffer, P. L. (2017a). Crystal Structures of Human GlyRα3 Bound to Ivermectin. Structure 25, 945–950.e2. doi:10.1016/j.str.2017.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Shaffer, P. L., Ayube, S., Bregman, H., Chen, H., Lehto, S. G., et al. (2017b). Crystal Structures of Human glycine Receptor α3 Bound to a Novel Class of Analgesic Potentiators. Nat. Struct. Mol. Biol. 24, 108–113. doi:10.1038/nsmb.3329

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, S., Cheng, X., Lee, J., Kim, S., Park, S. J., Patel, D. S., et al. (2017). CHARMM-GUI 10 Years for Biomolecular Modeling and Simulation. J. Comput. Chem. 38, 1114–1124. doi:10.1002/jcc.24660

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, S., Kim, T., Iyer, V. G., and Im, W. (2008). CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. J. Comput. Chem. 29, 1859–1865. doi:10.1002/jcc.20945

PubMed Abstract | CrossRef Full Text | Google Scholar

Jolliffe, I. T. (2002). Principal Component Analysis. doi:10.1007/b98835

CrossRef Full Text | Google Scholar

Kumar, A., Basak, S., Rao, S., Gicheru, Y., Mayer, M. L., Sansom, M. S. P., et al. (2020). Mechanisms of Activation and Desensitization of Full-Length glycine Receptor in Lipid Nanodiscs. Nat. Commun. 11, 17364. doi:10.1038/s41467-020-17364-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, T. M., Sivilotti, L. G., Colquhoun, D., Gardiner, R. M., Schoepfer, R., and Rees, M. (1998). Properties of Human glycine Receptors Containing the Hyperekplexia Mutation α1(K276E), Expressed inXenopusoocytes. J. Physiol. 507, 25–40. doi:10.1111/j.1469-7793.1998.025bu.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, J. W. (2004). Molecular Structure and Function of the glycine Receptor Chloride Channel. Physiol. Rev. 84, 1051. doi:10.1152/physrev.00042.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, J. W., Zhang, Y., Talwar, S., and Estrada-Mondragon, A. (2017). Glycine Receptor Drug Discovery. 1st ed. Elsevier. doi:10.1016/bs.apha.2017.01.003

CrossRef Full Text | Google Scholar

Martin, N. E., Malik, S., Calimet, N., Changeux, J.-P., and Cecchini, M. (2017). Un-Gating and Allosteric Modulation of a Pentameric Ligand-Gated Ion Channel Captured by Molecular Dynamics. PLoS Comput. Biol. 13 (10), e1005784. doi:10.1371/journal.pcbi.1005784

PubMed Abstract | CrossRef Full Text | Google Scholar

Monod, J., Wyman, J., and Changeux, J.-P. (1965). On the Nature of Allosteric Transitions: a Plausible Model. J. Mol. Biol. 12, 88–118. doi:10.1016/S0022-2836(65)80285-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Nemecz, Á., Prevost, M. S., Menny, A., and Corringer, P.-J. (2016). Emerging Molecular Mechanisms of Signal Transduction in Pentameric Ligand-Gated Ion Channels. Neuron 90, 452–470. doi:10.1016/j.neuron.2016.03.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Nosé, S. (1984). A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 52, 255. doi:10.1080/00268978400101201

CrossRef Full Text | Google Scholar

Orellana, L., Gustavsson, J., Bergh, C., Yoluk, O., and Lindahl, E. (2019a). EBDIMS Server: Protein Transition Pathways with Ensemble Analysis in 2D-Motion Spaces. Bioinformatics 35, 3505–3507. doi:10.1093/bioinformatics/btz104

PubMed Abstract | CrossRef Full Text | Google Scholar

Orellana, L. (2019). Large-Scale Conformational Changes and Protein Function: Breaking the In Silico Barrier. Front. Mol. Biosci. 6, 117. doi:10.3389/fmolb.2019.00117

PubMed Abstract | CrossRef Full Text | Google Scholar

Orellana, L., Rueda, M., Ferrer-Costa, C., Lopez-Blanco, J. R., Chacón, P., and Orozco, M. (2010). Approaching Elastic Network Models to Molecular Dynamics Flexibility. J. Chem. Theory Comput. 6, 2910–2923. doi:10.1021/ct100208e

PubMed Abstract | CrossRef Full Text | Google Scholar

Orellana, L., Thorne, A. H., Lema, R., Gustavsson, J., Parisian, A., Hospital, A., et al. (2019b). Oncogenic Mutations at the EGFR Ectodomain Structurally Converge to Remove a Steric Hindrance on a Kinase-Coupled Cryptic Epitope. Proc. Natl. Acad. Sci. U. S. A. 116, 110009. doi:10.1073/pnas.1821442116

CrossRef Full Text | Google Scholar

Orellana, L., Yoluk, O., Carrillo, O., Orozco, M., and Lindahl, E. (2016). Prediction and Validation of Protein Intermediate States from Structurally Rich Ensembles and Coarse-Grained Simulations. Nat. Commun. 7, 12575. doi:10.1038/ncomms12575

PubMed Abstract | CrossRef Full Text | Google Scholar

Parrinello, M., and Rahman, A. (1981). Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 52. doi:10.1063/1.328693

CrossRef Full Text | Google Scholar

Prevost, M. S., Sauguet, L., Nury, H., van Renterghem, C., Huon, C., Poitevin, F., et al. (2012). A Locally Closed Conformation of a Bacterial Pentameric Proton-Gated Ion Channel. Nat. Struct. Mol. Biol. 19, 642–649. doi:10.1038/nsmb.2307

PubMed Abstract | CrossRef Full Text | Google Scholar

Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., et al. (2013). GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 29, 845–854. doi:10.1093/bioinformatics/btt055

PubMed Abstract | CrossRef Full Text | Google Scholar

Saul, B., Kuner, T., Sobetzko, D., Brune, W., Hanefeld, F., Meinck, H. M., et al. (1999). Novel GLRA1 Missense Mutation (P250T) in Dominant Hyperekplexia Defines an Intracellular Determinant of glycine Receptor Channel Gating. J. Neurosci. 19, 869. doi:10.1523/jneurosci.19-03-00869.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

Smart, O. S., Neduvelil, J. G., Wang, X., Wallace, B. A., and Sansom, M. S. P. (1996). HOLE: A Program for the Analysis of the Pore Dimensions of Ion Channel Structural Models. J. Mol. Graph. 14. doi:10.1016/S0263-7855(97)00009-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Taly, A., Hénin, J., Changeux, J.-P., and Cecchini, M. (2014). Allosteric Regulation of Pentameric Ligand-Gated Ion Channels: An Emerging Mechanistic Perspective. Channels 8, 350–360. doi:10.4161/chan.29444

PubMed Abstract | CrossRef Full Text | Google Scholar

Tate, J. G., Bamford, S., Jubb, H. C., Sondka, Z., Beare, D. M., Bindal, N., et al. (2019). COSMIC: The Catalogue of Somatic Mutations in Cancer. Nucleic Acids Res. 47, D941. doi:10.1093/nar/gky1015

PubMed Abstract | CrossRef Full Text | Google Scholar

The UniProt Consortium (2021). UniProt: the Universal Protein Knowledgebase in 2021. Nucleic Acids Res. 49, D480–D489. doi:10.1093/nar/gkaa1100

PubMed Abstract | CrossRef Full Text | Google Scholar

Tirion, M. M. (1996). Large Amplitude Elastic Motions in Proteins from a Single-Parameter, Atomic Analysis. Phys. Rev. Lett. 77, 1905–1908. doi:10.1103/physrevlett.77.1905

PubMed Abstract | CrossRef Full Text | Google Scholar

Van den Eynden, J., Ali, S. S., Horwood, N., Carmans, S., Brône, B., Hellings, N., et al. (2009). Glycine and glycine Receptor Signalling in Non-neuronal Cells. Front. Mol. Neurosci. 2, 9–12. doi:10.3389/neuro.02.009.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, E. L., Cheng, X., Jo, S., Rui, H., Song, K. C., Dávila-Contreras, E. M., et al. (2014). CHARMM-GUI Membrane Builder toward Realistic Biological Membrane Simulations. J. Comput. Chem. 35, 23702. doi:10.1002/jcc.23702

CrossRef Full Text | Google Scholar

Yang, L.-W., Eyal, E., Bahar, I., and Kitao, A. (2009). Principal Component Analysis of Native Ensembles of Biomolecular Structures (PCA_NEST): Insights into Functional Dynamics. Bioinformatics 25, 606–614. doi:10.1093/bioinformatics/btp023

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, S. Z., and Bordey, A. (2009). GABA’s Control of Stem and Cancer Cell Proliferation in Adult Neural and Peripheral Niches. Physiology 24, 171. doi:10.1152/physiol.00002.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, H., Bai, X. C., and Wang, W. (2021a). Characterization of the Subunit Composition and Structure of Adult Human glycine Receptors. Neuron 109, 2707. doi:10.1016/j.neuron.2021.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Zhu, H., Lape, R., Greiner, T., Du, J., Lü, W., et al. (2021b). Mechanism of Gating and Partial Agonist Action in the glycine Receptor. Cell 184, 26. doi:10.1016/j.cell.2021.01.026

CrossRef Full Text | Google Scholar

Zhang, X., Zhang, R., Zheng, Y., Shen, J., Xiao, D., Li, J., et al. (2013). Expression of Gamma-Aminobutyric Acid Receptors on Neoplastic Growth and Prediction of Prognosis in Non-small Cell Lung Cancer. J. Transl. Med. 11, 102. doi:10.1186/1479-5876-11-102

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, W., and Auerbach, A. (2011). Decrypting the Sequence of Structural Events during the Gating Transition of Pentameric Ligand-Gated Ion Channels Based on an Interpolated Elastic Network Model. PLoS Comput. Biol. 7,e001046. doi:10.1371/journal.pcbi.1001046

CrossRef Full Text | Google Scholar

Zhu, F., and Hummer, G. (2010). Pore Opening and Closing of a Pentameric Ligand-Gated Ion Channel. Proc. Natl. Acad. Sci. U.S.A. 107, 19814–19819. doi:10.1073/pnas.1009313107

PubMed Abstract | CrossRef Full Text | Google Scholar