Symmetry-adapted Markov state models of closing, opening, and desensitizing in α7 nicotinic acetylcholine receptors

α7 nicotinic acetylcholine receptors (nAChRs) are homopentameric ligand-gated ion channels gated by acetylcholine. These receptors play crucial roles in controlling electrical signaling within the nervous system by facilitating the passage of cations across the membrane. Recent studies have resolved and functionally annotated closed, open, and desensitized states of α7 nAChRs, providing insight into ion permeation and lipid modulation. However, the process by which α7 nAChRs transition between states remains unclear. To understand gating and lipid modulation, we generated two ensembles of molecular dynamics simulations of the apo form of α7 nAChRs, with or without cholesterol. Using symmetry-adapted Markov state modeling, we developed a five-state gating model. Free energies recapitulated functional behavior, with the closed state dominating in the absence of agonists. Notably, the transition rate from open to a non-conductive intermediate (flipped) state corresponded to experimentally measured open durations around 0.1 ms. The introduction of cholesterol relatively stabilized the desensitized state and reduced barriers between desensitized and open states. These results establish plausible asymmetric transition pathways between functionally important states, they define lipid modulation effects in the α7 nAChR conformational cycle, and provide an ensemble of structural models that could be utilized for guiding rational design strategies to develop lipidic pharmaceuticals targeting these receptors.


