1. Computational and Systems Biology
  2. Structural Biology and Molecular Biophysics
Download icon

Markov state models of proton- and pore-dependent activation in a pentameric ligand-gated ion channel

  1. Cathrine Bergh
  2. Stephanie A Heusser
  3. Rebecca Howard
  4. Erik Lindahl  Is a corresponding author
  1. Science for Life Laboratory and Swedish e-Science Research Center, Department of Applied Physics, KTH Royal Institute of Technology, Sweden
  2. Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University, Sweden
Research Article
  • Cited 0
  • Views 489
  • Annotations
Cite this article as: eLife 2021;10:e68369 doi: 10.7554/eLife.68369

Abstract

Ligand-gated ion channels conduct currents in response to chemical stimuli, mediating electrochemical signaling in neurons and other excitable cells. For many channels, the details of gating remain unclear, partly due to limited structural data and simulation timescales. Here, we used enhanced sampling to simulate the pH-gated channel GLIC, and construct Markov state models (MSMs) of gating. Consistent with new functional recordings, we report in oocytes, our analysis revealed differential effects of protonation and mutation on free-energy wells. Clustering of closed- versus open-like states enabled estimation of open probabilities and transition rates, while higher-order clustering affirmed conformational trends in gating. Furthermore, our models uncovered state- and protonation-dependent symmetrization. This demonstrates the applicability of MSMs to map energetic and conformational transitions between ion-channel functional states, and how they reproduce shifts upon activation or mutation, with implications for modeling neuronal function and developing state-selective drugs.

Introduction

The family of pentameric ligand-gated ion channels (pLGICs), also known as Cys-loop receptors, controls electrochemical signal transduction in numerous tissues and cell types, from bacteria to humans. A rapid cycling between conducting and nonconducting conformations in response to chemical stimuli, such as neurotransmitter binding or changes in pH, is fundamental to their function. These channels are often found in the postsynaptic membrane of neurons and undergo allosteric conformational changes, where the pore in the integral transmembrane domain (TMD) opens for ion conduction upon selective neurotransmitter binding in the extracellular domain (ECD). Prokaryotic homologs, such as the proton-gated channel GLIC from the cyanobacterium Gloeobacter violaceus, share many topological features with eukaryotic pLGICs, and have been proposed to follow similar gating pathways (Figure 1A–C). For any pLGIC, channel function can ultimately be explained by its free energy landscape, where different state populations shift upon activation through agonist binding. Understanding these landscapes and how the free energies either of stable states or barriers change upon ligand-binding or mutation is thus crucial for full understanding of the gating process, with applications including the development of state-dependent drugs for better treatment of diseases related to channel malfunction.

Global architecture of GLIC, electrophysiology data and computational methodology.

GLIC in an open conformation shown from (A) the top and (B) the side in a POPC lipid bilayer. (C) Two opposing subunits highlighting the pore of the channel. Light colors represent the open conformation (PDB ID 4HFI) and dark colors the closed one (PDB ID 4NPQ). Arrows indicate important gating motions - the tilting of the M2 helices, beta expansion and ECD spread. Residues I233 and H235 on the pore-lining M2 helices were mutated in both simulations and electrophysiology experiments. (D) Electrophysiology data for wild-type GLIC, the gain-of-function I233T variant and the loss-of-function H235Q variant. (E) Simulation methodology: the eBDIMS method provides (1) coarse-grained seed structures along the transition pathway followed by (2) reconstruction of the atomistic detail. Atomistic structures were then (3) embedded in lipid bilayers and massively parallel unrestrained MD launched. Analysis involves (4) a feature transformation to account for the symmetry of the pentamer, followed by (5) dimensionality reduction with tICA and (6) MSM construction.

Recent advances in structural biology have enabled a steady increase in the number of available pLGIC structures. The GLIC model system is notable in this regard, accounting for over 40% of pLGIC entries in the protein data bank, including apparent closed and open states. However, the conformational diversity of these states is highly limited, leading to crude representations of the free-energy landscapes from experimental structures alone. Computational methods like molecular dynamics (MD) simulations can be used to sample more of the conformational landscape, and several studies have been conducted on GLIC to study short-timescale motions; such as simulations of the transmembrane domain only (Zhu and Hummer, 2010; Zhu and Hummer, 2012a; Zhu and Hummer, 2012b), studies of the ion permeation pathway through potential-of-mean-force calculations (Cheng et al., 2010; Fritsch et al., 2011), and steady-state simulations reaching 100 ns to 1 μs timescales (Nury et al., 2010; Prevost et al., 2012; Calimet et al., 2013), some also with additional ligands or modulations (Brannigan et al., 2010; Willenbring et al., 2011; LeBard et al., 2012; Murail et al., 2012; Laurent et al., 2016; Heusser et al., 2018; Faulkner and de Leeuw, 2020). Still, due to the large system size and relatively long timescales of the gating transitions, in practice it has not been feasible to sample complete gating transitions, especially if ligand-binding and unbinding events are involved (Chakrapani and Auerbach, 2005; Velisetty and Chakrapani, 2012; Gonzalez-Gutierrez et al., 2012; Laha et al., 2013; Gonzalez-Gutierrez et al., 2013; Velisetty et al., 2012; Menny et al., 2017). To bridge this gap, various enhanced sampling methods can be used, often by the application of biasing forces or presumed reaction coordinates. For instance, application of the string method with swarms of trajectories recently enabled the identification of local rearrangements in channel closure, including contraction of the upper pore, loosening of β-strand contacts in the lower ECD, and general expansion of the upper ECD (Lev et al., 2017). This provides precious information of structural rearrangements, but the choice of collective variables in combination with the timescales of individual simulations may influence what motions are sampled.

Markov state models (MSMs) show promise in modeling the thermodynamics and kinetics of biological systems without making prior mechanistic assumptions (Dämgen and Biggin, 2019; Prinz et al., 2011; Husic and Pande, 2018). By counting transitions, multiple simulations can effectively be stitched together to capture processes on timescales longer than any individual simulation, and probing more biologically relevant dynamics. However, in principle MSMs first require the whole equilibrium distribution to be sampled, which is difficult to achieve by starting simulations from experimental structures alone, since the stability required for crystal packing or cryo-EM data processing results in structures sometimes representing metastable states. Simulations started from these states then usually remain confined to the energy well for long time periods. In contrast, seeding simulations at regions in the free energy landscape that are not necessarily metastable enhance sampling without introducing any biasing forces, and the actual sampling is performed without limiting the system to any particular reaction coordinate.

Here, we have used such enhanced seeding approaches combined with MSMs to sample the GLIC opening-closing transition. Both the wild-type and two different variants with mutations along the pore-lining M2 helices - where we also showed one to yield gain-of-function similar to human homologs (Filatov and White, 1995; Kosolapov et al., 2000) and the other loss of function (Fourati et al., 2018) - were simulated in resting or activating conditions (neutral or low pH). From the resulting trajectories we were able to construct MSMs of the complete ion channel, quantitatively map free-energy landscapes, extract molecular details regarding the gating mechanisms, and show how the gating shift observed in electrophysiology recordings of the wild-type and variants is modeled correctly by the MSMs. Additionally, we present new evidence for the role of symmetry in gating.

Results

Enhanced sampling enables MSM construction of GLIC gating

To shed light on the gating mechanisms, we combined enhanced sampling of MD simulations at resting and activating pH with Markov state modeling. X-ray structures of GLIC crystallized at pH 7 (PDB ID 4NPQ) and pH 4 (PDB ID 4HFI) have been reported to represent closed and open states, respectively (Sauguet et al., 2014; Figure 1C). Because unbiased molecular dynamics simulations started from these states alone were not expected to thoroughly sample the activation process, we seeded simulations along the presumed gating transition which enabled simulations to run in a massively parallel fashion. To achieve this, the starting structures were simplified to Cα traces, and used to drive elastic network-driven Brownian dynamics (eBDIMS) (Orellana et al., 2016), where the protein, represented as an elastic network, was pushed from closed to open X-ray structures and vice versa to generate two initial pathways (Figure 1E1). Following side-chain reconstruction (Figure 1E2), this approach resulted in a set of 50 initial models interspersed in principal component space between the open and closed X-ray structure clusters, all with standard titration states representing neutral pH (deprotonated) that in experiments should result in a closed channel. In a duplicate set of initial models, a subset of acidic residues was modified to the most probable titration state under activating conditions (protonated), as previously described (Nury et al., 2011). Each resulting model was then subjected to unrestrained simulation in a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer for over 1 μs, producing 120 μs total sampling in the two conditions (Figure 1E3, Figure 2—figure supplement 1). Among simulations performed at each pH, we performed a feature transformation to account for the symmetry of the homopentameric protein (Figure 1E4), followed by dimensionality reduction by time-lagged independent component analysis (tICA) (Pérez-Hernández et al., 2013; Schwantes and Pande, 2013), to capture the slowest motions observed in the simulations (Figure 1E5). Further clustering in the resulting tICA space yielded kinetically meaningful states that could be used for MSM construction (Figure 1E6), validated by convergence of the main transition timescale along with other statistical tests of eigenvalues, eigenvectors and sampling (Figure 2—figure supplement 2; Figure 2—figure supplement 3).

To assess the ability of this computational approach to predict functional properties, we introduced the gain-of-function mutation I233T, located at the midpoint of the GLIC transmembrane pore, as well as the loss-of-function mutation, H235Q, a few residues away facing the intrasubunit helical bundle (Figure 1C). The midpoint position, 9' counting from the intracellular side, has been shown to constitute a hydrophobic gate that critically influences conduction properties in GLIC as in other family members (Gonzalez-Gutierrez et al., 2017), while substitution at the interfacial H235, or 11’ position, enabled crystallization of an intermediate, locally closed conformation (Prevost et al., 2012). Indeed, we confirmed by two-electrode voltage-clamp electrophysiology recordings in Xenopus oocytes that the I233T substitution enhanced proton sensitivity, shifting half-maximal activation by more than one pH unit, and producing moderately conductive channels even at neutral pH. Conversely, the histidine substitution decreased proton sensitivity, resulting in moderately conductive channels at low pH (Figure 1D). Note that maximal conduction in this plot does not necessarily indicate that all channels are open. These I233T and H235Q substitutions were introduced into additional replicate sets of initial simulation models in both deprotonated and protonated states as described above. The resulting mutants were prepared, simulated, transformed, and analyzed in the same way as the wild-type channels, achieving convergence on similar timescales (Figure 2—figure supplement 2A,D). Thus, we were able to produce six independent, statistically validated, MSMs representing distinct combinations of pH and gain- or loss-of-function mutant variations.