Introduction
Pentameric ligand-gated ion channels (pLGICs) play a crucial role in the transmission of signals within the nervous system [1][2][3].Among pLGICs, nicotinic acetylcholine receptors (nAChRs) constitute a family that can be activated by the neurotransmitter acetylcholine (ACh).These receptors are widely distributed across di erent species and can be found at both neuromuscular junctions and synaptic clefts [4].The nAChRs exhibit heterogeneity through the assembly of di erent combinations of subtypes, among which 16 have been identi ed within the human genome.These receptors demonstrate diverse physiological and pharmacological properties.Notably, the 7 subtype has the ability to form homomeric receptors and is abundant in the central nervous system, predominantly located pre-or peri-synaptically.Dysfunction of the 7 subtype can contribute to various neurological and in ammatory disorders [5][6][7][8][9].Consequently, this subtype has emerged as a prominent target for drug development.However, the lack of a comprehensive understanding of its gating mechanism has hindered the development of e ective therapeutic interventions in clinical settings [10,11].
The gating process of nicotinic acetylcholine receptors (nAChRs) primarily encompasses three functional states: closed, open, and desensitized [12].The resting, or closed, state predominates when no agonists are bound to the receptors.Upon binding of agonists in the extracellular domain (ECD), the receptors undergo a subsequent conformational change, leading to the opening of the hydrophobic pore in the transmembrane domain (TMD) and allowing cations to pass through.Following prolonged exposure to agonists, the pore closes again, entering a distinct non-conductive desensitized state.The 7 subtype also displays a distinct kinetic pro le compared to other subtypes [13].It undergoes desensitization signi cantly faster in response to acetylcholine (ACh), even faster than some of the solution exchange rates in electrophysiology recordings [13].This rapid desensitization poses challenges in achieving a consensus estimation of the open duration pro le using various techniques.
Recently, advances in cryogenic electron microscopy (cryo-EM) have enabled the resolution of all three states (closed, open, and desensitized) of the 7 subtype of nAChRs [14,15].We subsequently validated functionally annotations of these structures using molecular dynamics (MD) simulations [16].These studies have provided detailed insights into the structural basis of gating and the critical residues that govern the ion permeation pathway of 7 nAChRs.However, the precise manner in which these states are interconnected, whether symmetrically or in an asynchronous fashion, and the regulatory mechanisms underlying the dynamics of the gating process remain unclear.
Various techniques in electrophysiology have been employed since the 1970s [17] to investigate conformational changes in nicotinic acetylcholine receptors (nAChRs) [18].For instance, the phenomenon known as "nachschlag shutting"-brief closures during a burst of openings [19,20]-can be attributed to a non-conductive state preceding channel opening.It is distinct from the closed state and is referred to as " ipped" or "primed" in the literature [12].This ipped state has been functionally observed in several members within the pentameric ligand-gated ion channel (pLGIC) superfamily [21,22].Interestingly, the transition between the ipped state and the open state has been found to be independent of agonist binding.In the case of glycine receptors, another member of the pLGIC superfamily, both full agonists and partial agonists exhibit similar shut time distributions originating from the ipped state [21].A recent study on nAChRs, using kinetic modeling of an ancestral homopentamer, placed the ipped state at the center of the gating cycle [23].The authors further proposed that the agonists that agonists might function in derepression-while appearing as activation-the intrinsic pro le inherited from ancestral nAChRs [23].Despite these insights, the precise con guration of the ipped state and its relationship to other functional states remain unknown, and the transient nature of states like this can make them challenging to pursue with structural biology.
Lipids have the ability to modulate the gating of nicotinic acetylcholine receptors (nAChRs) through at least two mechanisms: altering bulk lipid properties or directly binding to the protein [24][25][26].
For 7, both these mechanisms may play a role.Structural and computational studies showed local membrane compression in the open state [14,16].We also identi ed potential binding sites for cholesterol within the 7 receptor TMD.Interestingly, PNU-120596 (PNU), a positive allosteric modulator (PAM), shared the same binding site as cholesterol at the protein-membrane interface.
Based on these ndings, we proposed a modulation mechanism for PNU in which it displaces cholesterol, thereby promoting receptor opening while preventing desensitization.These results highlight the insights provided by MD simulations in understanding the binding of lipidic modulators to various functional states and their short-range perturbation e ects.However, it should be noted that the timescale of the gating process, which occurs over milliseconds or longer, is beyond the reach of conventional MD simulations [27].
Markov state models (MSMs) may provide a useful approach for studying the kinetics of biological systems [28,29].MSMs are closely related to continuous-time kinetic modeling [30][31][32], but they represent the system as a discrete-time Markovian process, in which information can be described solely by the transition probability matrix, which de nes the probabilities of transitioning between di erent states over a discrete lag time, denoted as .The MSM framework has gained popularity as it enables the connection of multiple short MD simulations and helps bridge the timescale gap to experimental studies of kinetics [33].In an MSM, conformational changes are considered as Markovian or memoryless transitions between discrete states.To construct an MSM, one must determine an appropriate lag time and the number of microstates so that the relaxation time within each microstate is shorter than for the model to be valid.This requirement poses challenges when modeling large proteins [34,35].Modelling proteins as large as pLGICs often necessitates either a large number of microstates or a long lag time, which contradicts the advantages o ered by the MSM framework.
In this study, we utilized MSMs to investigate the gating mechanism of 7 nicotinic acetylcholine receptors after seeding simulations along three approximate gating transition pathways which encompassed experimentally resolved closed, open, and desensitized states.Taking advantage of the symmetrical nature of 7 nAChRs, we constructed MSMs using a reduced set of microstates.
By analyzing the MSMs, we were able to map the free energy landscape of the entire functional cycle of 7 nAChR, and in particular unravel the sequential order of conformational changes from the ECD to the TMD through the coupling region.The simulations allowed us to explored the role of symmetry in the gating process, and by repeating simulations of the same system in the presence of cholesterol we were able to quantitatively assess its modulating e ect on the gating of 7 nAChRs.

Initial pathway generation
Six pathways were generated using Climber [36] by conducting forward and reverse interpolations between closed (PDB ID 7KOO), open (PDB ID 7KOX), and desensitized (PDB ID 7KOQ) states.A total of 150 seeds were strategically chosen to provide a sparse sampling across the gating space.It's important to note that the initial pathways only need to adequately cover the necessary range, as we rely on unrestrained simulations for sampling and the MSMs to connect the states and modify the path.

Molecular dynamics simulations
All seeds were inserted into 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) lipid bilayers using MemProtMD [37].They were then solvated with water (TIP3P [38]) and supplemented with 0.15 M NaCl to simulate physiological conditions.In the case of the cholesterol (CHOL) system, ve cholesterol molecules were initially placed within the intersubunit binding sites for all seeds.Subsequently, seeds were embedded into bilayers composed of 33% cholesterol and 67% POPC, which mimics relevant physiological conditions.
For the simulations, approximately 400 POPC lipids were utilized for the apo system, while a combination of 300 POPC lipids and 150 cholesterol molecules were employed for the CHOL system.
The simulation boxes had dimensions of 12 nm ù 12 nm ù 18 nm.To describe the system, the CHARMM36m force eld (July 2020 version) [39] was employed.To ensure system stability, relaxation was performed in several steps.First, energy minimization was performed using the steepest descent algorithm for 5000 steps.Subsequently, a three-stage NPT relaxation was carried out, with each stage lasting 10 ns.Initially, restraints were placed on all heavy atoms, followed by only backbone atoms, and nally only on C-atoms of the protein.
The docked cholesterol molecules were also restrained during the rst two rounds of relaxation to maintain their positions, but unrestrained in the nal stage.
To maintain temperature at 300 K independently for the protein, lipid, and solvent components, the v-rescale thermostat [40] was employed.Additionally, the c-rescale barostat [41] was utilized to semi-isotropically maintain the pressure at 1 bar.Bond lengths were constrained using the LINCS [42] algorithm.
For production simulations, all restraints were removed, and systems simulated for 1-2 µs using GROMACS-2021 [43], while employing the same thermostats and barostats as during equilibration.To enhance sampling along pathways between the closed and open states in the apo system, an additional 50 simulations were performed.These simulations were conducted along the pathway to improve conformational space coverage in less sampled regions.In total, the apo system underwent 200 µs of production simulations, while the CHOL system was simulated for 277 µs.
These production simulations provided the necessary data for subsequent analyses, with data collection occurring at 1 ns intervals.

Feature extraction
For the initial model from each seed, including the starting models, a contact map was constructed by considering the interactions between all C ↵ atoms.A distance cuto of 1 nm was applied to determine if two C ↵ atoms were in contact.
Contacts with a peak-to-peak value (range) versus mean ratio of less than 0.2 were deemed insignificant and removed from the analysis.This helps lter out contacts that do not exhibit substantial variation in their interaction.After ltering, a total of 1610 signi cant contacts remained.These contacts were organized into blocks to maintain invariance to the permutation of subunits (see next section).The block arrangement ensures that the order of contacts remains the same regardless of how the subunits are arranged.For example, the rst block of contacts would encompass both intrasubunit contacts within subunit A and intersubunit contacts between subunits A and B.
To facilitate further analysis, reciprocal transformations were performed on the contact distances.

Symmetry and augmentation in feature space
The 7 nicotinic receptors are homopentameric proteins composed of ve identical subunits arranged in C5 symmetry around the central pore.This arrangement gives rise to rotational invariance, allowing us to understand the protein dynamics as a combination of ve coupled subsystems.
For example, a simpli ed two-state (C, O) model can be used to describe the conformational states of the subsystems.The model can be represented as follows: Considering the feature space, the observables that describe the dynamics of the system exhibit equivariance with respect to the symmetry group Z5 Írr 5 = eÎ.In other words, under the group actions, such as }, the Koopman propagation (K) of general observable functions (f ) commutes [44]: The simulation sampling, thus, can be augmented by applying permutations to the original feature array from simulations when they are grouped based on the system's symmetry, which e ectively improves the simulation sampling vefold.

Symmetry-adapted time lagged independent component analysis (SymTICA)
Time lagged independent component analysis (TICA) is a dimensionality reduction technique utilized used to extract the slowest dynamics present in a given feature space [45].It operates on a sequence of time series data denoted as X t = {x 1 , x 2 , ..., x N } t and computes the mean-free covariance matrix C 0 and the time-lagged covariance matrix C ⌧ given a speci c lag time .Under the assumption of reversible dynamics in TICA, the symmetry of C ⌧ is enforced numerically during the estimation process.The symmetric Koopman matrix, which encodes the kinetic information of the system, is then obtained.This matrix is further decomposed into a spectrum of eigenvalues and corresponding eigenvectors ⌫.These eigenvalues and eigenvectors provide insights into the dominant slow modes or collective motions of the system, allowing for a reduced-dimensional representation of the dynamics.
In the presence of symmetry within the dataset, the covariance matrix of the Koopman matrix exhibits a distinct block form.For instance, in the case of 5-fold symmetry, the covariance matrix (N ù N) can be represented as: Similarly, the Koopman matrix follows the same block structure, and it can be expressed as: As a result, when solving the eigenvalue problem of the Koopman matrix K⌫ = ⌫, the eigenvectors can be decomposed as To enhance data e ciency and reduce estimation error in a large kinetic model (see the toy model system in the Appendix 1), our aim is to understand the system dynamics under degeneracy.This involves seeking a symmetry-adapted mapping that remains invariant under the symmetry group.
In this case, we have , and the eigenvalue problem can be simpli ed to The features can be projected onto the dominant eigenspace (independent components, ICs) of K by direct summing over all ve subunits (sub ICs).Meanwhile, symmetry can be evaluated by analyzing the di erences between each subunit v X i , such as the standard deviation.Alternatively, symmetry-aware multidimensional scaling (MDS) [46] can be employed, where the adjacency matrix is calculated from the minimum Euclidean distances among all permutations, allowing for the estimation of symmetry.The distance matrix is then decomposed into a low-dimensional representation using classical MDS.The symmetry-aware MDS was performed using the scikit-learn package [47].