Free-energy landscapes capture effects of protonation and mutation

Free-energy landscapes obtained from the first eigenvector of each MSM clearly distinguished closed- and open-like regions in all six conditions (Figure 2). The tICA coordinates are presumed to represent the two slowest motions, where the largest eigenvector contributions are originating primarily from interactions between helices in the transmembrane domain (Figure 2—figure supplement 4). Along the first tIC – with the largest eigenvector contributions focused around the M2 helices (Figure 2—figure supplement 4A,C,E) – the closed X-ray structure consistently projected to the lowest free energy minimum in the landscape, suggesting a predominant population of closed channels under both deprotonated and protonated conditions. Due, in part, to its low conductance in single-channel recordings (Bocquet et al., 2007), the open probability of GLIC is not well established; however, recent cryo-electron microscopy studies support a predominance of nonconducting states even at pH below 4 (Rovšnik et al., 2021). In deprotonated conditions, the I233T substitution created a well-defined second local minimum in free energy along tIC1, centered near the open X-ray structures (Figure 2C), not seen in the wild-type or H235Q variants (Figure 2A,E). Protonated conditions further deepened this secondary minimum for both wild-type and I233T variants (Figure 2B,D), but for the H235Q variant this secondary minimum was notably displaced along tIC2 such that it did not overlap with the open X-ray structures (Figure 2F). Despite topological differences, the distribution of closed versus open X-ray structures was comparable along tIC1 in all six conditions, suggesting a common component for this slowest transition. Several so-called locally closed structures (e.g. PDB ID 3TLS) – featuring a closed-like TMD but open-like ECD – projected to a region intermediate along tIC1, but within the broad closed-state free energy well, suggesting this component primarily distinguishes TMD rather than ECD state. Interestingly, a few locally closed structures were separated out along tIC2 into a shallow secondary minimum of the closed state energy well for the deprotonated H235Q variant (Figure 2E). Conversely, a lipid-bound state (PDB ID 5J0Z) - in which the pore is expanded at the outward vestibule, although nonconductive at the hydrophobic gate - clustered with open channels. In contrast to tIC1, tIC2 did not distinguish X-ray structures in a consistent manner; the largest eigenvector contributions were interspersed between multiple transmembrane helices (Figure 2—figure supplement 4B,D,F) and higher values of this component may represent conformations sampled primarily in MD simulations.

Figure 2 with 4 supplements see all
Free-energy landscapes capture shifts upon protonation and mutation.

Free-energy landscapes projected onto the first two tICA coordinates for (A) deprotonated wild-type, (B) protonated wild-type, (C) deprotonated I233T mutant, (D) protonated I233T mutant, (E) deprotonated H235Q mutant, and (F) protonated H235Q mutant. Dots indicate experimental structures, with red representing closed X-ray structures at pH 7 (PDB ID 4NPQ), orange locally closed X-ray structures (PDB IDs 3TLS, 5MUO, 4NPP(B)), light blue modulated states (PDB ID 5J0Z) and blue open states (PDB IDs 4HFI, 3EAM, 3P4W, 4IL4, 3UU5, 4NPP(A)). Pink, magenta, and purple dots indicate closed cryo-EM structures at pH 7, 5, and 3 (PDB IDs 6ZGD, GZGJ, GZGK) respectively. Less than 100% of the channels are expected to adopt an open state even under protonated conditions (B, D, F). At protonated conditions the H235Q mutation results in an 'open'-state free-energy minimum distinct from experimental open structures (F).

Kinetic clustering distinguishes metastable open and closed states

To quantify population shifts in the evident closed- and open-like regions of the tIC landscape, we coarse-grained the MSMs into separate metastable states. By clustering according to the first dynamical MSM eigenvector, each free energy landscape could be divided into two macrostates corresponding to closed and open conformations, respectively (Figure 3A). In all but the deprotonated wild-type and H235Q landscapes, an evident free energy barrier indicated these states to be metastable (Figure 2). Based on this clustering, and consistent with the qualitative comparisons above, the fractional population of open-like macrostates moderately increased (6% to 12%) upon I233T substitution under deprotonated conditions, and increased further for both variants upon protonation (to 17% and 20% for wild-type and I233T, respectively). Upon H235Q substitution the open state population instead decreased slightly (6% to 3%) under deprotonated conditions, and reached a lower open probability (15%) compared to wild-type after protonation (Figure 3B). MSM kinetics indicated the I233T mutation accelerated opening transitions, decreasing the closed–open mean first-passage time by over one third (8.9 to 5.8 μs). Conversely, the H235Q mutation substantially slowed opening transitions by almost doubling the transition times (8.9 to 13.3 μs and 8.0 to 14.3 μs). Additionally, protonation appeared to slow closing transitions, increasing open–closed times nearly 3-fold for wild-type and I233T variants (0.5 to 1.4 μs and 0.6 to 1.7 μs for wild-type and I233T, respectively), while increasing five-fold for the H235Q mutant (0.4 to 2.1 μs) (Figure 3C).

Figure 3 with 3 supplements see all
Two-state clustering distinguishes metastable open and closed states.

(A) Two metastable states separated by the highest free energy barrier, with red representing closed-like states and blue open-like states. White dots indicate experimental structures, with labels marking closed (C), locally closed (LC) and open (O) clusters. (B) Populations of the closed (red) and open (dark blue) macrostates, with populations of hydrated conformations (>20 water molecules in the pore) marked in lighter blue. Protonated conditions consistently decreased closed-like and increased open-like populations, albeit to less than 100%. Relative to wild-type, the I233T substitution decreased closed-like populations, while H235Q decreased open-like populations, and even more so the population of hydrated conformations. Relative changes in open population reflect shifts in functional activity in Figure 1. (C) Transition rates (arrows) with numbers representing mean first passage times of crossing the highest free energy barrier. Transition rates to the closed state consistently decreased under protonated versus deprotonated conditions. In deprotonated conditions, the I233T mutation increased the rate of opening, while the H235Q mutation decreased opening rates in both conditions. (D) Hydration of the transmembrane pore. Side panels show the population of states with different hydration levels. The pore can be seen hydrating or dehydrating when crossing the main free energy barrier (white). The I233T mutation resulted in higher levels of hydration in both open and closed states, while the H235Q mutation was less hydrated in the open state.

To further validate the functional annotation of closed- and open-like macrostates, we plotted pore hydration across each free-energy landscape (Figure 3D) and extracted populations of hydrated conformations (Figure 3B). Under all conditions, the macrostate barrier corresponded to a dramatic shift in hydration levels along tIC1. All regions of tIC space sampled some dehydrated states, likely corresponding to transient, reversible, obstructions observed in individual trajectories (Figure 3—figure supplement 1; Figure 3—figure supplement 2; Figure 3—figure supplement 3). However, protonation increased hydration for all variants, particularly at larger (open-like) values of tIC1 (Figure 3D). As expected with a polar residue at the hydrophobic gate, the I233T substitution further increased hydration at all values along tIC1, suggesting more closed-like states might achieve ion conduction in this variant. Conversely, in addition to being displaced from experimental structures (Figure 2F), the open-like free-energy well was substantially less hydrated in the context of H235Q, indicating a reduced propensity for ion conduction. Thus, qualitative and quantitative comparisons of tIC landscapes supported reproducible state distinctions, and recapitulated functional effects of both protonation and mutation, supporting our model as a reasonable representation of GLIC gating.

Higher order clustering reveals conformational trends in gating

To identify conformations along the gating pathway more precisely, we reclustered each dataset according to a larger eigenvector set, obtaining models with four or five states each (Figure 4A). Note that these states are not necessarily metastable, but this clustering enables further studies of the different kinetically similar regions of the energy landscapes. Conformations corresponding to low values of both tIC1 and tIC2 (state I, Figure 4A) consistently comprised the most populated cluster (Figure 4B), with representative samples featuring visibly tilted M2 helices and a contracted pore (doi:10.5281/zenodo.5500174). Accordingly, closed as well as locally closed X-ray structures projected to the state-I cluster in all but the protonated wild-type and deprotonated H235Q systems (Figure 4A). Two more states with contracted pores clustered at low–intermediate values of tIC1, varying somewhat with system conditions; these states were distinguished along tIC2 (high for state II, low for state III; Figure 4A). Representative conformations in state II were visibly similar to state I (doi:10.5281/zenodo.5500174). For the protonated wild-type system, the closed X-ray structure projected in State II, with locally closed structures along the state I–II border. For the deprotonated H235Q variant, locally closed structures (including the experimental structure of this variant) clustered in state II, with closed states near the state I-II border (Figure 4A). State III corresponded visibly to a partly expanded pore and often high degrees of TMD subunit asymmetry (doi:10.5281/zenodo.5500174), not represented by any known X-ray structures (Figure 4A). Open and lipid-modulated X-ray structures projected to a cluster at high tIC1 and low tIC2 (state IV, Figure 4A), substantially populated only upon protonation and/or I233T substitution (Figure 4B); as expected, representative samples featured the most expanded pore relative to other states (doi:10.5281/zenodo.5500174). A final cluster at high values of both tICs (state V, Figure 4A) was the least populated in all conditions, except for the protonated H235Q variant where the secondary free energy minimum had shifted to state V rather than state IV (Figure 4B), again sampling intermediate values of pore expansion. Visual inspection of all macrostates (Figure 4A) and free-energy landscapes (Figure 2) indicated that the most likely transition pathway between closed- (states I/II) to open-like (state IV) would proceed via state III, or state II in the protonated H235Q model.