Markov state model
The SymTICA analysis was conducted using a lag time of 50 ns, and the rst two independent components (ICs) that separate the three starting states were selected.Other faster dimensions with a single Gaussian shape distribution in their sampling were excluded from further analysis [48].
After cross-validation based on the VAMP-2 score [49], 1000 microstates were clustered using the k-means algorithm [50] with k-means++ [51] initialization, implemented in the Deeptime package [52].Bayesian Markov state models [53] were estimated from the discrete microstate trajectories with a lag time of 100 ns after the Markovian property was validated using the corresponding implied timescales (ITS) [54] ranging from 50 ns to 500 ns.
Subsequently, a ve-state coarse-grained MSM was constructed using the PCCA+ algorithm [55,56], which partitions the microstates into kinetically distinct sets.The MMS was validated using the Chapman-Kolmogorov test [54].
To calculate membrane thickness, the lipophilic package [62] was employed.The thickness was determined by selecting phosphorus (P) atoms from POPC molecules and measuring the distance between the upper and lower lea ets of the membrane.
Pore hydration around residue 9 ® or 2 ® was de ned as the number of water molecules within a cylindrical region centered at the speci c residue and spanned ± 2 Å along pore axis.The tilts of the transmembrane domain helices were determined using the HELANAL algorithm [63].The tilt angle between each helix and the membrane normal was calculated, providing information about the orientation of the helices relative to the membrane.
All distributions plotted were reweighted by the stationary probability of each microstate.

Cholesterol density calculation with symmetry-corrected alignment
To evaluate the cholesterol density around the protein, conventionally, a single reference structure is used to rst align all the trajectories based on minimum root-mean-square deviation (RMSD) metric.The density is calculated with a grid-based approach.However, this approach is not suitable for asymmetric proteins from multiple trajectories because di erent trajectories may sample different asymmetric conformations but still belong to the same state.Therefore, we developed a symmetry-corrected alignment method to align the trajectories based on the symmetry of the protein.This method is an extension of the RMSD-based structure alignment algorithm implemented in MDAnalysis [61,64,65] to obtain the minimum RMSD value for each permutation and the corresponding rotation matrix.

Delineating degenerate gating dynamics with symmetry-adapted time-lagged independent component analysis
Recently validated cryo-electron microscopy (cryo-EM) structures of 7 in closed (PDB ID 7KOO), open (PDB ID 7KOX), and desensitized (PDB ID 7KOQ) states were used as initial endpoints in this study.To explore the conformational space of the protein, the Climber algorithm [36] was employed.This algorithm constructed initial pathways in all directions, resulting in a total of six pathways and a total of 150 seeds that mapped the conformational space of interest (Figure 1A, Figure S1).Each seed was embedded in a lipid bilayer, equilibrated, and further simulated without any restraints or ligands for 1-2 µs (apo system).
To address the statistical error and lack-of-sampling problem associated with exploring the entire gating space, which presumably involves asymmetric conformational changes in any or all of the ve subunits (Figure 1B), we developed symmetry-adapted time-lagged independent component analysis (SymTICA) to extract degenerate conformational changes.SymTICA is an extension of time-lagged independent component analysis (TICA) [45] that speci cally accounts for the symmetry present in the system (see Method section for detailed description).Similar to TICA, SymTICA aims to capture the slowest dynamics of the system by projecting the high-dimensional conformational space onto a low-dimensional space.Speci cally, it identi es low-dimensional features within each subsytem that exhibit symmetry-in this case, the C ↵ contacts within one subunit and with neighboring subunits.The resulting low-dimensional features, referred to as sub independent components (subICs), were then combined to construct the full independent components (ICs) by performing a direct sum operation (Figure 1C).This approach enables the extraction of collective conformational changes that respect the underlying symmetry of the system.
In our analysis, we projected all conformational snapshots obtained from our MD simulations onto the space de ned by the rst two ICs (IC1-IC2), which was constructed using 1610 C ↵ contacts (Fig-

ure S2
).This allowed us to visualize the conformational landscape and identify the dominant conformational changes associated with the gating process.IC1 primarily captured the separation between the closed state and the open/desensitized states, while IC2 primarily captured the transition between the desensitized state and closed/open states (Figure 1D-E, Figure S3).The observation that the projections along IC1 and IC2 were distinct from subIC1-IC2 indicated an asymmetric nature of the transitions (Figure 1D-E).If the transitions were symmetric, we would expect the two projections to be nearly identical.The transitions between all three states appeared to be wellconnected in the IC1-IC2 space, suggesting a continuous pathway between di erent gating states.
This provided evidence that we could proceed with constructing MSMs to quantitatively describe the dynamics of the entire gating process.

Markov state model de nes ve conformational states in the gating cycle
Following the discretization of the full ICs space into 1000 microstates (Figure S5), we performed a further coarse-graining step to identify ve metastable macrostates (Figure 2A, Figure S4).Three  This observation suggests that the SymTICA analysis and the resulting macrostate model were able to capture relevant conformational states, even when not used to seed the transitions.
The free-energy surface obtained from MSM analysis of the apo system contained distinct basins corresponding to di erent conformational states (Figure 2B).State C mapped to the deepest basin, while state O had a free energy 2.8 kcal/mol higher.Perhaps surprisingly, state D was also energetically higher than state C by 3.9 kcal/mol, as detailed in the nal section of results.State F and state I were higher in free energy by 2.9 kcal/mol and 1.9 kcal/mol, respectively, compared to state C Table 1.Thus as anticipated for the unliganded condition, the majority of channels resided in their closed state.
Kinetically, the MSM allowed us to characterize conformational transitions up to millisecond scales, which can also be illustrated as a virtual trajectory of stitched-together shorter trajectories (Supplemental video 1).The mean rst passage times (MFPTs) between di erent states shows how the system has a strong preference for the closed state C (Figure 2C, Figure S5).The 133-µs opento-ipped MFPT (Figure 2C-D) was notably comparable to the experimentally measured opening duration of 100 µs in the presence of agonists [66,67].Although opening duration time cannot be trivially obtained in the absence of agonist for 7, it is plausible that these values indicate an inherent opening rhythm of these receptor subtypes [23,68].Furthermore, MSM analysis allowed us to quantify the open burst duration, which represents the average length of spontaneous channel opening events (Figure 2E).Although there is limited electrophysiological evidence for spontaneous opening in 7-likely due to its short interval-the open burst duration obtained from MSM analysis aligned with the shortest component of the opening duration observed in previous singlechannel recordings [69].These ndings indicate that the MSM captured essential dynamics of the gating process.Consequently, the MSM was employed to investigate conformational changes in the receptor in more detail.

Sequential conformational changes in gating
We rst assessed the extent of pore hydration in each macrostate (Figure 3A).In state O, the pore was fully hydrated, consistent with our previous calculations that it allows Ca 2+ ions to permeate [16].As the channel transitioned to the desensitized state, the inner mouth of the pore became constricted.Additionally, in states I and C, the hydrophobic region at position 9 ® exhibited dewetting.Notably, state F exhibited heterogeneous pore hydration-never as hydrated as state O-and should thus be non-conductive.Sojourns in such a state could result in brief closures during a burst of opening [19,20].
Furthermore, we observed that the membrane thickness varied with the receptor state.As the system progressed from C to F the local membrane underwent compression, resulting in a subsequent decrease of approximately 4 Å in state O, which aligns with our previous ndings [16].The membrane thickness then gradually relaxed back from state O to state D and I, ultimately resembling the thickness observed in state C (Figure 3B-C).
The observed di erences in bulk membrane properties and pore pro les may be attributed to the  conformational variance among di erent states.Conformational changes in the ECD were transmitted to the pore through the coupling region between the domains (Figure 4A-C).Speci cally, the distance between the pre-M1 loop and the latch region [14] at the C-terminal end displayed a similar distribution for states C, I, and D, characterized by a loosening of contact (Figure 4 B).In contrast, both states F and O exhibited a tightening of this contact, which may account for subsequent conformational changes observed in the M1, M3, and M4 helices within the TMD (Figure 4D,F,G).
The tilt of these helices demonstrated distinct di erences between state C/D/I and state O/F.These conformational changes within the TMD could play a role in membrane compression and facilitating the opening of the pore through the tilting of the M2 helix (Figure 4E).
Interestingly, the contact between the 1-2 loop in the ECD and the M2-M3 loop in the TMD ex-