Higher order clustering of the GLIC free-energy landscapes.

(A) Each free-energy landscape can be further clustered into models with four or five macrostates that – despite not being metastable– allow for more fine-grained structural analysis of the energy landscape. White dots represent experimental structures marking closed (C), locally closed (LC) and open (O) clusters, respectively. (B) Populations for each macrostate. The different regions will be referred to as: red - State I, orange - State II, green - State III, light blue - State IV, and dark blue - State V. Conformations sampled from these states can be accessed at doi:10.5281/zenodo.5500174.

To further validate these multi-state spaces in context of past mechanistic models, we compared our clusters on the basis of conformational features implicated in channel gating.

As previously described, a key distinction between closed- and open-like states was expansion of the upper-TMD pore, quantified here by the radial distance of the upper M2 helix from the pore center-of-mass (M2 spread). On this basis, states I/II and IV/V were respectively contracted and expanded, with state III – and state V in the protonated H235Q variant – sampling intermediate values (Figure 5A; Figure 5—figure supplement 3). The expansion around the main hydrophobic gate at the 9’ position was generally expanded in the I233T mutant versus wild-type conditions (Figure 5B; Figure 5—figure supplement 4), as expected upon substitution of threonine for isoleucine. However, the H235Q mutant instead resulted in a more constricted pore with the most populated state V values approaching those of wild-type closed-state. All states were somewhat constricted at 9’ relative to the open X-ray structure (Figure 5B), though our previous measurements confirmed sampling of hydrated conformations (≥30 water molecules in the channel pore, Figure 3D) in open-like states in both conditions. Interestingly, a proposed secondary gate at the intracellular end of the pore (−2' radius, Figure 5—figure supplement 1A; Figure 5—figure supplement 9) implicated in desensitization (Gielen and Corringer, 2018) exhibited a bimodal distribution in radii, with state I largely retaining an expanded ≥5-Å radius similar to X-ray structures, while other states partially sampled a contracted radius ∼2 Å.

Figure 5 with 11 supplements see all
Probability distributions of a few variables proposed to be important in GLIC gating.

The left-most cartoons illustrate the definition of each variable, while data is presented as probability distributions with means and standard deviations plotted as bars. Colors represent the macrostates in Figure 4, and blue, red and purple horizontal lines represent the experimental structures 4HFI (Sauguet et al., 2013), 4NPQ (Sauguet et al., 2014) and 6ZGK (Rovšnik et al., 2021), respectively. The spread of the pore-lining M2 helices (A, Figure 5—figure supplement 3) is captured by the open and closed macrostates, with intermediate states taking intermediate values. However, the open-state minimum is more contracted for the protonated H235Q variant (state V). The radius from the pore center to the 9’ hydrophobic gate (B, Figure 5—figure supplement 4) captures clear differences between the three variants, where the I233T mutant has a more expanded pore than wild-type while the pore of H235Q is more contracted. The distance between the pore-lining M2 helix and the M1 helix of the neighboring subunit (C, Figure 5—figure supplement 5) is also correctly represented by the open and closed macrostates, with State III taking intermediate values. The beta expansion (D, Figure 5—figure supplement 6) yields distributions with expectation values of closed and open states aligning well with the experimental structures, while the intermediate state produces a bimodal distribution. Interestingly, the probability distributions of the closed-like states of the I233T mutation at deprotonated conditions show increased biomodality as well. The upper spread of the extracellular domain (ECD) (E, Figure 5—figure supplement 7) does not result in a clear separation of the macrostates, but a smaller pH-dependent shift can be observed. The Cα distance between E35 and T158 (FFigure 5—figure supplement 8) of primary and complementary subunits, respectively, capture large pH-dependent shifts in for all macrostates.

In parallel with expansion of the upper pore, gating transitions have been associated with contraction at TMD subunit interfaces, quantified here by the distance between proximal regions of the principal upper-M2 and complementary upper-M1 helices (M2+-M1 distance). Subunit interfaces in states I/II and IV/V were expanded and contracted respectively, with state III again sampling intermediate values (Figure 5C; Figure 5—figure supplement 5). Transitions at this interface have also been linked to relief of a helix kink, proximal to a conserved proline residue in M1. Indeed, this helix kink was more acute in states I/II, sampling even smaller angles than closed X-ray structures; conversely, the kink was largely relieved in state IV, with even larger angles than open structures (Figure 5—figure supplement 1B; Figure 5—figure supplement 10). State III again sampled intermediate values. Interestingly, kink angles for state V overlapped states I/II in deprotonated conditions, but shifted toward state IV in protonated conditions.

Gating transitions in the TMD are coupled to allosteric rearrangements in the ECD, particularly so-called β-expansion involving the first and last extracellular β-strands in each subunit. Proximal to the TMD interface, the cleft between these ECD strands is relatively expanded in closed X-ray structures but contracted in open structures, strengthening a salt bridge between β1-D32 and β10-R192 (Lev et al., 2017). Interestingly, values of β-expansion exhibited a bimodal distribution, with distinct peaks centered around closed and open structures (Figure 5D; Figure 5—figure supplement 6). In our MSMs, state-I samples were more expanded, while states IV/V were more contracted; distributions in states II/III featured distinct peaks at both closed- and open-like values.

More global gating motions in the ECD are thought to include blooming and twisting, with channel activation involving contraction and untwisting relative to the TMD. Surprisingly, no consistent state-dependent trends were noted in ECD spread (Figure 5E; Figure 5—figure supplement 7) or twist (Figure 5—figure supplement 1C; Figure 5—figure supplement 11). However, comparing for example the predominating state I in each condition revealed an overall pH dependence, with the ECD generally more contracted and twisted in protonated than deprotonated conditions (Figure 5E; Figure 5—figure supplement 1C). Time-series of ECD blooming and twisting motions from all trajectories showed rapid adaptation in response to change in pH regardless of initial seed ECD conformation (Figure 5—figure supplement 2). However, simulations started in an overall closed-like conformation displayed more flexible ECDs and sampled more broadly than simulations started in an open-like conformation, indicating state-dependent differences in ECD flexibility rather than average ECD blooming or twisting values. Interestingly, the ECD rarely contracted to the extent of open X-ray structures, nor twisted to the extent of closed X-ray structures in any condition, suggesting that crystal contacts may favor uncommon conformations in this domain. A recent cryo-EM structure of the closed state at pH 3 (Rovšnik et al., 2021) suggest a more compact ECD compared to the closed-state X-ray structure which align better with the expectation value of our data at low pH (Figure 5E), further supporting the idea that crystal contact might favor compact ECDs in the open state X-ray structure.

Finally, we studied E35, thought to represent the main proton sensor in GLIC (Nemecz et al., 2017; Hu et al., 2018). The probability distributions of the distance between E35 and T158 of the complementary subunit show clear shifts in response to changes in protonation, leading to local backbone rearrangements (Figure 5FFigure 5—figure supplement 8). Indeed, our simulations sampled even tighter contacts between these residues than observed in experimental structures. Thus, our models capture pH-dependent changes around E35, further supporting the role of E35 as an important proton sensor.

Symmetry analysis reveal protonation- and state-dependent differences

Although each GLIC molecule is composed of five identical subunits and exhibits fivefold symmetry in the context of a crystal, it remains unclear how symmetry is retained or broken in the course of channel gating. Previous simulations suggested a role for conformational asymmetry, particularly in the TMD, in facilitating closing transitions (Nury et al., 2010; Mowrey et al., 2013b). To estimate subunit symmetry, we quantified pairwise RMSDs between homologous atoms in neighboring and opposing subunits, and plotted the resulting symmetry value for each simulation frame to its corresponding position in tIC space. In this representation, regions of higher pairwise RMSD correspond to lower symmetry. To enable identification of domain-specific changes, analyses were also performed independently for the TMD and ECD in each condition (Figure 6, top bars). In both protonated datasets, pairwise RMSDs of the open state were significantly lower compared to deprotonated conditions, suggesting that protonation plays a role in the symmetrization of this state. These high symmetry levels could be deduced from both the ECD and TMD (Figure 6B, D and F, top bars). Upon channel closure, the overall symmetry decreases, which is primarily driven by the symmetry loss seen in the ECD. The TMD on the other hand seems to recover some symmetry in the closed state. At deprotonated conditions, pairwise RMSDs were generally higher across the board, and neither ECD or TMD were able to reach the low RMSD levels of the open state at protonated conditions. The ECD symmetry still seemed to decrease upon channel closure, but the differences were diminished compared to protonated conditions (Figure 6A, C and E, top bars).

State- and protonation-dependent differences in ion channel symmetry.

Heatmaps show pairwise RMSDs between all subunits of the channel, measuring the conformational symmetry of the pentamer. The two top bars show pairwise RMSDs of the transmembrane (TMD) and extracellular (ECD) domains separately, represented as stacked histograms along tIC 1. At deprotonated conditions wild-type (A) and the two mutations I233T (C) and H235Q (E) show decreased levels of symmetry in both TMD, ECD and overall. At protonated conditions wild-type (B), the I233T (D), and H235Q mutations (F) display high levels of symmetry in the open state coming from both ECD and TMD. Notably, the symmetry of the ECD is decreasing during channel closure and total closed state symmetry is coming mainly from the TMD.

Discussion

We have constructed Markov state models of a pentameric ligand-gated ion channel that enabled quantitative modeling of protonation and mutation effects, identification of intermediate states and characterization of the effect of symmetry in channel gating. Our free-energy landscapes showed deepening of the open state free energy well upon protonation, destabilization of the closed state upon I233T mutation of the hydrophobic gate (Figure 7A), and shifting of the open state toward conformations with more constricted and less hydrated pores upon H235Q mutation (Figure 7B). Our models captured shifts in free energies in agreement with our electrophysiology recordings, although only a fraction of channels would be open even under protonated conditions. In addition to capturing features of the gating mechanism already proposed for pLGICs, our MSMs allowed for further exploration of features that correlate with gating. Here, we focused on the effect of conformational symmetry between subunits and found that the open state displayed higher levels of symmetry, which was particularly enhanced after protonation (Figure 7C).

Proposed models for the free energy landscape and symmetrization in GLIC gating.

Sketches of the free-energy landscapes for protonated and deprotonated (A) wild-type and I233T variants, and (B) wild-type and H235Q variants. An open-state free energy well is formed when the channel is protonated, but only a small fraction of channels will be open at any point. The I233T mutation destabilizes the closed state (A), while the H235Q mutation results in open channels being trapped in a state with a more constricted pore behind a higher energy barrier (B). (C) Conformational symmetries of GLIC are affected by the protonation state. Upon protonation the open state displays a high level of symmetry, which is reduced to an intermediate level in the closed state. When deprotonated, the open state achieves an intermediate level of symmetry, which is reduces to low levels of symmetry in the closed state. This suggests that protonation is important for symmetrization of the open state.

Thermodynamic properties calculated from the present models were largely consistent with functional recordings, showing a shift toward relative stabilization of open versus closed states upon protonation of acidic residues or polar substitution at the 9’ hydrophobic gate (Figure 7A). Interestingly, along with a modest decrease in population, the open-like state itself changed to represent conformations with more constricted and less hydrated pores upon H235Q substitution (Figure 7B). Free energies of the wild-type gating transition have previously been estimated using string method optimization of a few collective variables assumed to be important in the gating transition (Lev et al., 2017). The variable associated with the largest barrier height in the work of Lev et al. is pore hydration, which is an integral part of the gating transition, thus suitable for comparison to the maximal energy barrier height from the MSMs. At low pH the string method hydration yielded a barrier of 1.5 kcal/mol, compared to 1.0–1.5 kcal/mol for the protonated MSM. At neutral pH the barrier from string method hydration gave 2.5 kcal/mol, while the MSM resulted in 1.5–2.0 kcal/mol height of the energy barrier. This indicated that the wild-type MSMs found transition pathways with comparable free-energy barriers.

Two-state clustering suggested a value of Popen well below 100% for GLIC, even at activating conditions. Although Popen has not been determined for GLIC, recent cryo-EM structures of closed-state GLIC solved at activating conditions point toward a significant closed population even at activating conditions (Rovšnik et al., 2021; Sauguet et al., 2014). Additionally, recent small-angle neutron scatting (SANS) experiments estimated low-pH Popen to 18% by fitting linear combinations of closed and open X-ray structures to the SANS curve, obtained from the channel population-average in solution and at room temperature (Lycksell et al., 2021). This fits well with our estimate of 17%. Other channels in the pLGIC family have also been shown to attain a range of maximal open probabilities; 10–40% for GABAARs (Pierce et al., 2019; Germann et al., 2019), 20–80% for 5-HT3s (Lambert et al., 1989; Mott et al., 2001), 0.2–3% for nAChRs (Pesti et al., 2014; Williams et al., 2012), and 90–100% for GlyRs (Ivica et al., 2021; Mangin et al., 2003), when saturated with their respective natural agonist. Given the spread in open probability between pLGIC subtypes it seems reasonable that a more distant bacterial homolog like GLIC could have a unique energy landscape. Additionally, since the open probabilities for GLIC are less than 100%, protonation could function like partial rather than full agonism.

Higher order clustering enabled more detailed investigation of the different regions of the energy landscapes, including intermediate and alternative open-like and closed-like conformations. In all cases, open and modulated crystal structures projected in state IV, while closed and locally-closed states projected into state I, or on the border between state I and state II for the protonated wild-type and deprotonated H235Q datasets. Locally closed structures, characterized by an open-like ECD and closed-like TMD, projected in the same free energy basin as closed-state 4NPQ but closer to the activation free energy barrier, indicating that the locally closed state could serve as a pre-activating state in the gating pathway. Prevost et al. solved crystal structures of multiple locally closed states trapped by various mutations, and out of these only one was capable of opening in electrophysiology experiments. Additionally, a range of conformations in the M2-M3 region were sampled in the cryo-EM structures (Prevost et al., 2012), indicating that the locally closed conformations might not represent a separate metastable state in wild-type GLIC. Further, the H235Q mutation has been shown to crystallize in a locally closed conformation at pH 4 (Fourati et al., 2018). Surprisingly, our MSMs of protonated H235Q resulted in only a modest deepening of the free energy minimum around the projected locally closed structures, but also in a heightened free-energy barrier between open and closed states, potentially facilitating a single state to be captured in experiments. This, in combination with our other observation that ECD compaction seems to be pH- rather than state-dependent, means that the most probable conformation for the H235Q variant at low pH has a closed-like TMD and more open-like ECD, similarly to the locally closed state.

Structural studies of the open, closed and locally closed crystal structures have revealed several conformational changes associated with GLIC gating (Sauguet et al., 2014). The first steps are thought to occur in the ECD through blooming and twisting motions of the entire domain. Surprisingly, our simulations did not capture state-dependent but rather pH-dependent differences in average ECD spread and twist values (Figure 5E; Figure 5—figure supplement 1; Figure 5—figure supplement 2) although larger fluctuations were observed in closed-like states, suggesting state-dependent differences in ECD flexibility rather than average spread and twist values. The recent closed-state cryo-EM structure solved at pH 3 (Rovšnik et al., 2021) displays intermediate values in ECD spread (Figure 5E), but a more twisted ECD compared to the closed X-ray structure (Figure 5—figure supplement 1C), while otherwise agreeing with values associated with the closed state TMD. While these ECD spread and twist results are somewhat counter-intuitive, the more contracted ECD could hint toward adaptation of the ECD in response to pH regardless of the TMD state. Conversely, we did observe transitions in the loops connecting the inner and outer sides of the ECD β-sandwich, which have been shown to be important for pLGIC channel function. In particular, crystal structures show breakage of the D32-R192 salt bridge in closed-state GLIC and mutational studies of D32 along with the neighboring sensor-residue E35, reveal loss of function (Bertozzi et al., 2016; Nemecz et al., 2017; Hu et al., 2018). In addition to capturing state-dependent differences in agreement with experimental data, our analysis also showed bimodality of the probability distributions at neutral pH for the closed states (Figure 5D). This effect was further enhanced upon mutation of the 9’ gate, indicating that there could be an allosteric pathway between the center of the transmembrane pore and the D32-R192 salt bridge. This is also supported by previous computational models based solely on the apparent open structure (Mowrey et al., 2013a).

Gating motions in the TMD are particularly characterized by the tilting of the pore-lining M2 helices, leading to constriction around the 9’ hydrophobic gate followed by pore dehydration. This process has previously been observed in early simulations of the TMD by Zhu and Hummer, 2010. Our models successfully captured how these M2 helix motions closely correlated with changes in 9’ radius and pore hydration levels (Figure 5—figure supplement 3Figure 5—figure supplement 4Figure 3D). Even though we capture state-dependent differences in the 9’ pore radius, and individual simulations sample pores that are as wide as the open-state structures, the most probable conformations of our open state clusters do not exhibit pores as expanded as in the majority of open crystal structures (Figure 5B). These structure have typically been co-crystallized with a hydrophobic plug of detergent molecules between the 9’ region and the top of M2, associated with a more expanded, funnel-shaped form, which is hardly sampled in simulations after plug removal. However, I233T substitution generally increased both 9’ radius to overlap better with those from the open X-ray structure, while the H235Q mutation led to overall more constriction around 9’ and less hydration, potentially impacting the ability of this variant to conduct ions. Other variables that have been associated with GLIC channel gating, including kinking of the M1 helix (Sauguet et al., 2014) and interactions in the TMD-TMD subunit interface (Fourati et al., 2018), were also captured by our MSMs (Figure 5—figure supplement 1BFigure 5C). Notably, our observations were largely consistent with (Lev et al., 2017), validating the use of these features as collective variables. Interestingly, the protonated alternative open-like state V displayed particular constriction at the −2’ gate incompatible with ion conduction (Figure 5—figure supplement 1A). This is a typical feature of desensitized states in the pLGIC family (Gielen and Corringer, 2018), although a desensitized state for GLIC has not yet been resolved. Energetically, open and desensitized states would be expected to be separated by an energy barrier, which is not the case in our models, but it is possible state V represents conformations relevant for pre-desensitization. Since the time constant for GLIC desensitization is assumed to exceed 10 s for pH values over 3.5 (Velisetty and Chakrapani, 2012) and due to the lack of structures covering this state, sampling the desensitization process without applying an enhanced seeding procedure remains a difficult task.

Analysis of the conformational symmetry of GLIC revealed a particularly symmetric open state at protonated conditions (Figure 6). At deprotonated conditions the symmetry levels were lower overall, indicating that protonation is important for channel open-state symmetrization. In all cases, significantly less symmetric ECDs could be observed in the closed state compared to the open state, although the difference was further enhanced at protonated conditions due to the more symmetric open state. The vast majority of GLIC structures cover the open state and display high levels of ECD symmetry, while the fewer structures covering the closed state are relatively poorly resolved, possibly indicating heterogeneity particularly in the ECD (Sauguet et al., 2014; Rovšnik et al., 2021). These results suggests that ECD symmetry could serve as an entropic driving force, where perturbation to the symmetric structure of the open state ECD facilitates channel closure. At activating conditions the ECD becomes protonated which facilitates high-level symmetrization characteristic for the open state. Furthermore, at protonated conditions the TMD displayed increased levels of symmetry both in the open state and closed state with less symmetric conformations in between. At deprotonated conditions, the symmetry level of the TMD transition was largely consistent with that of the protonated conditions, but the TMDs of the closed-like conformations were less symmetric. So far consensus has not been reached regarding whether homomeric channels transition symmetrically or asymmetrically, but evidence points toward asymmetric transitions being common. Although most structures derived from X-ray crystallography and cryo-EM are symmetric, partly due to the use of symmetric restraints during model building, there are a few examples of pLGICs with asymmetric structures. Recent cryo-EM structures of the 5-HT3A receptor in lipid nanodiscs revealed an asymmetric closed state, whereas a symmetrized open state was stabilized by serotonin molecules bound at all five ECD sites (Zhang et al., 2021). Additionally, a cryo-EM structure of the resting state Torpedo nAChR with toxins bound to two of the five subunits displayed asymmetries in both the ECD and TMD (Rahman et al., 2020), and cryo-EM structures of the GABAA receptor displayed ECD asymmetry in the resting state compared to the desensitized where ligands were bound (Masiulis et al., 2019). When it comes to computational work, long (15–20 μs) simulations of 5-HT3A have shown evidence of asymmetric closure of the TMD pore upon channel pre-activation (Guros et al., 2020). Mowrey et al. also proposed that asymmetric propofol binding could create unbalanced forces such that symmetry breaking would facilitate channel conformational transitions (Mowrey et al., 2013b). In our simulations, ligands (protons) were bound symmetrically across all subunits which further indicates that asymmetric transitions could be important regardless of symmetric or asymmetric ligand-binding. Understanding these asymmetries is nonetheless important for understanding the gating mechanism, but could also bring clarity regarding differences between homomeric and heteromeric pLGICs, as well as effects of different permutations of subunit assembly.

One limitation of our simulations was that protonation states were fixed throughout the simulations and selected according to consensus from previous studies. Indeed, protonation of some of titratable residues are expected to have a large effect on the results, while others are expected to be less important. Mutational studies of all titratable residues in GLIC identified that mutation of E26, D32, E35, D122, E222, H235, E243, and H277 altered GLIC activation, and out of these E35 was identified as the main proton sensor (Nemecz et al., 2017; Hu et al., 2018). In our simulations, protonation of E35 indeed had a large effect on the local backbone conformation, suggesting importance for proton-sensing (Figure 5F). Previous pKa calculations generated the following most probably protonation pattern at pH 4.6; E26, E35, E67, E75, E82, D86, D88, E177, E243, and H277 doubly protonated (Nury et al., 2011), which is also the pattern used in this paper. Other simulation studies have used similar patterns with slight modifications to model low pH conditions; E69 protonated (Nury et al., 2010), E69 protonated and three or four among five randomly chosen symmetrical copies of D/E residues protonated (Cheng et al., 2010), H127 doubly protonated (Prevost et al., 2012), E69 protonated and H127 doubly protonated (Calimet et al., 2013), E69 protonated, H127 doubly protonated, and D86 and D88 deprotonated (Lev et al., 2017), and D13, D31, D55, D91, D97, D145, D153, D154, D161, D178, D185, E14, E69, E147, E163, E181, E272, and H127 doubly protonated (Fritsch et al., 2011). Although these simulations were primarily covering the open state, we see no large discrepancies in results and can therefore expect our simulations to be generally insensitive to the permutations in protonated residues mentioned above. Once the method has matured, constant pH simulations will likely solve the combinatorial problem that protonation state selection creates (Chen et al., 2014).

In summary, our Markov states models are able to predict shifts in free energies upon change in activating stimulus (pH) as well as gain-of-function and loss-of-function mutations. The models predict relatively low values of maximal open probability, with relative differences in agreement with electrophysiology recordings. Our simulations also captured state-dependent differences in previously proposed mechanistic variables as well as state- and pH-dependent differences in channel conformational symmetry. Due to the possibility of direct comparison and validation to electrophysiology, these models are able to predict high-level statistical properties together with the details of channel conformational changes. This enables further exploration of changes that correlate with gating to further the fundamental understanding of pLGIC gating and, potentially, the development of state-selective drugs. We expect this methodology to be transferable to other channels with potential to predict properties from electrophysiology, even in absence of full equilibrium sampling.

Materials and methods

Pathway construction using eBDIMS

Request a detailed protocol

To enhance sampling, 50 seeds along the GLIC closed–open gating pathway were obtained along forward and reverse elastic-network driven Brownian dynamics (eBDIMS) simulations (Orellana et al., 2016; Orellana et al., 2019). Here, the channel was represented as an elastic network model using the Cα representation of apparent closed (PDB ID 4NPQ) (Sauguet et al., 2014) and open (PDB ID 4HFI) (Sauguet et al., 2013) structures, and driving transitions in both directions by progressively minimizing the difference in internal distances between the current and the target states. Langevin dynamics with implicit solvent and harmonic forces was used to model the dynamics, with a 12 Å intra-subunit and 8 Å inter-subunit cutoff distance, respectively. The force constant of the elastic network was 10 kcal/(mol·Å2), as previously used in Orellana et al., 2016.

Model reconstruction

Request a detailed protocol

The atomistic detail of the seeds was reconstructed using Modeller version 9.22 (Sali and Blundell, 1993) in two steps. First, side chain atoms from the template X-ray structure (PDB ID 4HFI) were added to each model, followed by a cycle of refinement with all Cα atoms restrained. Restraints on Cα atoms were then substituted with restraints on backbone hydrogen bonds, taken from helix and sheet annotations in the template PDB file, for another cycle of refinement to ensure proper secondary structure.

Molecular dynamics simulations

Request a detailed protocol

The reconstructed seed models were embedded in a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer, solvated by water and 0.1 M NaCl. Activating conditions were modeled by protonation of a subset of acidic residues (E26, E35, E67, E75, E82, D86, D88, E177, E243; H277 doubly protonated) to approximate the probable pattern at pH 4.6, as previously described (Nury et al., 2011). All systems were energy minimized with steepest descent for 10,000 steps. NPT equilibration was carried out in four steps, restraining all heavy atoms in the first cycle, then only backbone atoms, then Cα atoms, and finally the M2 helices, for a total of 76 ns. Equilibration and production runs were performed using the Amber99SB-ILDN force field with Berger lipids, together with the TIP3P water model. Temperature coupling was achieved with velocity rescaling (Bussi et al., 2007) and pressure coupling with the Parrinello-Rahman barostat (Parrinello and Rahman, 1981). The simulations were prepared and run with GROMACS versions 2018.4 and 2019.3 (Abraham et al., 2015) for 1.2 µs each, allowing for each individual trajectory to substantially deviate from the starting seed conformation and collectively sample space around the transition pathway broadly (Figure 2—figure supplement 1).

Markov state models

Request a detailed protocol

The ion channel was described with a set of 1585 features, including interatomic distances within and between subunits. Since GLIC is a symmetric homopentamer, we introduced a feature transformation where all distances originating from one subunit were scored as an independent trajectory; thus each trajectory contained information as to how one subunit moved in relation to all others. Dimensionality reduction was achieved using tIC analysis, previously shown to form a good basis for discretization of conformational space into Markov states, as the tIC eigenvetors linearly approximate the MSM eigenvectors (; Schwantes and Pande, 2013). Typically, MSMs are constructed with the assumption of exhaustive sampling of the equilibrium distribution, and inclusion of as much kinetic information as possible. Although suitable for peptide-sized systems, we found it practically unfeasible for the large-scale motions in ion channels since the MSM tends to optimize for slow but undersampled processes that may not be of primary interest. Accordingly, we discretized the state space into a few, but meaningful, tIC dimensions by omitting faster dimensions where the data was represented as a single Gaussian. Hyperparameter optimization is commonly solved by maximizing the variational approach for Markov porcesses (VAMP) (Wu and Noé, 2019) score with cross-validation to avoid overfitting; however, in our case, we found that VAMP favored exploration of slower and undersampled processes rather than convergence of a few timescales of interest, in the absence of of exhaustive equilibrium sampling. Instead, we relied on a simpler elbow approach to select appropriate hyperparameters, by calculating the open probabilities using PCCA+ (Röblitz and Weber, 2013) for various hyperparameter combinations. Based on the resulting plots (Figure 2—figure supplement 3A), we selected 300 clusters for kMeans clustering from the ’elbow’ of the plot to ensure convergence and avoid overfitting. The tIC lag times yielded consistent results within a 5–25 ns range for all datasets (Figure 2—figure supplement 3B), so a 20 ns lag time was selected, with kinetic mapping (Noé and Clementi, 2015). Instead of variationally optimizing all eigenvalues (as in VAMP), we plotted the slowest implied timescale of greatest interest using different hyperparameter combinations (Figure 2—figure supplement 3C). Indeed, our previous selection was optimal for the deprotonated conditions, so we used these parameters for all six datasets for consistency. We also note that the timescales of the selected model and the variationally optimal one were almost within the error margin for the protonated conditions. The MSM lag time was determined to 20 ns from the implied timescales, by selecting the shortest possible Markovian lag time to maximize resolution of the MSM processes (Figure 2—figure supplement 2). Two- and four-or-five-state macrostate models were obtained through PCCA+ clustering of the MSM eigenvectors (Röblitz and Weber, 2013). For the two-state model, the number of macrostates was selected in accordance to the number of observed energy minima (Figure 2) so that macrostates would represent metastable states. Higher-order clustering was done to enable more fine-grained analysis of conformations at more kinetically different regions of the phase space, and particularly to enable extraction of an intermediate cluster (Figure 4, State III). The choice of number of macrostates for these models were dependent on the structure of the MSM eigenvectors, and optimized to enable comparability between datasets and avoid macrostates with too small populations. One simulation from the protonated wild-type dataset was rejected from analysis due to obstruction of the tIC landscape with a slower but undersampled process, leaving 58.8 μs total sampling; MSMs of the other three conditions (deprotonated wild-type, deprotonated and protonated I233T, deprotonated and protonated H235Q) contained 60 μs sampling each. Sampling and convergence of each MSM were assessed through convergence of the slowest implied timescales (Figure 2—figure supplement 2A,D), Chapman-Kolmogorov tests (Prinz et al., 2011; Figure 2—figure supplement 2C,F) and by assessing the level of reversible sampling achieved (Figure 2—figure supplement 2B,E). Markov state modeling was done with PyEMMA version 2.5.7 (Scherer et al., 2015).

Electrophysiology

Request a detailed protocol

Two-electrode voltage-clamp electrophysiology was performed as previously described (Heusser et al., 2018). Briefly, GLIC cDNA subcloned in vector pMT3 was modified using commercially synthesized primers (Invitrogen, Stockholm, Sweden) and the GeneArt site-directed mutagenesis system (Thermo Fischer, Waltham, MA). Plasmids were amplified using a HiSpeed Plasmid Purfication Midi kit (Qiagen, Hilden, Germany), and verified by cycle sequencing (Eurofins Genomics GmbH, Ebersberg, Germany). Nuclei of isolated stage V–VI Xenopus laevis oocytes (EcoCyte BioScience, Dortmund, DE) were injected with 0.5–3.0 ng cDNA and stored in incubation medium (88 mM NaCl, 10 mM HEPES, 2.4 mM NaHCO3, 1 mM KCl, 0.91 mM CaCl2, 0.82 mM MgSO4, 0.33 mM Ca(NO3)2, 2 mM sodium pyruvate, 0.5 mM theophylline, 0.1 mM gentamicin, 17 mM streptomycin, 10,000 u/L penicillin, pH 8.5) for 2–7 days. Glass electrodes were pulled and filled with 3 M KCl to reach an initial resistance of 0.1–0.5 MΩ. Expressing oocytes were clamped at −70 mV using an OC-725C voltage clamp (Warner Instruments, Hamden, CT, USA), and perfused with running buffer (123 mM NaCl, 10 mM HEPES, 2 mM KCl, 2 mM MgSO4, 2 mM CaCl2, pH 8.5) at a flow rate of 0.35 mL/min. Activation buffers contained 10 mM MOPS or citrate in place of HEPES, adjusted in 0.5 pH unit increments, and were exchanged for running buffer via a VC3-8 valve controlled pressurized perfusion system (ALA Scientific Instruments, Farmingdale, NY, USA). Currents were digitized at a sampling rate of 5 kHz with an Axon CNS 1440A Digidata system using pCLAMP 10 (Molecular Devices, Sunnyvale, CA, USA). Each pH response was measured as the peak current after 1 min exposure to activation buffer, normalized to the same oocyte's response to the lowest pH tested. Each reported value represents the mean from six to eight oocytes, ± standard error of the mean. Proton concentration-dependence curves were fit by nonlinear regression with bottom and top restraints using Prism 8.0 (GraphPad Software, La Jolla, CA). Measurements from variant H235Q, collected by equivalent methods, are reproduced here from our previously published data (Fourati et al., 2018).

Data availability

Additional data including simulation parameters, Markov state models, sampled conformations and full trajectories can be accessed at https://doi.org/10.5281/zenodo.5500174.

The following data sets were generated
    1. Bergh C
    2. Heusser SA
    3. Howard R
    4. Lindahl E
    (2021) Zenodo
    Markov State Models of Proton- and Pore-Dependent Activation in a Pentameric Ligand-Gated Ion Channel.
    https://doi.org/10.5281/zenodo.5500174

References

    1. Filatov GN
    2. White MM
    (1995)
    The role of conserved leucines in the M2 domain of the acetylcholine receptor in channel gating
    Molecular Pharmacology 48:379–384.

Decision letter

  1. Toby W Allen
    Reviewing Editor; RMIT University, Australia
  2. José D Faraldo-Gómez
    Senior Editor; National Heart, Lung and Blood Institute, National Institutes of Health, United States
  3. Frédéric Poitevin
    Reviewer

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

This article presents state-of-the-art molecular dynamics simulations of the pH-gated pentameric ion channel GLIC, which has been the subject of many structural and functional studies. GLIC can be considered as a model system for pentameric ligand-gated ion channels that are responsible for fast chemical-electrical communication between cells in animals. The findings include the solution of open- and closed-like channel forms, intermediates and a "pre-desensitised" state. The approach reproduces modulation by pH and opposing mutations, correctly reproducing loss or gain of function, representing convincing proof of the success of the approach. Overall, the sampling of channel dynamics is significant and the description of state interconversions sheds new light on pLGIC mechanisms. The reviewers were convinced by the substantial changes and additions to the manuscript, as well as the new insight into the roles of pH and mutations in GLIC function. Concerns over convergence, sampling and methods descriptions have been allayed and the manuscript is now suitable for publication.

Decision letter after peer review:

Thank you for submitting your article "Markov State Models of Proton-and Gate-Dependent Activation in a Pentameric Ligand-Gated Ion Channel" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Toby Allen as Reviewing Editor and Reviewer #1, and the evaluation has been overseen by José Faraldo-Gómez as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Frédéric Poitevin (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

While the reviewers agree on the strengths and importance of the findings, the following changes or additions are required:

1. Analyse one the numerous mutants that are blocked in the so-called locally-closed state, thereby helping the community to understand the nature of this often-encountered phenotype.

2. Demonstrate convergence and sampling, including proof of lack of dependence on the initial path.

3. Better describe of the nature of the observed states (especially the locally closed state) and the tICA gating coordinates.

4. Additional methods description and clarification is needed, with clearer relation to past simulation studies to avoid misapprehensions.

Reviewer #1 (Recommendations for the authors):

The authors state that the problem with some past methods is a presumed pathway (which should be distinguished from a presumed finite set of coordinates/collective variables, as in Ref.20). Yet the MSM is not without any assumptions. The seeding of simulations in any MSM, including here, stems from an initial presumed gating transition from some simplified model. See above comments.

In the beginning of the Results section (page 5) we see a typical representation of free energy within a tICA1/2 space. See above comment. Can the key tICA1 vector be visualised, say with a vector mapping (movements of residues…) on the structure to allow it to be more readily interpreted? Moreover, it is not clear what tIC2 represents by analysis of figure 2 and later figures, and its importance to gating not immediately obvious. Is there a physical interpretation for it? Overall, I find the discussion around figure 2 to be wanting of more physical interpretation.

The authors find predominantly closed channels, although I233T (which shifts pH50 to allow conduction even at neutral pH in experiments) revealed a second open-like free energy minimum. Protonation revealed a much broader second open minimum (though still by far predominantly closed, even at low pH). This predominant closure was not seen in in past string solutions of Lev et al., and it is not clear how this is consistent with the observed pH50 for wildtype GLIC experimentally. See above comment on this, relating to the preprint in Ref.33, which appears to be the main evidence used by the authors for a low probability of a conducting form.

Regarding other states, I feel the intermediate states (such as those of states I-V with partial TMD/ECD change) have not been well characterised, or their implications explained as well as could be. Also, the authors could better visualise the location of any LC forms on the maps and explain their roles in gating (see above comments). In figure 2 the LC states are said to lump in with closed state structures (page 6). The dots in figure 4 have no colours. I guess some of the left dots are LC? See my comments above about confusion over LC. It is particularly interesting that one state with -2' constriction may relate to the desensitised state (although the authors see no barrier towards it, and so its relation to a desensitised state is not clear). Why the method cannot sample the actual desensitised state is not well explained – though may suggest a limitation in the analysis based on microsecond trajectories.

Analysis of states against past variables (as in Lev et al., Ref 20) mostly goes as expected, although the ECD behaviour does not (page 10). Though pH dependence is seen, the MSM does not reveal clear bloom and twist diffs O-C states and the authors suggest it is an artefact of x-ray crystals. While twist was less obvious, ECD spreading did differ between states in the past model based on strings by Lev et al., Is this because the finite space of Lev et al., or becoming trapped near those values? Why did the ECD spread effect disappear in this current model? Was it present after the initial seeding but vanished in subsequent MD libraries, or did this come from the initial seeding procedure?

Symmetry of the ECD decreases in closed state, and at high pH. This is as to be expected based on the higher structural diversity (and symmetry loss) of the ECD seen in Xray of the closed/high pH structure. A lot has been made about the role of asymmetry. Can the authors better explain to the reader why asymmetry in gating really matters?

In comparison to the previous string method on page 11, the authors write that barriers like 1-1.5 and 1.5-2 are lower than 1.5 and 2.5 kcal/mol, suggesting the current model extracted lower barriers – presumably being better sampled. However, given the errors I would just say they are consistent.

The statement on page 4 that Lev et al., only used short picosecond timescale simulations in their string method could be misleading for the readers. The string method used by Lev et al., makes use of large libraries of 5 and 20ps trajectories to sample the directions of change, but are carried out in parallel and repeated for many hundreds of iterations, such that the timescale sampled is much longer than that suggested by the authors here. More critically, the string method is an enhanced sampling approach such that timescales of motions cannot be judged by these values, as they adaptively sense the directions of change, allowing the approach to explore configurational space well beyond the actual simulation times. In fact, one might suggest the MSM approach used here, based on libraries of 1 microsecond free trajectories, may be limited in the states they can explore for processes that occur on much longer timescales. While it is always possible to have more sampling of random trajectories in the string method, I would say the main advantage of the present MSM is the lack of presumption of order parameters, which could be further explored in string methods. I do note, however, that one of the conclusions here (on page 12) is that the variables used by Ref.20 were consistent with the findings here.

Regarding the methods, overall the techniques used could be better explained, both for the MD/MSM expert and the general reader, avoiding jargon and relying on packaged methods, and better justifying the choices made. Things like eBDIMS, TIC etc will not mean much to many readers in the current form. The motivations for the choices and their details are important to this study.

Reviewer #2 (Recommendations for the authors):

The authors should take the opportunity of their working methodology to explain to the community the role of the LC state: is it an on-pathway form in which some of the mutants get stuck? or is it something different (off-pathway)?

– To this effect, an in-depth study of one of the many known LC mutants is needed.

– p. 2 There is a rather bold and questionable general statement in p. 2 about experimental structures, stating that "the stability required for crystal packing or cryo-EM data processing results in structures mostly representing metastable states".

I strongly suggest replacing the word "mostly" by "sometimes".

Reviewer #3 (Recommendations for the authors):

My main recommendations to the authors are very general and do not require to be answered for acceptance of publication.

Regarding the apparent subjectivity in the clustering choices and the difficulty/impossibility of exhaustive sampling for this kind of system, it would be interesting to motivate the choices made (5 states, PCCA+, …). For example, would the authors recommend their approach for simulations involving other large systems – can it be transferred – or was it thought in reaction to the particular sampling obtained here? Elements of answers to this question have the potential to help overcome the problem of robust analysis of undersampled datasets, which could benefit greatly the simulation community.

It is somewhat unclear to me what the new insights exactly are:

– The idea that protonation is mainly a driver of compaction, not of the gating is interesting. Unfortunately, there is no discussion of the combinatorial problem that protonation poses or mention of constant pH simulations that would be more realistic than applying fixed protonation that were proposed in the past – exploring the effect of a few selected protonation state changes would be interesting to see if the results are robust to a few perturbations.

– Discussion about desensitized state is interesting; it would have been nice to hear more from the authors about their take on the possibility to access this regime with simulations. In principle MSM have the potential to yield kinetic models that could be directly compared to the ones extracted from electrophysiology; what is the authors take on this? Will that ever be possible? When? Would a more experimental approach be desirable, instead of post-hoc comparison?

https://doi.org/10.7554/eLife.68369.sa1

Author response

Essential revisions:

While the reviewers agree on the strengths and importance of the findings, the following changes or additions are required:

1. Analyse one the numerous mutants that are blocked in the so-called locally-closed state, thereby helping the community to understand the nature of this often-encountered phenotype.

As we communicated earlier, this was indeed a very worthwhile effort despite the efforts required. We have completed an additional 120 µs simulations of the H235Q mutant known to crystallize in a locally closed state at low pH. We constructed two additional MSMs of this variant at deprotonated and protonated conditions and analyzed the models in a similar fashion as our four previous MSMs. We have updated Figures 1-7, Figure 2—figure supplement 1-4, Figure 3—figure supplement 3, Figure 5—figure supplement 1, and Figure 5—figure supplement 3-11 with data related to this variant. We have also extended the results and Discussion sections to include these new results. While some caveats apply (the tIC coordinates we use primarily describe the normal gating, and they are not identical for different mutants), there are also clear differences related to the local minima vs. transition barrier and hydration that we discuss in the text. Thank you for suggesting this; it turned out to be a very nice extension for the work!

2. Demonstrate convergence and sampling, including proof of lack of dependence on the initial path.

Good point. We want to be a little bit careful about absolute claims, since no method is able to truly reach convergence or be completely independent of starting conditions for systems of this size, but we have added three new figures that enable the reader to assess both sampling and how far away from the initial path we sample: Figure 2—figure supplements 1-3 together with relevant discussion in the methods section.

3. Better describe of the nature of the observed states (especially the locally closed state) and the tICA gating coordinates.

We have added more descriptions of conformations not represented by experimental structures in the results and Discussion sections. Additionally, we have added a paragraph discussing the role of the locally closed state. Regarding the tICA coordinates we added Figure 2—figure supplement 4 to clarify the meaning of the tICA coordinates.

4. Additional methods description and clarification is needed, with clearer relation to past simulation studies to avoid misapprehensions.

Very reasonable point; there are a number of groups that deserve credit in this field! We have substantially rewritten the methods section to better describe our methodological choices. We have also added three additional references to GLIC simulation studies and better clarified their similarities or differences in the introduction and Discussion sections.

Reviewer #1 (Recommendations for the authors):

The authors state that the problem with some past methods is a presumed pathway (which should be distinguished from a presumed finite set of coordinates/collective variables, as in Ref.20). Yet the MSM is not without any assumptions. The seeding of simulations in any MSM, including here, stems from an initial presumed gating transition from some simplified model. See above comments.

Very fair point; we were too brief when justifying our choice of MSMs. While they have a (neat) ability to move pathways across local maxima, and allowing multiple paths, this is of course very dependent on sampling, which we hope is now better covered with our updated figures and discussion, as well as the introduction covering more of earlier studies and alternative approaches.

In the beginning of the Results section (page 5) we see a typical representation of free energy within a tICA1/2 space. See above comment. Can the key tICA1 vector be visualised, say with a vector mapping (movements of residues…) on the structure to allow it to be more readily interpreted? Moreover, it is not clear what tIC2 represents by analysis of figure 2 and later figures, and its importance to gating not immediately obvious. Is there a physical interpretation for it? Overall, I find the discussion around figure 2 to be wanting of more physical interpretation.

We have followed the reviewer’s suggestion and projected vector mappings onto the structure in Figure 2—figure supplement 4.

The authors find predominantly closed channels, although I233T (which shifts pH50 to allow conduction even at neutral pH in experiments) revealed a second open-like free energy minimum. Protonation revealed a much broader second open minimum (though still by far predominantly closed, even at low pH). This predominant closure was not seen in in past string solutions of Lev et al., and it is not clear how this is consistent with the observed pH50 for wildtype GLIC experimentally. See above comment on this, relating to the preprint in Ref.33, which appears to be the main evidence used by the authors for a low probability of a conducting form.

Thanks for raising this important point! It is certainly interesting that the string method and MSMs gave different results related to open probabilities. Although our electrophysiology recordings are not able to reveal the fraction of open channels at maximal conduction, the shifts in open probabilities seen in our MSMs are consistent with all shifts observed in the electrophysiology recordings, also for the two additional MSMs of the loss-of-function mutation added in this revision. Additionally, these results are also in line with a new SANS study from our lab with an experimental estimation of the open probability (ref. 44) that was published after our initial submission. We have now added this reference to our discussion:

“Additionally, recent small-angle neutron scatting (SANS) experiments estimated low-pH Popen to 18% by fitting linear combinations of closed and open X-ray structures to the SANS curve, obtained from the channel population-average in solution and at room temperature [44]. This fits well with our estimate of 17%.”

Regarding other states, I feel the intermediate states (such as those of states I-V with partial TMD/ECD change) have not been well characterised, or their implications explained as well as could be.

We agree that the features of intermediate conformations could be more extensively described. First, we have described more features of state III in text, particularly on pages 6-9 and 12. Regarding the locally closed state, we have extended the results and discussion with new simulations of the H235Q mutation to better address the questions of ECD spread in relation to the locally closed state.

Also, the authors could better visualise the location of any LC forms on the maps and explain their roles in gating (see above comments). In figure 2 the LC states are said to lump in with closed state structures (page 6). The dots in figure 4 have no colours. I guess some of the left dots are LC? See my comments above about confusion over LC.

We again regret that our figures were unclear. We have therefore updated Figure 2, Figure 3 and Figure 4. We have also extended the discussion about the locally closed state.

It is particularly interesting that one state with -2' constriction may relate to the desensitised state (although the authors see no barrier towards it, and so its relation to a desensitised state is not clear). Why the method cannot sample the actual desensitised state is not well explained – though may suggest a limitation in the analysis based on microsecond trajectories.

We too are definitely excited about this, while some caution should be applied! We have added a sentence to the discussion detailing the unfeasibility of sampling the desensitization process. We also wish to clarify that the limitation lies in sampling and not analysis:

“Since the time constant for GLIC desensitization is assumed to exceed 10s for pH values over 3.5 [17] and due to the lack of structures covering this state, sampling the desensitization process without applying an enhanced seeding procedure remains a difficult task.”

Analysis of states against past variables (as in Lev et al., Ref 20) mostly goes as expected, although the ECD behaviour does not (page 10). Though pH dependence is seen, the MSM does not reveal clear bloom and twist diffs O-C states and the authors suggest it is an artefact of x-ray crystals. While twist was less obvious, ECD spreading did differ between states in the past model based on strings by Lev et al. Is this because the finite space of Lev et al., or becoming trapped near those values? Why did the ECD spread effect disappear in this current model? Was it present after the initial seeding but vanished in subsequent MD libraries, or did this come from the initial seeding procedure?

The point about ECD behavior is an important one and we regret this was not clearer in our initial submission. We also wish to clarify that we do not think differences in ECD spread and twist are artifacts in X-ray crystals: we do see differences, but not as extreme as in most structures. From our MSMs we see pH-dependent rather than state-dependent differences in ECD spread and twist values, suggesting that changes in these variables are a result of pH-sensing. To support this argument, we added Figure 5-supplement 2 with time-series of the ECD spread and twist values for all seeded wild-type simulations with starting ECD spread and twist values interspersed between those for 4HFI and 4NPQ. From this figure we see rapid adaptation (50 ns) of the ECD spread and twist values to averages distinct for each pH value. However, we do see a difference in the fluctuations of these variables based on the starting conformation, but the average is generally not affected.

Regarding comparison to the results from Lev et al., ECD spread and twist differences did not disappear from the MSMs, but represent a fast motion as the system equilibrates to the pH value of the stationary distribution. Since the free energy landscapes from our MSMs represent the stationary distribution at a certain pH value after long times, the ECD spread and twist do not show up as part of the gating mechanism in the stationary distribution, but rather as a prerequisite to reach that stationary distribution, as indicated by pH-dependent shifts in averages between protonated and deprotonated datasets. Since Lev et al., used collective variables rather than kinetic models, the fast equilibrating motions of the ECD spread might very well show up in the beginning of the gating mechanism. We extended the discussion on page 12 to include more detail about this.

Symmetry of the ECD decreases in closed state, and at high pH. This is as to be expected based on the higher structural diversity (and symmetry loss) of the ECD seen in Xray of the closed/high pH structure. A lot has been made about the role of asymmetry. Can the authors better explain to the reader why asymmetry in gating really matters?

We have followed the advice of the reviewer and added the following sentence to the discussion:

“Understanding these asymmetries is nonetheless important for understanding the gating mechanism, but could also bring clarity regarding differences between homomeric and heteromeric pLGICs, as well as effects of different permutations of subunit assembly.”

In comparison to the previous string method on page 11, the authors write that barriers like 1-1.5 and 1.5-2 are lower than 1.5 and 2.5 kcal/mol, suggesting the current model extracted lower barriers – presumably being better sampled. However, given the errors I would just say they are consistent.

We definitely agree with this (although 1kcal/mol is a 5-fold difference in kinetics!) and now refer to the free energy barriers as “comparable” to those of Lev et al., rather than “slightly lower”:

“At neutral pH the barrier from string method hydration gave 2.5 kcal/mol, while the MSM resulted in 1.5-2.0 kcal/mol height of the energy barrier. This indicated that the wild-type MSMs found transition pathways with comparable free energy barriers.”

The statement on page 4 that Lev et al., only used short picosecond timescale simulations in their string method could be misleading for the readers. The string method used by Lev et al., makes use of large libraries of 5 and 20ps trajectories to sample the directions of change, but are carried out in parallel and repeated for many hundreds of iterations, such that the timescale sampled is much longer than that suggested by the authors here. More critically, the string method is an enhanced sampling approach such that timescales of motions cannot be judged by these values, as they adaptively sense the directions of change, allowing the approach to explore configurational space well beyond the actual simulation times. In fact, one might suggest the MSM approach used here, based on libraries of 1 microsecond free trajectories, may be limited in the states they can explore for processes that occur on much longer timescales. While it is always possible to have more sampling of random trajectories in the string method, I would say the main advantage of the present MSM is the lack of presumption of order parameters, which could be further explored in string methods. I do note, however, that one of the conclusions here (on page 12) is that the variables used by Ref.20 were consistent with the findings here.

Sorry; this was not meant as a comment about the relative merits of the methods, but rather an attempt at very briefly summarize approaches. Indeed, just as MSMs enable sampling of kinetic processes far beyond the timescales of individual trajectories, so does the string method! We have removed “short ps-timescale” so that the sentence in the introduction now reads:

“This provides precious information of structural rearrangements, but the choice of collective variables in combination with the timescales of individual simulations may influence what motions are sampled.”

Regarding the methods, overall the techniques used could be better explained, both for the MD/MSM expert and the general reader, avoiding jargon and relying on packaged methods, and better justifying the choices made. Things like eBDIMS, TIC etc will not mean much to many readers in the current form. The motivations for the choices and their details are important to this study.

Check. We have added explanations to terms like tICA and eBDIMS in the beginning of the Results section and rewritten the methods section to better motivate choices.

Reviewer #2 (Recommendations for the authors):

The authors should take the opportunity of their working methodology to explain to the community the role of the LC state: is it an on-pathway form in which some of the mutants get stuck? or is it something different (off-pathway)?

– To this effect, an in-depth study of one of the many known LC mutants is needed.

In hindsight, we appreciate being pushed to do this! We have run additional simulations of the H235Q mutant. Regarding its relation to the LC state, we added a paragraph to the discussion on page 12.

– p. 2 There is a rather bold and questionable general statement in p. 2 about experimental structures, stating that "the stability required for crystal packing or cryo-EM data processing results in structures mostly representing metastable states".

I strongly suggest replacing the word "mostly" by "sometimes".

We have followed the suggestion from the reviewer so that the sentence now reads:

“since the stability required for crystal packing or cryo-EM data processing results in structures sometimes representing metastable states.”

Reviewer #3 (Recommendations for the authors):

My main recommendations to the authors are very general and do not require to be answered for acceptance of publication.

Regarding the apparent subjectivity in the clustering choices and the difficulty/impossibility of exhaustive sampling for this kind of system, it would be interesting to motivate the choices made (5 states, PCCA+, …). For example, would the authors recommend their approach for simulations involving other large systems – can it be transferred – or was it thought in reaction to the particular sampling obtained here? Elements of answers to this question have the potential to help overcome the problem of robust analysis of undersampled datasets, which could benefit greatly the simulation community.

Good idea. We have expanded the motivations of various choices in the methods section and particularly included motivations for the number of PCCA+ clusters used in our analyses. We also added a sentence with our opinion on the potential transferability of the methodology to similar systems:

“We expect this methodology to be transferable to other channels with potential to predict properties from electrophysiology, even in absence of full equilibrium sampling.”

It is somewhat unclear to me what the new insights exactly are:

– The idea that protonation is mainly a driver of compaction, not of the gating is interesting. Unfortunately, there is no discussion of the combinatorial problem that protonation poses or mention of constant pH simulations that would be more realistic than applying fixed protonation that were proposed in the past – exploring the effect of a few selected protonation state changes would be interesting to see if the results are robust to a few perturbations.

The reviewer raises an important point regarding the effect of protonation state selection on our conclusions. We agree that it would be interesting to test different combinations of protonated residues, but the computational cost of such a combinatorial study would be many times high, and it would still suffer from the limitations of fixed (integer) protonation states similar in all subunits and during the simulations. We are also working on better constant-pH simulation methodology, in particular for interacting sites, but is unfortunately far beyond the realistic scope of this paper. However, we have added a paragraph in the Discussion section (pages 13-14) comparing our protonation protocol to two experimental studies and six simulation studies. We also agree with the reviewer regarding the potential of constant pH simulations and added a sentence and additional reference (59) to the discussion:

“Once the method has matured, constant pH simulations will likely solve the combinatorial problem that protonation state selection creates [59]”

– Discussion about desensitized state is interesting; it would have been nice to hear more from the authors about their take on the possibility to access this regime with simulations. In principle MSM have the potential to yield kinetic models that could be directly compared to the ones extracted from electrophysiology; what is the authors take on this? Will that ever be possible? When? Would a more experimental approach be desirable, instead of post-hoc comparison?

The questions about directly comparing MSMs to electrophysiology are interesting and we have added a sentence at the end of the discussion giving our opinion on this:

“Due to the possibility of direct comparison and validation to electrophysiology, these models are able to predict high-level statistical properties together with the details of channel conformational changes. This enables further exploration of changes that correlate with gating to further the fundamental understanding of pLGIC gating and, potentially, the development of state-selective drugs. We expect this methodology to be transferable to other channels with potential to predict properties from electrophysiology, even in absence of full equilibrium sampling.”

https://doi.org/10.7554/eLife.68369.sa2

Article and author information

Author details

  1. Cathrine Bergh

    Science for Life Laboratory and Swedish e-Science Research Center, Department of Applied Physics, KTH Royal Institute of Technology, Solna, Sweden
    Contribution
    Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7540-5887
  2. Stephanie A Heusser

    Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University, Solna, Sweden
    Present address
    Department of Drug Design and Pharmacology, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3224-4547
  3. Rebecca Howard

    Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University, Solna, Sweden
    Contribution
    Conceptualization, Data curation, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Erik Lindahl

    1. Science for Life Laboratory and Swedish e-Science Research Center, Department of Applied Physics, KTH Royal Institute of Technology, Solna, Sweden
    2. Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University, Solna, Sweden
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    erik.lindahl@dbb.su.se
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2734-2794

Funding

Knut och Alice Wallenbergs Stiftelse

  • Erik Lindahl

Vetenskapsrådet (2017-04641)

  • Erik Lindahl

Vetenskapsrådet (2018-06479)

  • Erik Lindahl

Vetenskapsrådet (2019-02433)

  • Erik Lindahl

Swedish Research Council

  • Rebecca Howard
  • Erik Lindahl

Horizon 2020 (BioExcel (823830))

  • Erik Lindahl

Swedish National Infrastructure for Computing (2020/3-37)

  • Erik Lindahl

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by grants from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (2017–04641, 2018–06479, 2019–02433), the Swedish e-Science Research Centre, and the BioExcel Center of Excellence (EU 823830). Computational resources were provided by the Swedish National Infrastructure for Computing. The authors are grateful to Joseph Jordan, Brooke Husic, and Lucie Delemotte for helpful feedback and discussions.

Senior Editor

  1. José D Faraldo-Gómez, National Heart, Lung and Blood Institute, National Institutes of Health, United States

Reviewing Editor

  1. Toby W Allen, RMIT University, Australia

Reviewer

  1. Frédéric Poitevin

Publication history

  1. Preprint posted: March 12, 2021 (view preprint)
  2. Received: March 13, 2021
  3. Accepted: October 14, 2021
  4. Accepted Manuscript published: October 15, 2021 (version 1)
  5. Version of Record published: December 1, 2021 (version 2)

Copyright

© 2021, Bergh et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 489
    Page views
  • 87
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Genetics and Genomics
    Jeffrey V Wong et al.
    Feature Article

    Making the knowledge contained in scientific papers machine-readable and formally computable would allow researchers to take full advantage of this information by enabling integration with other knowledge sources to support data analysis and interpretation. Here we describe Biofactoid, a web-based platform that allows scientists to specify networks of interactions between genes, their products, and chemical compounds, and then translates this information into a representation suitable for computational analysis, search and discovery. We also report the results of a pilot study to encourage the wide adoption of Biofactoid by the scientific community.

    1. Computational and Systems Biology
    2. Medicine
    Abel Torres-Espín et al.
    Research Article Updated

    Background:

    Predicting neurological recovery after spinal cord injury (SCI) is challenging. Using topological data analysis, we have previously shown that mean arterial pressure (MAP) during SCI surgery predicts long-term functional recovery in rodent models, motivating the present multicenter study in patients.

    Methods:

    Intra-operative monitoring records and neurological outcome data were extracted (n = 118 patients). We built a similarity network of patients from a low-dimensional space embedded using a non-linear algorithm, Isomap, and ensured topological extraction using persistent homology metrics. Confirmatory analysis was conducted through regression methods.

    Results:

    Network analysis suggested that time outside of an optimum MAP range (hypotension or hypertension) during surgery was associated with lower likelihood of neurological recovery at hospital discharge. Logistic and LASSO (least absolute shrinkage and selection operator) regression confirmed these findings, revealing an optimal MAP range of 76–[104-117] mmHg associated with neurological recovery.

    Conclusions:

    We show that deviation from this optimal MAP range during SCI surgery predicts lower probability of neurological recovery and suggest new targets for therapeutic intervention.

    Funding:

    NIH/NINDS: R01NS088475 (ARF); R01NS122888 (ARF); UH3NS106899 (ARF); Department of Veterans Affairs: 1I01RX002245 (ARF), I01RX002787 (ARF); Wings for Life Foundation (ATE, ARF); Craig H. Neilsen Foundation (ARF); and DOD: SC150198 (MSB); SC190233 (MSB); DOE: DE-AC02-05CH11231 (DM).