Symmetry in gating transitions
The SymTICA methodology utilizes the dynamic variances of sub independent components (sub ICs) to capture the symmetric information present in the system.The analysis of symmetry-aware multidimensional scaling (MDS) mapping and the standard deviation of sub IC1 revealed distinct conformational characteristics of the di erent states.In state C, MDS mapping and sub IC1 standard deviation analysis indicated a conserved symmetric conformation (Figure 5A, Figure S6).The majority of conformations in state O and D also exhibited symmetry, although some heterogeneity was observed.On the other hand, state F and I predominantly adopted asymmetric conformations.
To further investigate the conformational di erences, we identi ed the major local free-energy minimum (with a population > 10%) using In eCS [70] C ↵ residues provided insights into the structural variations (Figure 5B).In states C and D, only one near-symmetric conformation was observed, while state O exhibits a mixture of near-symmetric conformations.In contrast, states F and I predominantly displayed a heterogeneous distribution of asymmetric conformations, which was also re ected in the edge length distributions (Figure 5C).

Free-energy landscapes capture stabilizing e ect of cholesterol
Finally, we conducted additional computational experiments to examine the e ect of cholesterol (CHOL) by re-embedding the initial conformations into a 33% CHOL membrane system (CHOL system).Five cholesterols were docked near the predicted intersubunit cholesterol binding site [16] (Figure 6A) to ensure su cient sampling of bound states.In the desensitized state, this binding site was previously shown to be selectively accessible to cholesterol, indicating a speci c interaction between cholesterol and the receptor in this state.Conversely, in the open state, the same binding site can be occupied by the positive allosteric modulator PNU [16].A total of 277 microseconds of simulations were performed for the CHOL system, and analyzed as for the apo system, including the identi cation of ve metastable states (C, F, O, D, I) (Figure S8).
The resulting free-energy surface exhibited signi cant di erences compared to the apo system.These ndings indicate that cholesterol not only stabilizes the desensitized state but also facilitates the desensitisation transition.
CHOL demonstrated a preference for binding to its presumed intersubunit binding sites (Figure 6A,E,F in state I and state F it shifted asymmetrically into the intersubunit pocket between these helices (Figure 6E,F) from time to time.CHOL rarely accessed this binding site in states O.

Discussion
Symmetry is a pervasive phenomenon in biological systems [72].Comparable principles informed by physics have been employed in the modeling of intricate kinetics in dynamical systems, e.g. in uid dynamics and control engineering [44,73].These principles aim to extract coherent structures from high-dimensional observable time-series data by imposing speci c constraints.The incorporation of such information can signi cantly improve data-driven approaches that are susceptible to noise and data size, and work around inaccuracies.As shown in this study, explicitly accounting for symmetry in Markov State Models by introducing SymTICA, we are both able to characterize the full gating cycle of the 7 nicotinic acetylcholine receptors from molecular simulations with kinetics up to millisecond scale, but it also provides valuable insights into symmetry aspects of the transitions.
Recent advances in the decomposition of large molecular systems into independent subdomains have shown promise in unraveling the dynamics of complex systems [74][75][76].A prerequisite for such decomposition is that each subsystem should ideally be weakly coupled or completely decoupled, displaying Markovian properties on their own.In contrast, SymTICA in this study takes a di erent approach by not assuming independence among the subsystems.Instead, it focuses on symmetric systems where the global dynamics rely on asymmetric pathways between states.This approach avoids the need to describe the transitions from each asymmetric pathway, which would require an exponential number of microstates.Nonetheless, it faces challenges in constructing a global Markov state model (MSM) capable of describing all kinetics during gating.Future e orts should aim to integrate these approaches in order to enhance the construction of a more comprehensive MSM model.Deep learning frameworks have demonstrated promising capabilities in incorporating constraints, such as utilizing a matrix with a speci c manifold [77], during the training process to e ectively optimize solutions.Guided by the variational approach for Markov processes (VAMP), neural networks, such as VAMPNET [78], have been trained to maximize VAMP scores, enabling the capture of slow and nonlinear kinetics.We anticipate the principles applied employed in this study can also be applied as constraints during training in VAMPNET, thereby enhancing its performance while e ectively capturing symmetry.
We have successfully constructed an MSM that establishes connections between three functional states of the 7 nicotinic receptor: closed, open, and desensitized.Furthermore, this allowed is to identify two additional metastable states, which we suggest correspond to the ipped [12,21,22] (transitioning between closed and open) and intermediate (transitioning between closed and desensitized) functional states.Notably, the transition between the open and ipped states exhibits the characteristic nachschlag shuttings observed in single-channel recording studies [19,20], a phenomenon that is shared across members of the nicotinic receptor family.The congruence between transition times quanti ed in simulations without agonists and those observed in electrophysiology experiments conducted with varying agonist concentrations support the hypothesis that the ancestral nature of spontaneous ipping is independent of agonism [23,68].Nonetheless, because of the eeting nature of these intermediate states, it is nearly impossible to structurally depict them e.g. with cryo-EM techniques.It is also important to acknowledge that the intermediate states observed in these simulations may not necessarily fully align with the sole ipped state responsible for the electrophysiology recordings.
Our study also o ers direct structural models for asymmetric states during the gating process.We nomena of asymmetric opening are not only evident in heteropentameric channels [79,80] but also in homopentamer forms [81].This observation may be correlated with the fact that only one agonist molecule is required to initiate the opening of the 7 channel [82].
Lipids serve as important endogenous modulators for membrane proteins [24,83,84].Expanding upon our previous hypothesis [16], we investigated the role of cholesterol in the modulation of the desensitized state.Our MSM analysis supports the notion that cholesterol exhibits a differential binding pro le, e ectively stabilizing the desensitized state.Additionally, we observed that cholesterol facilitates the transition from the open to the desensitized state, potentially by introducing coupling between di erent domains in the transmembrane region.Previous studies examining the e ect of cholesterol modulation on other subtypes of nicotinic acetylcholine receptors (nAChRs) have proposed similar mechanisms, wherein the absence of cholesterol leads to an uncoupled state, rendering the channel unable to function [24].Another study highlighted the role of cholesterol in the ↵7 subtype by modifying desensitization kinetics, possibly by disrupting lipid rafts [85].Our ndings suggest that cholesterol may also modulate channel function through direct binding, further emphasizing its impact on nAChR dynamics.
Markov State Models like the ones reported here can serve as powerful frameworks to extract representative snapshots from simulations.By leveraging the Boltzmann weights, they enable the generation of structural ensembles that can be utilized in drug design studies, for example by docking novel compounds to selectively stabilize speci c states [86], which could also include stabilization of states not yet determined exxperimentally.We hope that future work can e ectively harness this pipeline to design improved modulators to enable more detailed pharmaceutical pro les for nAChRs.

Zhuang
et al. | bioR iv | December 4, 2023 | 1-33 It's important to note that within this model, there are multiple semi-degenerate asymmetric states, such as (C, C, C, C, O) and (C, C, C, O, C).
Here, C d (N_5 ù N_5) represents the diagonal block of the covariance matrix; C o1 (N_5 ù N_5) represents the o -diagonal block of the covariance matrix between features in subsystem i and i+1; C o2 (N_5 ù N_5) represents the o -diagonal block of the covariance matrix between features in subsystem i and i+2.

ZhuangFigure 1 .
Figure1.Symmetry-Adapted Time-Lagged Independent Component Analysis (SymTICA) to model nAChR gating.(A) Initial pathways generated using Climber[36] between three structural models: closed (PDB ID 7KOO), open (PDB ID 7KOX), and desensitized (PDB ID 7KOQ).(B) Asymmetric pathways during transitions between di erent states of a symmetric pentamer.These pathways exhibit the same kinetic properties.(C) Illustration of SymTICA applied to a symmetric pentamer system.Trajectory features are arranged in a block that respects the system's symmetry.Four additional pseudo-trajectories are generated by permuting the blocks.Each block is then decomposed into sub-independent components (subICs), which combine to form the full ICs through direct summation.(D) Projection of all conformational snapshots onto the subIC1-IC2 space.The three structural models are projected on the density map (grey: closed; blue: open; yellow: desensitized).(E) Projection of all conformational snapshots onto the IC1-IC2 space.The three structural models are projected on the density map (grey: closed; blue: open; yellow: desensitized).

Figure 2 .
Figure 2. MSM and kinetics of the apo system.(A) Mapping of a coarse-grained MSM with ve macrostates onto IC1-IC2 coordinates.Four structural models are projected onto the plot.(B) Free energy surface of gating projected onto IC1-IC2 coordinates.The color map represents the free energies in di erent regions.Four structural models were projected.(C) Mean rst passage times (indicated by arrows) between consecutive macrostates, with circle size indicating the relative stability of each macrostate.(D) Histogram of rst passage times between states O and F (blue) and between states O and D (orange), estimated from the MSM.(E) Histogram of open burst durations estimated from the MSM.

Figure 3 .
Figure 3. Characteristics of Macrostates: pore hydration and membrane compression.(A) Snapshots of pore hydration in each macrostate, visualizing water molecules as spheres.The inset shows a histogram of the number of water molecules around the 9 ® residue versus the 2 ® residue.This analysis was performed within a cylindrical region centered at each residue, extending ±2 Å along the pore axis.(B) Membrane thickness pro les in each macrostate were obtained by radially averaging the membrane thickness centered at the pore axis.(C) Illustration of the protein embedded within the membrane, with phosphorus (P) atoms depicted as spheres.The membrane undergoes gradual compression as the protein transitions from state C to state F and then to state O.It relaxes back from state O to state D, and further transitions through state I back to state C.

Figure 4 .
Figure 4. Sequential conformational changes of coupling region and transmembrane domain.(A) Visualization of two consecutive subunits, highlighting the contact between the pre-M1 loop and latch helix (shown in tan) and the contact between the 1-2 loop and M2-M3 loop (shown in blue).The M1, M2, M3, and M4 helices are labeled with dashed lines.(B) Histogram showing the distribution of contacts between the pre-M1 loop and latch helix for each macrostate.(C) Histogram showing the distribution of contacts between the 1-2 loop and M2-M3 loop for each macrostate.(D, E, F, G) Histograms showing the distribution of tilt angles for the M1, M2, M3, and M4 helices, respectively, in each macrostate.All distributions are weighted by the MSM weights.

Figure 5 .
Figure 5. Asymmetric intermediate states.(A) Symmetry-aware multidimensional scaling (MDS) projected sub Independent Components (IC1) onto the MDS1-MDS2 space.The zero point represents perfect symmetry with mean value, while a shifted distribution suggests sampling of asymmetric conformations.(B) Representative snapshots of the 9 ® C ↵ pentagon captured in each macrostate at each local free energy minimum.The thickness of each pentagon indicates the relative population of that conformation within the corresponding macrostate.(C) Distribution of edge lengths of the 9 ® C ↵ pentagon in each macrostate.The distribution is weighted by the MSM weights.
While state C still occupied the deepest basin , the desensitized State D was considerably stabilized in the presence of cholesterol (Figure6B) and the transition time from state O to state F remained around 100 µs (Figure6C-D).Notably, the energy barriers between states F, O, and D were attened, leading to a faster transition from state O to state D (200 µs) compared to the apo system (1.8 ms).

,Figure S7 )FFigure 6 .
Figure S7) in the desensitized state compared to other states.In state D, CHOL interacted with M253, a residue previously shown to in uence PNU binding in both experimental studies and simulations[16,71].In state C, CHOL super cially bound at the periphery of the M1 and M3 helices; observed that the closed conformation displays lower exibility compared to the open and desensitized states; higher exibility in the latter states could facilitate ion permeation.Additionally, both intermediate states exhibit predominantly asymmetric con gurations.Corresponding phe-Zhuang et al. 2023 | Markov State Models of ↵7 Nicotinic Acetylcholine Receptors bioR iv | 15 of 33

Table 1 .
MSM free energy di erences between states in apo and CHOL system.