# A Structural Model of the Human α7 Nicotinic Receptor in an Open Conformation

Letizia Chiodo, Thérèse E. Malliavin, Luca Maragliano, Grazia Cottone, Giovanni Ciccotti

► To cite this version:

Letizia Chiodo, Thérèse E. Malliavin, Luca Maragliano, Grazia Cottone, Giovanni Ciccotti. A Structural Model of the Human α7 Nicotinic Receptor in an Open Conformation. PLoS ONE, 2015, 10 (7), pp.e0133011. 10.1371/journal.pone.0133011 . pasteur-01399209

HAL Id: pasteur-01399209

https://pasteur.hal.science/pasteur-01399209v1

Submitted on 18 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# A Structural Model of the Human α7 Nicotinic Receptor in an Open Conformation

Letizia Chiodo¹, Thérèse E. Malliavin², Luca Maragliano³, Grazia Cottone⁴,⁵*, Giovanni Ciccotti⁶,⁵

1 Center for Life Nano Science @Sapienza, Istituto Italiano di Tecnologia, Rome, Italy, 2 Institut Pasteur and CNRS UMR 3528, Unité de Bioinformatique Structurale, Paris, France, 3 Center for Synaptic Neuroscience, Istituto Italiano di Tecnologia, Genoa, Italy, 4 Department of Physics and Chemistry, University of Palermo, Palermo, Italy, 5 School of Physics, University College Dublin, Dublin, Ireland, 6 Department of Physics, University of Roma “La Sapienza”, Rome, Italy

* grazia.cottone@unipa.it

## Abstract

Nicotinic acetylcholine receptors (nAchRs) are ligand-gated ion channels that regulate chemical transmission at the neuromuscular junction. Structural information is available at low resolution from open and closed forms of an eukaryotic receptor, and at high resolution from other members of the same structural family, two prokaryotic orthologs and an eukaryotic GluCl channel. Structures of human channels however are still lacking. Homology modeling and Molecular Dynamics simulations are valuable tools to predict structures of unknown proteins, however, for the case of human nAchRs, they have been unsuccessful in providing a stable open structure so far. This is due to different problems with the template structures: on one side the homology with prokaryotic species is too low, while on the other the open eukaryotic GluCl proved itself unstable in several MD studies and collapsed to a dehydrated, non-conductive conformation, even when bound to an agonist. Aim of this work is to obtain, by a mixing of state-of-the-art homology and simulation techniques, a plausible prediction of the structure (still unknown) of the open state of human α7 nAChR complexed with epibatidine, from which it is possible to start structural and functional test studies. To prevent channel closure we employ a restraint that keeps the transmembrane pore open, and obtain in this way a stable, hydrated conformation. To further validate this conformation, we run four long, unbiased simulations starting from configurations chosen at random along the restrained trajectory. The channel remains stable and hydrated over the whole runs. This allows to assess the stability of the putative open conformation over a cumulative time of 1 µs, 800 ns of which are of unbiased simulation. Mostly based on the analysis of pore hydration and size, we suggest that the obtained structure has reasonable chances to be (at least one of the possible) structures of the channel in the open conformation.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Chiodo L, Malliavin TE, Maragliano L, Cottone G, Ciccotti G (2015) A Structural Model of the Human α7 Nicotinic Receptor in an Open Conformation. PLoS ONE 10(7): e0133011. doi:10.1371/journal.pone.0133011

**Editor:** Oleg Y. Dmitriev, University of Saskatchewan, CANADA

**Received:** April 3, 2015

**Accepted:** June 22, 2015

**Published:** July 24, 2015

**Copyright:** © 2015 Chiodo et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was supported by Science Foundation of Ireland (SFI, grant 08-IN.1-1869, G. Ciccotti, http://www.sfi.ie/), Istituto Italiano di Tecnologia (IIT, seed project grant N.259SIMBEDD, G. Ciccotti, http://www.iit.it/it/home.html), and University of Palermo (FFR program 2012/2013, www.unipa.it). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
**Competing Interests:** The authors have declared that no competing interests exist.

## Introduction

The nicotinic acetylcholine receptors (nAChRs), belonging to the Cys-loop super-family of ligand-gated ion channels (or LGICs), are membrane proteins present in neurons and at neuro-muscular junctions [1]. They interact with a variety of ligands and are involved in various physiological processes [1–4] and pathologies [3, 5, 6]. The overall structure presents a fivefold symmetry, and it is composed of an extracellular domain (ECD, also called ligand binding domain, LBD, in eukaryotes), a transmembrane domain (TMD) and an intracellular part. The nAChRs comprise five subunits that can be identical or different and are arranged symmetrically to create the channel pore. Despite the large variability on the single components of the subunits, giving rise to the specificity of functions and location of these proteins, all subunits display the same architecture: the extracellular domain is composed mainly of β strands, while the transmembrane domain contains only α helices. The pore is formed by the alignment of one helix from each subunit, called M2. The orthosteric binding site of agonist and antagonist ligands is located in the LBD, and involves amino acid residues from two adjacent subunits, called principal and complementary. The C-loop, which is part of the principal subunit, caps the orthosteric binding site in the presence of a bound ligand. The nAChR channel pore, located in the TMD, opens following the binding of agonist ligands in the orthosteric site of the LBD. Conversely, the channel is mostly closed in the resting and in the inactive states, when either no ligand is present or an antagonist is bound, as well as in the desensitized, agonist-bound state.

Clearly, a detailed knowledge of the conformational transition and its dependence on ligand binding is critical to develop new pharmacological approaches influencing the receptor's biological function. Unfortunately, despite enormous efforts devoted to their crystallization, as nAChRs are membrane proteins, little high-resolution information is available on their atomic structure [7] and on the gating mechanism. Partial knowledge is obtained from low resolution experimental data, from physiological studies coupled to mutagenesis [8] and/or from structures of proteins that share some functional or structural homology.

At present, insight on nAChRs atomic structure and mechanism comes from the following Cys-loop family members structures:

i. two Electron Microscopy (EM) models of *Torpedo* acetylcholine receptor with closed [9] and open [10] pore;

ii. the X-ray crystallographic structures of two pH-gated prokaryotic homologues of nAChRs, one from *Erwinia chrysanthemi* (called ELIC) in closed form and one from *Gloeobacter violaceus* (called GLIC), initially assumed in open conformation [11–20].

iii. X-ray crystallographic structures of a water-soluble homologue of the LBD of nAChRs, the pentameric acetylcholine binding protein (AChBP), which has been co-crystallized with a number of agonists and antagonists [21–34]. The AChBP structures are the ones for which the best information on ligand-protein interactions within the orthosteric site is available for the moment [35–40]. Transferring conclusions to full-length channels, however, is not straightforward.

iv. the X-ray structures of some related eukaryotic proteins as the glutamate-gated chloride channel (GluCl) from *Caenorhabditis elegans* (open in presence of the allosteric agonist ivermectin [41] and in apo-closed conformation [42]), the human GABA<sub>A</sub> receptor [43] (open) and the mouse 5HT3 receptor [44] (whose channel conductive state is not clear from experimental data).
Each of these structures displays structural and functional features which make them not fully relevant for the homology modeling of human nAChRs. However, the model obtained from the EM imaging of *Torpedo* acetylcholine receptor has been extensively used in the past for the first attempts of nAChRs homology modeling [45–50]. Normal mode analysis has been performed on top of coarse grain models of the protein. The normal mode of lowest frequency has been exploited to force the channel in the open state, both for homopentameric [45, 47, 48] and heteropentameric models of nAChRs [50], providing the first understanding of the conformational events leading to the channel opening. More recently, the X-ray structure of the eukaryotic GluCl raised hopes to be used as template for homology models of open channels. GluCl has been already used to model the full-length α1β2γ2 GABA receptor [51], and human glycine receptors bound to agonist and antagonist [52]. Unfortunately, MD simulations [53–55] pointed out an instability of the native GluCl structure which displays a tendency to collapse once the agonist ivermectin is removed [53, 55]. What is even worse, from another study [54] it appears that the presence of ivermectin barely affects the conformation of the TMD pore, which becomes dehydrated on a time scale of one microsecond.

As for the other eukaryotic structures listed above [43, 44], they are quite recent and no MD investigation on them has been made public up to now, to the best of our knowledge.

In the present work, we provide a first step towards obtaining an all-atom model of the human α7 full-length receptor in an open conformation (Fig 1), bound to the agonist epibatidine (S1 Fig). Our approach is based first on homology modeling of the channel using different high-resolution templates, and then on extended MD simulations of the protein structure solvated and embedded in a lipid bilayer. In order to obtain a consistent full structure in an open conformation we decided to build a chimera using as templates the agonist bound AChBP structure [26] for the LBD and the open GluCl structure [41] for the TMD. However, similarly to previous studies, in an unbiased simulation we observe a collapse of this structure toward a dehydrated conformation, even in presence of epibatidine. We then employ a restraint to avoid closure, and obtain in this way a stable hydrated conformation, which we eventually refine via extended, unbiased simulations.

The stability of the structures here obtained is assessed based on a cumulative simulation time of 1.2 µs, 800 ns of which sample the unbiased open conformation. All structures were analyzed using several descriptors that have been already invoked to discriminate between open and closed channel conformations [50, 53–55], such as the pore profile, the water density inside the pore, the rigid-body motions of the TMD helices and the LBD/TMD quaternary twist. Analysis of the ion dynamics in the full-length channel, of the LBD-TMD interface and of the agonist conformations at the binding pockets contribute to the assessment of the obtained structures.

Consistent with previous findings [46, 50], our results show, in particular, that tilting of the helices and quaternary twist are not sufficient to properly describe the complex structural features of the channel. At variance, the analysis of the pore motif and hydration is much more informative: indeed, based on the features of the pore profile and the stationarity of the hydration level, we are able to discern two classes of stable structures from our simulations, and in particular to identify those representative of a possibly active state of the channel, which we refer to as putative open.

## Materials and Methods

### Homology modeling

As templates for homology modeling of the α7 nAChR we used the following X-ray structures: (i) AChBP from *Aplysia californica* bound to epibatidine (PDB entry: 2BYQ [26]), (ii) GluCl
<figure><img src="image_4.png" /><figcaption>**Fig 1. The system.** Protein (in cartoon representation) embedded into fully solvated POPC lipids slab (yellow lines) and in physiological solvent (water, blue dots); Na and Cl ions and epibatidine ligands are not shown.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g001

from *Caenorhabditis elegans* (PDB entry: 3RIF [41]). The two template structures have been manually handled to build a chimera protein in which the ECD domain (residues 0–203 in the chain A of 2BYQ) except the Cys-loop (residues 127–140 in the chain A of 2BYQ) was taken from AChBP, and the TMD domain (residues 209–339 in the chain A of 3RIF) as well as the Cys-loop (residues 130–144 in the chain A of 3RIF) were taken from GluCl. The motivation for this cutting of the structures is given in Fig 2, where the chains A of 3RIF and 2BYQ are superimposed. Indeed, in the 3RIF structure, the Cys-loop (in yellow on Fig 2C) interacts with the M2-M3 loop of TMD, whereas the corresponding region in 2BYQ (in blue on Fig 2C) moves apart from the M2-M3 loop. As the reciprocal position of the Cys-loop and M2-M3 loop is known to play an important role in the conformational transition of LGICs [53, 56–59] we decided to use the chimera protein described below as a template for the homology modeling. The homology model of α7 nAChR was then built using the chimera template and a sequence alignment of α7 and the chimera. The percentage of identity between the two sequences is 25.37%, similarly to what obtained in previous works using the same template [60, 61]. Models built from this degree of identity have a Ca RMSD of 2-3 Å with respect to the true
<figure><img src="image_5.png" /><figcaption>**Fig 2.** The chimera. A,B) X-ray crystallographic structures of the chains A and E of AChBP from *Aplysia californica* (PDB entry: 2BYQ [26]), and GluCl from *Caenorhabditis elegans* (PDB entry: 3RIF, [41]), respectively. The structures are drawn in cartoons, and the regions used for superimposition are colored in red. The residues of the superimposed structures are: (i) 55–62, 77–80, 107–118, 91–93, 138–147 for 2BYQ, and (ii) 54–61, 76–79, 110–121, 91–93, 141–150 for 3RIF. C) Superimposed subunits from 2BYQ (green) and 3RIF (in blue) using the residues defined above. The Cys-loop regions are colored in blue for 2BYQ and in yellow for 3RIF, and the cysteines of the Cys-loop are drawn in magenta for 2BYQ and in red for 3RIF. The superimposition and figure were realized using pymol [70].</figcaption></figure>

doi:10.1371/journal.pone.0133011.g002

structure [62, 63]. In these cases, a proper alignment of the sequences before building the model becomes extremely important. The software *t_coffee* 9.01 [64] was used to align the full GluCl sequence, the truncated GluCl sequence present in 3RIF, the α7 nAChR, the sequences of human α2, α3 and α4 of nAChR, and the chimera sequences, to determine the cytoplasmic region of α7. This region is thought to be folded as an α helix from the models determined using the EM map recorded on *Torpedo* acetylcholine receptor [9]. As this region has been replaced by a tripeptide in the GluCl sequence used for the determination of the 3RIF structure, we decided to replace the equivalent region of the α7 nAChR by the same tripeptide. We also truncated the N and C terminal extremities which contained many insertions in sequences alignment. The alignment of the truncated α7 nAChR sequence along the chimera sequence is described in Fig 3. Our alignment shows only 6 gaps, while functional regions such as the C-loop, the Cys-loop, as well as many conserved residues, are perfectly matched. Results are in full agreement with the sequence alignment between 2BYQ and the LBD of α7 realized by Gratzioso et al. [61] (compare Fig 1 in [61] with Fig 3). The homology modeling was realized using the script *build_model.py* from Modeller 9.10 [65], with the refinement protocol *very_slow* and 10 refined models.
<figure><img src="image_6.png" /><figcaption>**Fig 3. Sequence alignment.** Alignment, produced by t_coffee 9.01 [64], of the α7 nAChR sequence after replacing the cytoplasmic region by a tripeptide linker and of the sequence of the chimera protein used as template. The N and C terminal regions of α7 nAChR have been removed due to poor alignments (data not shown) with other human α subunits of nAChR. The UniProt sequence P36544 was used for the α7 nAChR.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g003

The best homology model was determined as the one showing the smallest percentage of residues in the generously allowed and disallowed region of the Ramachandran diagram, calculated using Procheck 3.5.4 [66]. The model was analyzed by the Molprobity server [67] molprobity.biochem.duke.edu: (i) to add the hydrogens, and (ii) to flip residues in order to optimize the sidechain interactions. The Ramachandran plot of the model output of Molprobity [68] contained 16 outliers, all located in structure loops. The model was then analyzed using the QMEAN server swissmodel.expasy.org/qmean [69], obtaining a QMEAN norm score of 0.5. This score value is satisfying as it is only slightly smaller than the corresponding score values of 0.62 and 0.53 obtained respectively on the crystallographic X-ray structures 2BYQ and 3RIF which were the initial points for the template building.

## System set-up and Molecular Dynamics simulations

Protonation state of ionizable residues at physiological pH (7.4) was predicted by using the web-based implementation of a method for pKa calculations based on continuum electrostatic model [71] (http://biophysics.cs.vt.edu/H++[72]).

MD Simulations were performed with the NAMD2.7b2 package [73, 74], with the CHARMM27 force field for proteins, ions, and water and the CHARMM36 force field for phospholipids. Topology for the epibatidine (S1 Fig) was generated by Antechamber [75]; ligand partial atomic charges were obtained by us from *ab initio* calculations with Gaussian 03 [76]. Epibatidine parameters are given as Supplementary material, see file S1 Text. The protein has been inserted into a pre-equilibrated 110 × 110 Å lipid bilayer originally composed by 356
palmitoyl-2-oleoyl-sn-glycerol-phosphatidylcholine (POPC) lipids and 5605 waters molecules, obtained from CHARMM-GUI at www.charmm-gui.org. The protein-lipid system underwent 10000 steps of conjugate gradient minimization; the temperature of the system was then increased gradually to 300 K over 2.4 ns MD simulation, by a sequence of 200 ps long simulations at constant pressure (1 atm) and temperature, each with an increase of 25 K and with random reassignment of the atomic velocities at each temperature increase, to let the membrane equilibrate around the protein. At this stage, harmonic positional restraints were applied to all protein atoms, with a force constant of 5 kcal/mol/Å². Water was then added to fully solvate the protein/lipid system; Na⁺ and Cl⁻ ions, corresponding to 100 mM solution were added to neutralize the net system charge. The total number of atoms is 142720 (26313 protein atoms, 27360 water molecules, 255 POPC lipids, 157 ions), see Fig 1. The system was again minimized for 10000 steps; harmonic positional restraints were applied to all lipid and protein atoms with a force constant of 5kcal/mol/Å². The temperature of the system was then increased gradually to 300 K over 0.6 ns MD simulation, by a sequence of 100ps long simulations at constant pressure (1 atm) and temperature, each with an increase of 50 K and with random reassignment of the atomic velocities at each temperature increase, to let the water and the ions equilibrate around the protein and the membrane. Further 5 ns were then run in the NPT ensemble at 300 K and 1 atm, by gradually removing the positional restraints on all protein atoms (force constant from 5 to 1 kcal/mol/Å², 1 ns for each force constant value), this time allowing lipids, water and ions to move freely.

Periodic boundary conditions were applied, with particle-mesh Ewald long-range electrostatics [77], using a grid spacing of 1 Å along with a sixth order B-spline charge interpolation scheme. A cutoff of 1.2 nm for Lennard-Jones potential, with a smooth switching function starting at 1.0 nm, was used. Bonds involving hydrogen atoms were constrained to their equilibrium length using the SHAKE/RATTLE algorithm [78], with time step of 1 fs.

We first ran two different simulations, using as starting condition the last conformation from the equilibration trajectory (see Table 1). In the first simulation, which we call *unrestrained*, the system is simulated for 200 ns in the NPT at 310 K and 1 atm with the protein unrestrained. We anticipate here that over the first few ps of this trajectory we observe a collapse of the pore channel toward a non-hydrated structure. The implications are thoroughly discussed in the Results section. This provided us with what we call in the following the *col-lapsed* channel, and prompted us to perform a second simulation, which we call *restrained*. The second simulation was run for 200 ns in the NPT at 310 K and 1 atm; with the aim at obtaining an open conformation, in this run a restraint has been applied on the two hydrophobic rings made of Leucine and Valine residues in position 16' and 13' that constitute the constriction point of the pore channel. We designed a “flat-bottom quadrati” potential acting on the distances between the centers of mass of Leu16' and Val13' on pairs of non-adjacent subunits [79]. The potential is a function of the distance that discourages values smaller than a threshold via a quadratic expression, and it is zero otherwise. The threshold value corresponds to the initial (open) structure; a force constant of 10 kcal/mol/Å² has been used. Also the restrained simulation was run with NAMD, the restraint being implemented via Tcl routines linked in the NAMD code itself. Time series of the restrained distances are shown in S2 Fig. The analysis of their characteristics described below, suggest that the conformations obtained in this trajectory are representative of the channel structure in an open form.

Four more simulations, in the following called *Free*₁, ..., *Free*₄ were also run in the NPT at 310 K and 1 atm, each one for 200 ns (see Table 1), starting from four initial configurations picked along the restrained trajectory at times where the restraint was not active. The similarity of the three-dimensional structures of these four initial configurations was assessed by measuring the Root Mean Square Distance of all *C*α atoms, after optimal rigid body superposition of
<table><caption>Table 1. Simulations of the α7 pentamer.</caption><thead><tr><th>Name</th><th>Starting point</th><th>Components</th><th>Duration (ns)</th></tr></thead><tbody><tr><td>unrestrained</td><td>equilibrated homology model</td><td>α7 channel, EPI</td><td>200</td></tr><tr><td>restrained</td><td>equilibrated homology model</td><td>α7 channel, EPI</td><td>200</td></tr><tr><td>Free<sub>1</sub></td><td>restrained trajectory at 20ns</td><td>α7 channel, EPI</td><td>200</td></tr><tr><td>Free<sub>2</sub></td><td>restrained trajectory at 31ns</td><td>α7 channel, EPI</td><td>200</td></tr><tr><td>Free<sub>3</sub></td><td>restrained trajectory at 34ns</td><td>α7 channel, EPI</td><td>200</td></tr><tr><td>Free<sub>4</sub></td><td>restrained trajectory at 70ns</td><td>α7 channel, EPI</td><td>200</td></tr></tbody></table>

Details of the MD trajectories. EPI: abbreviation for Epibatidine.

doi:10.1371/journal.pone.0133011.t001

pairs of structures (i.e. *Free*<sub>1</sub> vs *Free*<sub>2</sub>, i,j = 1, 4). The RMSDs values are in the range 1.3 Å - 2.2 Å, pointing out they are representative of the same conformation. No restraint has been applied along these four trajectories.

Summarizing, we have three sets of trajectories (*unrestrained*, *restrained* and *Free*<sub>1,2,3,4</sub>), which provided us with two classes of conformations: the collapsed channel (from the *unrestrained* trajectory) and the putative open channel (from the *restrained* and *Free*<sub>1,2,3,4</sub> trajectories). Therefore, the stability of the putative open state here presented has been assessed over a cumulative simulation time of about 1.0 µs, 800 ns of which are unbiased.

## Trajectory analysis

To assess the overall stability of the new homology model, we calculated: i) the $C_α$ atoms Root Mean Square Deviation (RMSD), from the configuration at the end of the equilibration stage, separately for the LBD and the TMD of the five subunits (in the following called P<sub>1</sub>, ..., P<sub>5</sub>); ii) the $C_α$ atoms Root Mean Square Fluctuations (RMSFs) with respect to the average positions calculated over the final 100 ns of each trajectory.

To describe the channel pore we calculated along all trajectories: i) the distances between specific residues of the M2 helices in non-adjacent subunits (crossed distances); ii) the time evolution of the number of water molecules in selected regions (the full pore channel, as enclosed by the five M2 helices between Gly2' and Glu20', 30 Å long; the constriction point between Leu9' and Leu16', 10 Å long, Fig 4D); iii) the pore radius, computed using the web-based implementation of a method to compute molecular channels [80] (http://bioinfo3d.cs.tau.ac.il/MolAxis/) in the configuration obtained by averaging over the last 100 ns portion of each trajectory; iv) the orientation of the five M2 helices and the quaternary twist angle between the LBD and TMD of each subunit [53].

Given that the M2 helices in the α7 do not undergo bending, we describe their orientation in terms of the tilt angle that the overall helix ($C_α$ atoms of 23 residues, from Gly2' to Glu20') forms with the protein axis. The comparison of the available X-ray structures of the putative closed and open forms of nicotinic receptors points out that in the opening/closure conformational transition the M2 helices undergo both a radial tilting and a tangential motion. For this reason, the M2 helix tilt has been usually decomposed into polar and azimuthal components [48]. As in Ref. [53], we use a reference frame centered on the center of mass of the M2 helix, with the z axis parallel to the protein symmetry axis and the x axis pointing outward along the radial direction. The polar tilt is then the angle between the z axis and the projection of the helix axis on the xz plane. The azimuthal tilt is the angle between the x axis and the projection of the helical axis on the xy plane. The twist angle is calculated as the angle between the projections of the vectors from the whole protein center of mass to the center of mass of the LBD and
<figure><img src="image_7.png" /><figcaption>**Fig 4. Pore structure and hydration.** A) Top view of M2 helices and orientation of the pore-delimiting hydrophobic side chains in the TMD of the collapsed channel (black) and open channel (red) in the structures averaged over the last 100ns of the unrestrained and restrained trajectories, respectively. The initial homology model structure is in grey. The pore size in the open channel enlarges due to side chains swinging away from the center of the pore (see e.g. in Ref. [50]). B) Snapshots of the M2 helices in the collapsed channel; C) in the putative open channel along the restrained simulation. The M2 helix of the subunit in the front is removed to reveal the water molecules inside the pore, in vdW representation. D) Hydrophobic (Leu16', Val13', Leu9') and hydrophilic (Glu20', Gly2') residues rings in the pore. The crossed distances between pairs of non-adjacent subunits for Leu at 16' and Val at 13' have been restrained with a flat-bottom potential to drive the pore towards the open state.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g004

TMD respectively (only the $C_α$ atoms are considered), on the plane perpendicular to the symmetry axis of the protein.

Quantities assessing the LBD and the LBD-TMD interface were analyzed: the C-loop/backwall distance ($d_C^{intra}$), defined as the distance between the center of mass of the C-loop ($C_α$ atoms of residues 179–188), capping the epibatidine molecules, and the center of mass of $C_α$ atoms of residues 139–140, in the backwall [40]; ligand-protein distances; intra-protein hydrogen bonds; the M2-M3 loop distance with respect to the Cys-loop and the $β1-β2$ loop, respectively, at the LBD-TMD interface. Distributions of cations and anions in the LBD region and TMD channel have been calculated as well.

## Results

In the following sections we discuss the stability of our model, and we assess its conformation via comparative inspection of pore geometry, hydration behavior, the LBD-TMD twist angle and the M2 helices orientation. The dynamics of ions within the channel is analyzed as well.
Finally, the pattern of intra-protein hydrogen bonds and the interaction of the agonist ligands with the LBD are discussed. We compare data obtained in the three sets of trajectories described in the Methods section (*restrained, unrestrained* and *Free*<sub>1,2,3,4</sub>). Results suggest that conformations extracted from these trajectories are representative of the collapsed and putative open forms of the channel.

## Model stability assessment

In this section we report mostly on the analysis of the unrestrained and restrained simulations, performed to probe the stability of the homology model and based on RMSD and RMSFs calculations (S3 Fig). The time behavior of the Cα RMSD over 200ns indicates that the conformation of the channel bound to epibatidine is stable on this time scale, both at the LBD and TMD level. The overall drift is slightly more pronounced (by 0.5 Å) in the restrained simulation than in the unrestrained one, remaining however below 4 Å in the LBD domain and below 3 Å in the TMD domain. The time behavior of the Cα RMSD over the restrained trajectory, calculated with respect to the average collapsed conformation (data not shown), is stationary around 3.3 Å for the full channel, 2.7 Å for the LBD and 1.9 Å for the TMD (averaged over the five subunits), indicating that the two conformations are globally different. Local differences are extensively discussed below.

The RMSFs along protein sequence is similar to the one obtained from previous simulations of homology models of human receptors based on different protein templates [45]. For example, in both simulations we observe a large peak in the region of helix M4 around residue 300. The high flexibility of this region is a known consequence of the structural uncertainty associated to the loop M3-M4 [81]. From a comparison of the RMSFs in the restrained and unrestrained trajectories, we clearly see no differences in peak positions, but differences in peak amplitudes. In particular, the A-loop and the F-loop are more mobile along the restrained trajectory. Otherwise, the C-loop, that caps the epibatidine molecule, and the Cys-loop, that drives the interaction with the TMD, show similar behavior in the two systems.

As for free trajectories, the time behavior of Cα RMSDs with respect the initial condition and the Cα RMSFs are both reported in Supplemental Material (S4 Fig). Results point out a substantial stability of the structures along the four free trajectories; the overall RMSD drift is less pronounced than in the unrestrained simulation (S3 Fig). Together with the observation that the four initial conditions are representative of the same putative open form (see “System set-up and Molecular Dynamics simulation”), this allow the assessment of the putative open channel at least on a cumulative time of 800 ns.

## Pore geometry and hydration

**Unrestrained trajectory.** In this simulation, a spontaneous conformational change is observed at the equilibration stage in few ps, as soon as the protein positional restraints are removed, and it is irreversible over the following 200 ns time scale of the MD trajectory. This change is characterized by a sudden decrease of the inter-residues crossed distances between M2 helices (S5 Fig, dotted lines). The drop in crossed distances is a consequence of an overall anticlockwise rigid movement of the M2 helices, which causes rotation of sidechains of hydrophobic residues Leu9', Val13', Leu16' towards the inner part of the channel (Fig 4A). Distributions of the distances between couples of residues pertaining to the same hydrophobic ring are shown in Fig 5.

The distributions were determined by merging the time series of the distances collected for the individual pairs; this to include the contributions of conformational transitions from individual subunits that would be otherwise blurred by averaging. It is evident from the multimodal
<figure><img src="image_8.png" /><figcaption>**Fig 5. Distribution of the crossed distances.** Distributions of crossed distances between pairs of non-adjacent subunits for Leu at 16', Val at 13', Leu at 9', Glu at 20' and Gly at -2'; data collected from the five subunits. The vertical dashed blue line represents the value in the starting homology modeled configuration. Black bars: unrestrained; red bars: restrained; green bars: free simulations (data collected from the four free trajectories).</figcaption></figure>

doi:10.1371/journal.pone.0133011.g005

distributions that different subunits can undergo slightly different behaviors, highlighting an intrinsic asymmetric behavior of the pentameric protein [45, 82]. The value obtained for our homology model is also reported as reference (Fig 5, blue dashed vertical line). In the collapsed channel (Fig 5, black bars), the Leu16', Val13', Leu9', Glu20' rings undergo a contraction, while the Gly2' ring delimiting the intracellular part of the pore undergoes a widening with respect to the starting structure.

This new configuration creates a hydrophobic steric hindrance that inhibits water passage, and therefore the ion permeation along the channel. Fig 6 shows the time series of the number of water molecules inside the channel pore. The number of waters in the full channel (30 Å length) is stationary around 60 over the full simulation (Fig 6, upper panel, black curve). The number of waters in the region of Leu9'-Leu16' (10 Å length) is of the order of few molecules (4 on average), with (unsuccessful) attempts of rewetting [83] at ~ 5ns, ~ 42ns, ~ 82ns, ~ 117ns, ~ 150ns, ~ 175ns and ~ 190ns along the 200 ns trajectory (Fig 6, lower panel, black curve). A representative snapshot of the M2 helices and water inside the pore in the collapsed form is shown in Fig 4B.

**Restrained trajectory.** The restraint applied on the crossed distances in the rings formed by Leu16' and Val13' residues induces a completely different hydration state in the channel. Distributions of the crossed distances are shown in Fig 5, red bars. Their time evolution, over the first 20ns of the trajectory, is shown in S5 Fig (solid lines). Remarkably, the restraint has
<figure><img src="image_9.png" /><figcaption>**Fig 6. Hydration behavior.** Upper panel: time evolution of water count in the pore lumen lined by the M2 helices. Lower panel: time evolution of water count in a region of 10 Å centered at the constriction point 13' (Val246). Black curve: unrestrained; red curve: restrained; green curve: free trajectories (data averaged over the four free trajectories).</figcaption></figure>

doi:10.1371/journal.pone.0133011.g006

non-local structural consequences: the crossed distances at the rings Glu20', Leu9', and Gly2', not affected by it, show also a structural change (Fig 5, red bars). In particular the Glu20' and Leu9' rings undergo an increase, while the Gly2' ring slightly shrinks. This suggests that the distance restraints induce a collective motion of the M2 helices. Indeed, as shown in Fig 4A, the M2 helices undergo a clockwise rotation with respect to the initial structure, with the hydrophobic sidechains pointing outside the pore.

As already observed in a study of channel opening of models of α7 based on normal mode analysis [47], the M2 rotation cause the pore to open up, allowing the water passage. In fact, as shown in Fig 6, the number of waters grows in almost 10ns, following the widening of the V-ring and L-ring due to the applied restraint. After this transition, the overall filling of the channel is almost doubled, becoming stable at around 115 water molecules (Fig 6, upper panel, red curve). Also the hydrophobic gate is now completely hydrated, with almost 40 molecules in the region corresponding to the constriction point (Fig 6, lower panel, red curve). A representative snapshot of the M2 helices and water inside the pore in this putative open form is shown in Fig 4C.

**Free trajectories.** Our most important result is that the restraint applied to the crossed distances between residues of the hydrophobic gate induced a reorganization of residues positions and, as a consequence, a new network of interactions capable of generating a conformation that is quite different from the collapsed one. To further validate and refine this conformation we submitted it to four independent, 200 ns long, free MD simulations.

In particular, the channel hydration observed in the restrained simulation is retained along the free trajectories (S6 Fig). The time evolution of the number of water molecules, averaged over the four trajectories, is shown in Fig 6 (green curve). Distributions of the crossed distances averaged over the four free trajectories show profiles similar to the one in the restrained trajectory (Fig 5, green bars).

Hydrations relaxes at stable values slightly lower than in the restrained simulation (95 waters in the full channel, 30 at the hydrophobic girdle), although still much higher than in the collapsed channel, both in the whole channel and at the hydrophobic gate. Events of sudden
<figure><img src="image_10.png" /><figcaption>**Fig 7. Hydration behavior.** Left and central panels: time evolution of water count within 3 Å from any atom of Ser2', Ser10', Thr6' and Thr12' in the pore channel. Black curves: unrestrained; red curves: restrained; green curve: trajectory Free₄. Right panel: percentages of hydrogen bond formation between water molecules and protein residues along the trajectories. Black bars: unrestrained; red bars: restrained; green bars: trajectory Free₄. The hydrogen bond is considered to be formed when the minimum distance among all possible donor/acceptor distances, between the two partners, is less than 2Å.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g007

dewetting followed by rewetting are observed in one of the four trajectories (S6 Fig, violet curve, at ~ 115 ns, ~ 135 ns, ~ 145 ns, ~ 180 ns).

The observation that water filling slowly follows channel opening after restraint application (10ns at least) suggests that the water pressure is not the primary cause for opening the pore, but instead it is a (meta-)stable conformation induced by the presence of the restraint, allowing water filling but also physiological dewetting-rewetting events. Notably, the same initial structure was used for both the restrained and unrestrained simulations, built using as template for the TMD the allegedly open GluCl structure (3RIF). After the equilibration stage, the number of water molecules in this structure is ~ 100 in the full channel and ~ 40 in the hydrophobic gate region. If water pressure was the primary factor contributing to the stabilization of the pore in an open form, then we should have not observed pore collapse in the simulation without restraints on protein atoms.

To investigate a possible role of water in stabilizing the pore from collapse successive to the effect of the restraint, we analyzed the water dynamics at the level of the hydrophilic rings in the pore, where a significant and persistent water population is expected. First, the number of water molecules found within 3 Å from any atom of Ser2', Ser10', Thr6' and Thr12' rings has been counted along the unrestrained, restrained and the Free₄ trajectory, selected as representative of the four free trajectories in which the restraints were not applied. The time behavior of the selected water count is shown in Fig 7, left and central panels. In addition, residue-water hydrogen bond analysis has been performed. The hydrogen bond is considered to be formed when the minimum distance among all possible donor/acceptor distances, between water and protein residue, is less than 2.0Å. The persistence of the selected hydrogen bond along the simulations, defined as the percentage of time along a full trajectory that the bond is formed, is also shown in Fig 7, right panel (data averaged over the five subunits).

Results show that the number of water molecules surrounding the Ser10' and Thr6' rings discriminates among the collapsed and putative open protein conformations, while similar behavior among the restrained, unrestrained and free trajectory is observed for water at Ser2' (located at the channel bottom) and at Thr12'. In particular, the number of water molecules in the Ser10' domain significantly increases when the pore is opened by the restraints (Fig 7, red curve) with respect to the collapsed trajectory (Fig 7, black curve); the water filling slowly follows channel opening (~ 10ns), as already observed in the time series of the total count of
water in the pore (see Fig 6). The filling is likely induced by the enlargement of the Leu9' ring, a consequence of the presence of the restraint (see above, "Restrained trajectory"). The number of water at Ser10' is stationary in the free trajectory, after the restraints are turned off (Fig 7, green curve). Representative snapshots of water molecules inside the pore at Ser10' and Thr 6' respectively, taken at different times along the Free₄ trajectory, are shown in S7 Fig.

Water/residue hydrogen bonds are persistently present along the restrained and free trajectories to a greater extent than in the collapsed one.

By looking at the time behavior of the number of water molecules in the pore in the various trajectories (Fig 6), we first observe a sudden dewetting at the early steps of both the restrained and unrestrained trajectories (t < 10ns). In the former however, between 10 and 20ns water is able to re-enter and stably hydrate the pore. From a careful inspection of the trajectories it emerges that this event is favored by the presence of the restraint, which creates accessible volume in the pore by pushing apart the sidechains of the hydrophobic residues. Once inside the pore, the water molecules settle in a (meta-)stable conformation (see e.g. water rings in S7 Fig) and actively contribute to keep the pore open even when the restraint is off (portions of restrained trajectory and free trajectories), by preventing the residues to change conformation and in particular the hydrophobic residues to rotate back towards the pore and consequently the helices to collapse.

## Pore profile

Results so far reported depict two different situations in the collapsed and putative open forms (both from the restrained and the free trajectories), in what concerns the pore geometry and the consequent hydration behavior. Fig 8 shows the pore radius along the channel axis (region corresponding to the M2 helices). We compare our homology model, the collapsed channel and the putative open channel obtained from the restrained simulation, together with the profiles from the structures of GluCl [41] (X-rays, open in the presence of ivermectin), GLIC [15] (X-rays, open), ELIC [17] (X-rays, closed), and Torpedo nAChR [9] (EM, closed). The pore profiles calculated in the free trajectories are shown in S8 Fig. The collapsed structure has one

<figure><img src="image_11.png" /><figcaption>**Fig 8. Pore profile.** Pore radius profiles along the pore region in the homology model configuration (blue curve), in the collapsed (black) and open (red) structures averaged over the last 100ns of the unrestrained and restrained trajectories, respectively. Also reported the profiles for open channels (GluCl and GLIC from X-rays) and closed channels (Torpedo nAChR from EM and ELIC from X-rays). Pore radius computed with MolAxis [80]. MolAxis is a channel finding algorithm, in which molecular channels are represented using "corridors", i.e. probable routes taken by a small molecule passing through a channel. MolAxis uses an algorithm that allows fast identification of corridors in the complementary space of the molecule. The black arrows, from left to right, indicate Glu20', Leu9', Thr6' and Gly2'.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g008
constriction point, corresponding to the Leu9' rings, with a minimum radius of ~ 2 Å. This constriction point at Leu9' is in common with the ELIC profile. It has been also observed in the relaxed native GluCL both in the presence and in the absence of ivermectin [54], and in a model of a Glycine receptor built on GLIC [84] (radius equal to 1.4 Å at 9'), both in the presence or absence of ethanol, included in the simulation with the aim to stabilize the Glycine receptor in an open form.

On the other hand, the radius of our collapsed channel is about 3 Å in the hydrophobic gate. The putative open channel has instead a quite large radius (~ 5 Å, comparable to the ~ 5 Å found in GLIC [15]) at the hydrophobic gate, allowing water permeation, and a constriction point of radius ~ 3 Å at the intracellular entrance (in correspondence to Gly2', near the selectivity filter (Glu1')).

For comparison, closed and putative open channel models of human α4β2 nAChR obtained by homology modeling and MD simulations [50] have radii of 2.7 Å and 3.4 Å at the extracellular entrance (corresponding to our Glu20'), respectively. A constriction point, observed for the homology model at Thr6' (Fig 8, blue curve), is present in GluCl, GLIC and ELIC, whereas it has disappeared in the putative open channel (red curve). The collapse of the homology model might be induced by the presence of this constriction point, whereas the action of the restraints induces the opening of the Thr6' ring and the consequent water filling already observed in Fig 7.

## Large scale motions: quaternary twist and M2 orientation

Normal mode analysis on full-length models of nicotinic receptors pointed out that the lowest-frequency mode involves a global twisting motion of the LBD domain relative to the TMD domain, in which the domains rotate in opposite direction around the pore axis [47, 49, 85-87]. This conformational change has been related to the channel gating mechanism [47, 49]. The value of the twist angle gives a measure of gate locking when the channel is bound to antagonists or in the absence of ligands, and it is expected to be larger for a closed channel than for an open one. The distributions of the twist angle values in our trajectories are shown in Fig 9, left panel. Data are averaged over the five subunits. In the inset, the distributions of data collected from the five subunits are shown; this to highlight the contribution of individual subunits that would be otherwise blurred by simple averaging.

In agreement with the literature [53, 84], we observe in all trajectories a large change of the twist angle with respect to the homology model: the twist angle averaged over the trajectories relaxes to values higher than the one observed in the starting structure (Fig 9, left panel), but the protein undergoes a smaller twist along the restrained trajectory with respect to the unrestrained one, with average values of 21° and 24°, respectively.

Over the four free trajectories (distributions not shown), the twist angle is lower than in the restrained simulation in two cases (peak at ~ 20°), similar to the restrained value in one case (~ 21°) and close to the unrestrained value in the fourth case (~ 24°). On average, as shown in Fig 9 (left panel, green bars), the twist angle is decreasing from the collapsed to the open conformations sampled both in the restrained and the free trajectories.

For sake of comparison, the value for simulated native GLIC and GluCl (in the presence of ivermectin) is about 15°, while in GluCl simulated in the absence of ivermectin the twist angle increase to 20°[53]. The value for the simulated closed ELIC structure is 23°. In a model of glycine receptor built on GLIC [84] and bound with ethanol to stabilize the open form, the twist angle is ~ 14°, slightly smaller than the value of ~ 16° in the form obtained in the absence of ethanol. To summarize, the variation of twist observed in the present work agrees with the observations previously made in MD simulations: the twist angle increases with the tendency of the channel to close.
<figure><img src="image_12.png" /><figcaption>**Fig 9. Twist and tilt angles.** Left panel: distributions of quaternary twist angle; central panel: distributions of the M2 helix polar tilt angle; right panel: distributions of M2 helix azimuthal tilt angle. Values averaged over the five subunits. In the insets, distributions calculated collecting data from the five subunits. The vertical dashed blue line represents the value in the starting homology modeled configuration. Black bars: unrestrained; red bars: restrained; green bars: free simulations (data averaged over the four free trajectories).</figcaption></figure>

doi:10.1371/journal.pone.0133011.g009

The orientation behavior of the M2 helices is more complex. The distributions of polar and azimuthal tilt angle components, calculated by averaging over the five subunits, are shown in Fig 9, central and right main panels, respectively. The distributions of data collected from the five subunits are also shown in the inset of each panel, highlighting the asymmetrical, non-certed motion of the helices [10, 47, 48, 50, 88]. For the polar tilt, the central main panel shows a change from slightly negative values in the collapsed channel (black bars, peak at -2°) to slightly positive values (red bars, peak at 1.1°) in the putative open conformation obtained from the restrained simulation. This variation is in agreement with the channel pore shape going from a straight to a V-shaped form (Fig 4B and 4C). The values averaged over the four independent free trajectories have a distribution with higher weight for positive polar tilt angles. For sake of comparison, the polar tilt angle increases from ~ 7° in native GluCl simulated in the absence of ivermectin to ~ 9° in GluCl stabilized in the open form with ivermectin [53]. In simulations of prokaryotes [53], the polar tilt angle was ~ 1° in the closed channel of ELIC and ~ 7° in the GLIC open channel. As for the azimuthal tilt, in the restrained trajectory the distribution averaged over the five subunits (right main panel, red bars) shifts toward positive values with respect to the starting homology model, with a large variability among the five subunits (see inset). This high variability is even more evident in the four free trajectories. Typical values for the azimuthal angle found in other simulations [53] are -8.5° in closed ELIC and ~ 2° in open GLIC; ~ 3° and ~ 5° in native GluCl in the absence and in the presence of ivemectin respectively. Overall, the slight increase of polar tilt angle from the closed to open conformations is in agreement with the observations previously made in the literature in X-ray structures and MD trajectories, although no specific values can be deciphered up to now.

## Ions dynamics

The α7 nAChR is cation selective. To further assess the α7 model proposed in this study, we performed an analysis of the distribution of both cations and anions in the channel pore. Fig 10 shows the distribution of the z-coordinate of ions along the restrained and unrestrained
<figure><img src="image_13.png" /><figcaption>**Fig 10. Ions distribution.** Blue bars: cations, that can permeate the nicotinic channel. Red bars: anions. Top panels show the channel pore profiles along the full-length protein. Central panel: restrained; lower panel: unrestrained trajectory. The cyan region including the cyan rectangle in the top panel) highlights the channel between Glu20' and Leu9'. The vertical black lines on the pore profile highlight the cytoplasmic and the extracellular pore limits, from Gly2' to Glu20'.</figcaption></figure>

doi:10.1371/journal.pone.0133011.g010

trajectories, from the intracellular to the extracellular exits, both for the sodium (cations) and chlorine (anions) ions. Ion populations calculated along the free trajectories are shown in S9 Fig. The overall profile along the full-length channel is in agreement with results reported in literature for nicotinic receptors, though obtained with higher ions concentration and/or via the application of a transmembrane potential along Brownian Dynamics simulations (see e.g Fig 7 in Ref. [89]). Starting from the bottom of the protein (-30 Å), the anions populate more the cytoplasmic extremity of the channel (-30 Å<z<-20 Å), but the cations could enter more into the channel, in the portion comprised between -20Å and 0 Å; this in particular in the unrestrained trajectory and, on average, along the four free trajectories. The hydrophobic gate, in the range 0 Å<z<20 Å is not populated. Indeed, given the low, physiological salt concentration (100 mM) used here to neutralize the total system charge, and, most importantly, the time scale of our MD simulation (200ns for each trajectory), it is highly improbable to detect any conductive event crossing the hydrophobic gate, where the free energy barrier to ion permeation is located. However, despite the low concentration here used with respect to other studies devoted to ion motion [89, 90], along one of the free trajectories (Free₄) we observed the simultaneous presence of two cations at the same time in the region just below the Leu9'
hydrophobic ring (-10 Å<z<0Å) [91]. This is a precursor event for the ion passage through the hydrophobic gate [89]. Ions population in the LBD shows an alternation of anions and cations. We observe in particular a cation excess located at z ~ 30Å and z ~ 50Å, in correspondence to charged residue rings composed by Asp38, Glu39 and Asp91, Glu92 respectively. These rings corresponds to the ones identified in Ref. [90] to be associated with prolonged cation dwell time.

## LBD-TMD interface

So far, we focused our attention mostly on the TMD portion of the protein. The above reported analyses of the pore structure, the hydration behavior, the M2 orientation and the ion dynamics point out that the structures we obtain, both in the restrained and free simulations, are representative of a stationary, active state of the channel. It is likewise important to assess how the remainder of the homology model we propose behave along the simulated trajectories.

The analysis of the nicotinic receptors gating mechanism pointed out that the TMD/LBD interface region is involved in the transmission of the signal from the LBD towards the channel region, through a cascade mechanism [8, 20, 46, 48]. Actually, the coupling of the agonist binding to channel gating has been explained on the basis of the concerted interaction of the M2-M3 loop in the TMD with the Cys-loop and the β1-β2 loop in the LBD [56, 92] (Fig 11).

In particular, MD simulations have shown how the sliding motion of the M2-M3 loop from the β1-β2 loop towards the Cys-loop could be responsible for the initiation of pore widening in the closed-to-open transition in response to ligand binding [46, 48]. In the present work, we calculated the distance between i) the Cα atom of Pro256 in the M2-M3 loop and the center of mass of the Cα atoms in the Cys-loop (residues 122–136); ii) the Cα atom of Pro256 and the

<figure><img src="image_14.png" /><figcaption>**Fig 11. The LBD-TMD interface.** Left panel: Cys-loop (orange), β1-β2 loop (yellow) and M2-M3 loop (red) in one subunit in the initial configuration. Pro256 belonging to the M2-M3 loop is also shown (red) in vdW representation. Central panel: distribution of the distances between the Cα atom of Pro256 on the M2-M3 loop and the center of mass of the Cys loop (Cα atoms of residues 122–136); right panel: distribution of the distances between the Cα atom of Pro256 on the M2-M3 loop and the Cα of residue Lys40 on the β1-β2 loop. Data collected from the five subunits. The vertical dashed blue line represents the value in the starting homology modeled configuration. Black bars: unrestrained; red bars: restrained; green bars: free simulations (data averaged over the four free trajectories).</figcaption></figure>

doi:10.1371/journal.pone.0133011.g011
Cα atom of Lys40 in the β1-β2 loop. This choice has been done by careful comparison of our protein primary sequence with other models in the literature. Indeed, in the case of native GluCl [53], it was monitored just the distance between the Cα of Pro268 on the M2-M3 loop and the Cα of VAL45 on the β1-β2 loop; in the case of the α4β2 human nAChR modeled in Ref. [50], two distances were monitored, between the Cα of Asp268 on the M2-M3 loop and the Cα of Arg48 on the β1-β2 loop; between the Cα of Lys274 on the M2-M3 loop and the Cα of Asp140 on the Cys loop, respectively. Distributions of the distance values are shown in Fig 11, panel B and C respectively, for both the collapsed and open forms. In the latter case, the M2-M3 loop detaches from the β1-β2 loop (Fig 11C), and stably resides closer to the Cys-loop (Fig 11B). In good agreement with results for open forms reported in Refs. [50, 53], results point out that both the restrained and free simulations provide conformations representative of the full-length channel in an open form.

## Intra protein hydrogen bonds

A transition from a *tense* to a more *relaxed* form of the channel has been invoked recently [10] in the description of the closed-to-open mechanism of *Torpedo* nicotinic receptor; in particular it was proposed that, due to the release of interhelical contacts, the M2 helices become free to adopt a configuration more asymmetrical in the open form, with respect to the closed structure. In that respect, intra protein hydrogen bond analysis, performed along all MD trajectories, provided a quite interesting qualitative picture. A hydrogen bond is considered to be formed when the donor-acceptor distance is lower than 3.0Å. We detected many more hydrogen bonds along the unrestrained trajectory than in the restrained or free ones, both in the LBD and TMD. In particular two hydrogen bonds are found between the TMD M2-M3 loop and the β1-β2 loop in subunit P4; three hydrogen bonds between the M2-M3 loop in the subunit P1 and the β10 in the subunit P2; two hydrogen bonds between the M2-M3 loop in the subunit P1 and the Cys-loop in the subunit P2. Two other intersubunit hydrogen bonds are formed between the M3 helix in the subunit P3 and the Cys-loop in the subunit P4. Quite interestingly, four intra and one intersubunit hydrogen bond were detected involving the TMD helices: one hydrogen bond between the helices M1 and M3 in the subunits P1; three hydrogen bonds between the helices M2 and M4 in the subunit P3; and one hydrogen bond between the helix M1 in the subunit P5 and the helix M3 in the subunit P1. All the hydrogen bonds listed above are persistently formed in time (percentage >60%) in the collapsed structure, while totally absent or scarcely persistent (<5%) in the open structures. Accordingly, the collapsed structure results more tense than the open forms.

## C loops conformation and interaction with the agonist epibatidine

In nicotinic receptors the ligand pocket, present in each of the subunits, is delimited by the so-called loop C and extends at the interface between subunits. The ligand type influences the degree of C-loop closure against the protein core [25, 93, 94], as the loop arrangements cluster into three groups: i) the agonist-bound “open” conformation; ii) the apo “intermediate” conformations; iii) the antagonist-bound “closed” conformations. To further assess our α7 model we analyzed the ligand pocket by describing the behavior of the C loop capping the epibatidine molecules using the distance ($d_{C}^{intra}$) between the C loop center of mass and the center of mass of a couple of residues in the backwall. Protein-ligand interactions have been also monitored along the MD trajectories of the complexed system. According to crystallography [26], we selected and monitored three electrostatic interactions: the first between the epibatidine bridge ring amine and the Tyr87 side chain oxygen; the second between the epibatidine pyrimidine amine and the Tyr189 side chain oxygen; the third between the epibatidine pyrimidine amine
<figure><img src="image_15.png" /><figcaption>**Fig 12. Protein-ligand interaction.** A) Epibatidine and key residues (in licorice) in the binding pocket, in the homology model. The primary and complementary subunits are shown in blue and yellow cartoon representation, respectively. Epibatidine is in green. Chlorine (purple) and nitrogen (blue) atoms on epibatidine in ball representation. B) Representative snapshot along the Free<sub>3</sub> trajectory with epibatidine in axial position, with the chlorine atom down; C) in equatorial position; D) epibatidine and key residues in the complementary subunit. E) distribution of d<sup>intra</sup><sub>C</sub>. Data collected from the five subunits. The vertical dashed blue line represents the value in the starting homology modeled configuration. Black bars: unrestrained; red bars: restrained; green bars: free simulations (data averaged over the four free trajectories).</figcaption></figure>

doi:10.1371/journal.pone.0133011.g012

and the Trp143 carbonyl oxygen. We also considered two other electrostatic interactions, those between the epibatidine aromatic chlorine and each of the two backbone nitrogens of Leu103 and Leu113 from the complementary subunit. Finally, we analyzed three van der Waals interactions: between the epibatidine pyrimidine ring and the Trp143 indole; between the epibatidine bridge ring and the Tyr189 aromatic ring; between the epibatidine bridge ring and the Tyr182 aromatic ring. The distance between the involved charged atom pairs has been calculated along the whole trajectories; the van der Waals interactions have been monitored through the distance between the center of mass of the bridge ring or of the pyrimidine ring and the center of mass of the residues aromatic rings. The epibatidine and the key residues above mentioned are shown in Fig 12, panel A. Distances averaged along the six trajectories are reported in Table 2. Results show that, on average, the ligand is kept in place during the restrained and free simulations.

The comparison with the unrestrained run suggests that the most of ligand/binding pocket interactions weaken if the channel collapses. This is confirmed for almost all subunits by an

<table><caption>**Table 2. Protein-ligand distances.**</caption><thead><tr><th>Distance</th><th>Unrestrained</th><th>Restrained</th><th>Free<sub>1</sub></th><th>Free<sub>2</sub></th><th>Free<sub>3</sub></th><th>Free<sub>4</sub></th></tr></thead><tbody><tr><td>EPI-N1/TRP143-O</td><td>6.2 ± 1.7</td><td>4.5 ± 0.9</td><td>4.7 ± 1.0</td><td>4.4 ± 0.9</td><td>5.1 ± 1.5</td><td>4.5 ± 0.9</td></tr><tr><td>EPI-N1/TYR189-OH</td><td>6.5 ± 2.4</td><td>5.3 ± 1.1</td><td>5.9 ± 1.4</td><td>4.9 ± 0.5</td><td>6.2 ± 2.1</td><td>5.3 ± 1.1</td></tr><tr><td>EPI-N2/TYR87-OH</td><td>6.2 ± 1.3</td><td>5.9 ± 2.0</td><td>6.0 ± 1.3</td><td>5.9 ± 1.2</td><td>6.6 ± 1.8</td><td>5.9 ± 2.0</td></tr><tr><td>EPI-CL/LEU103-N</td><td>8.9 ± 5.2</td><td>4.2 ± 0.8</td><td>4.1 ± 0.8</td><td>4.3 ± 0.8</td><td>5.4 ± 2.6</td><td>5.1 ± 1.2</td></tr><tr><td>EPI-CL/LEU113-N</td><td>8.2 ± 3.0</td><td>5.1 ± 1.2</td><td>5.3 ± 1.3</td><td>5.8 ± 1.2</td><td>6.5 ± 2.1</td><td>5.6 ± 1.6</td></tr><tr><td>EPI-pyrimidine/TRP143-indole</td><td>6.1 ± 1.3</td><td>5.0 ± 0.7</td><td>5.1 ± 0.8</td><td>5.3 ± 0.8</td><td>5.2 ± 0.9</td><td>4.9 ± 0.7</td></tr><tr><td>EPI-bridge ring/TYR189-ring</td><td>6.3 ± 2.5</td><td>5.3 ± 1.1</td><td>5.0 ± 0.7</td><td>4.9 ± 0.5</td><td>6.0 ± 2.2</td><td>5.5 ± 1.7</td></tr><tr><td>EPI-bridge ring/TYR182-ring</td><td>5.3 ± 0.7</td><td>5.2 ± 0.4</td><td>5.3 ± 0.4</td><td>5.4 ± 0.4</td><td>5.8 ± 1.2</td><td>5.6 ± 0.8</td></tr></tbody></table>

Distances (in Å) between epibatidine and selected residues (see Section “C loops conformation and interaction with the agonist epibatidine”). Values averaged along 200ns, using data collected from the five subunits.

doi:10.1371/journal.pone.0133011.t002
<table><caption>Table 3. Protein-ligand interaction energies.</caption><thead><tr><th>Subunit</th><th>Trajectory</th><th>van der Waals</th><th>electrostatic</th></tr></thead><tbody><tr><td rowspan="2">P1</td><td>Restrained</td><td>-15.7</td><td>-1.2</td></tr><tr><td>Unrestrained</td><td>-10.9</td><td>-1.6</td></tr><tr><td rowspan="2">P2</td><td>Restrained</td><td>-13.6</td><td>-2.3</td></tr><tr><td>Unrestrained</td><td>-11.8</td><td>-0.1</td></tr><tr><td rowspan="2">P3</td><td>Restrained</td><td>-14.9</td><td>-4.1</td></tr><tr><td>Unrestrained</td><td>-13.7</td><td>-1.7</td></tr><tr><td rowspan="2">P4</td><td>Restrained</td><td>-12.8</td><td>-3.8</td></tr><tr><td>Unrestrained</td><td>-18.6</td><td>-1.5</td></tr><tr><td rowspan="2">P5</td><td>Restrained</td><td>-12.3</td><td>-1.5</td></tr><tr><td>Unrestrained</td><td>-11.2</td><td>-0.9</td></tr></tbody></table>

Interaction energies (in kcal/mol) between epibatidine and selected residues (see Section “C loops conformation and interaction with the agonist epibatidine”). Values averaged along the 100–200ns portion of the unrestrained or restrained trajectory. The error bar is within 3% in all cases.

doi:10.1371/journal.pone.0133011.t003

analysis of van der Waals and electrostatic interaction energies between the ligand and the pocket residues above mentioned (Table 3).

As shown in Fig 12, panel E, the C-loop/backwall distance is, on average, 14.6 Å in the homology model; it increases to values between 15 Å and 16 Å in the open conformations. In particular, for the restrained trajectory the distribution shows peak maxima at 15.2, 15.9 and a shoulder at 16.5 Å; the open free trajectories provide a bimodal distribution, with peaks at 15.1 and 16.1 Å. This behavior corresponds to the different positions of the ligand in the binding pocket, as also observed in other bound ligand receptors [95, 96]. Indeed, in some subunits, the epibatidine undergoes several transitions from an axial to an equatorial orientation, i.e. pointing the chlorine atom up toward the apical side, as it was in the homology model (Fig 12, panel A), or flipped down (as in Fig 12, panel B) or radially toward the pore center (Fig 12, panel C). In the distributions of $d_{C}^{intra}$ the peaks at 15.2 and 15.1 Å are related to the ligand axial orientation, and the ones at 16.5 and 16.1 Å to the equatorial orientation. The shoulder at 15.9 Å in the restrained distribution is due to a position intermediate between the two. In one case, at the end of the free$_{3}$ trajectory, the epibatidine detaches from the primary subunit, sliding toward the complementary subunit, where interactions with charged and aromatic side chains (e.g. Tyr26, Ser53, Thr55, Gln155) compete with the ones in the primary subunit to keep the ligand in place at the interface (Fig 12, panel D).

Discussion

The knowledge of few high resolution crystal structures together with atomistic simulations, has recently started to yield valuable insight on the gating mechanism of nAChRs in response to ligand binding [13, 53–55, 97]. However, open-pore conformations of nAChRs are still scarcely available. In this respect, the X-ray structure determination of the prokaryotic pentameric GLIC channel, claimed in open form [15], has contributed to the understanding of the structure/dynamics/function relationships in nAChRs. However, stability assessment via full atomistic MD pointed out that GLIC undergoes a spontaneous closure [79, 97] and that pH variations are needed to stabilize the channel in the open form [11, 13]. In one study only, GLIC was simulated for 200ns without observing pore collapse [98].

With respect to the electron microscopy structures used in the past for modeling hetero- and homopentamers like the human α7, the X-ray structure of the glutamate receptor GluCl
[41, 42] is undoubtedly a better template. Besides the higher atomic resolution, at variance with GLIC, GluCl is activated by ligands distinct from the solvent, and therefore much more suitable to study ligand-gated ion channels. Concerns have been raised however on the open channel stability of the X-ray structure of GluCl: recent, independent MD simulations [53–55] point out that the experimental open structure collapses towards a non-conductive, dehydrated channel; it has been shown that sustaining the open form is very difficult, already on the tens of nanoseconds time scale. This in principle could hamper its use for molecular modeling of the nicotinic receptors transition. Normal mode analysis has been used [55] in conjunction with an elastic network model of the protein to produce an open structure based on native GluCl; the putative open conformation was assessed along 100ns of MD. In this model, the M2 helices motion in the TMD controls the channel hydration, through an overall backbone and side-chain reorientation; the polar and azimuthal tilt angles were found to be 5.25° and 1.0°, respectively. The radius of the obtained open pore was 5 Å.

The present work adds to the field, as, to the best of our knowledge, it is the first to use GluCl as template for a human nAChR. We created a chimeric model using the structure of GluCl [41] to model the TMD of the channel and the Cys-loop region, plus a known structure of AChBP to model the ECD of the channel bound to the agonist epibatidine [26]. The structure obtained can be considered as a model of the full-length human α7 pentamer in the active state both at the TMD and LBD level.

We have presented a combination of homology modeling and full-atomistic, explicit solvent MD simulations of the human nicotinic receptor α7. It is well known that the accuracy of a computational model is determined mostly by two factors [62]: i) the sequence identity between the target and the template, and ii) the quality of the template. Although the percentage of identity between our two sequences is low (25.37%), we obtained an optimal alignment of the functional regions and only few gaps. As for the quality of the template, we used crystallographic structures with the highest resolution available (2BYQ: 3.4 Å and 3RIF: 3.35 Å), as described in Materials and Methods. Unbiased MD in a lipid bilayer/water system was performed to assess the stability of the model and to refine it at high-resolution, in light of a well-established approach of combined static modeling and dynamical optimization [99–101]. Thus, our refinement is hampered only by the known limitations of MD simulations, namely the quality of the force field and the restricted conformational sampling [102]. Results show that the structure rapidly evolve towards a collapsed state, in agreement with the recent literature [53–55]. At variance with using normal modes to force the closed channel to assume the open form [55], in the present study we performed MD simulations by perturbing the chimeric model to maintain the channel pore open, using a restraint on the distances between M2 helices designed to switch on only when the helices tend to come close to each other. Remarkably, such local restraint on few specific, functional residues let us observe the spontaneous change of other, not directly affected regions towards the stable structure with a wide pore.

Overall, such structure was preserved for a cumulative simulation time of about 1 µs, 800ns of which are unbiased, free trajectories. Based mostly on the analysis of the level of hydration, pore size and profile (see Fig 8 and S8 Fig) and supported from the comparison with other putative open forms reported in the literature [15, 50, 55], we suggest that the obtained structure has reasonable chances to be (at least one of the possible) structures of the channel in the open conformation.

An unsolved issue about channel receptors concerns the different states of the receptor described by the electro-physiological experiments: open, closed, resting, desensitized. There is no clear correspondence between these experimental observations and simple structural criteria, beside the C loop opening/closure in the LBD. This structural feature was considered in the past as related to the type of ligand, but this is in contradiction with recent results on a AChBP
structure bound to dihydro-beta-erythroidine (DHβE)[103], which is an antagonist, but so small that the C loop does not open more than in presence of an agonist.

It is known that the M2 helices play a role in channel gating but parameters describing rigid-body motions of the M2 helices such as helices tilting, or the quaternary LBD/TMD twist, are not sufficient to describe the complex motions in the pore. Discerning open from closed structures based on these structural quantities can be an ambiguous task, as already observed for both heteropentameric [50] and homopentameric [46] models of human nAChRs. In particular, it has been reported that the M2 helices undergo a non- concerted motion, with no correlation between the radial and lateral tilting [50], in which channel opening is driven by the changes in the sidechain orientations of few hydrophobic pore-lining residues, similarly to what we find in this work. In this respect, the analysis of the pore profile and hydration is much more informative. Our α7 nAChR we model has a pore consistent with an open channel, slightly larger than the pore in the open model of α4β2 nAChR [50], the modeled system closest to the one modeled here. Furthermore, the comparison with hydrophobic pore models and GLIC [83, 88] points out that the stationary hydration level in our putative open structure is consistent with the one of a conductive channel. Atomic level details of the TMD-LBD interface, the protein hydrogen bonds network, and the epibatidine conformations at the binding sites contribute to and support the characterization of the active channel both in the extracellular and transmembrane regions.

To conclude, this work provides a first step towards obtaining a reliable and stable all-atom model of the channel in an open form. Recently it has been observed the existence of multiple conformations for a given state, in GLIC [104], Torpedo nAChR [10] and in ELIC [105]. In this respect, the structure we found could be considered a good candidate for a distribution of conformations adequate to the functional open state. Our results potentially open multiple lines of investigation; indeed, the new structure could be used for studying the nAChR gating transition in path-sampling simulations, where the reliability of end-point structures is a fundamental requirement [106–108], or in virtual screening studies to dock different ligands.

## Supporting Information

**S1 Fig. The epibatidine.** Schematic structure of the epibatidine molecule. (TIF)

**S2 Fig. Crossed distances.** Time evolution of the crossed distances between Val13' (left panels) and Leu16' (right panels) centers of mass, in facing subunits, for each pair of subunits, along the restrained simulations. The horizontal red lines mark the values calculated on the initial structure, i.e. the threshold values used in the implemented flat-bottom quadratic restraint. (TIF)

**S3 Fig. Stability assessment.** Left panels: Root Mean Square Deviation (in Å) of individual subunits calculated from the starting conformations (after equilibration) as

$$ RMSD(t_j) = \sqrt{\frac{\sum_{C_\alpha=1}^{N_\alpha} (\mathbf{r}_{C_\alpha}(t_j) - \mathbf{r}_{C_\alpha}^{\text{ref}})^2}{N_\alpha}}, $$ where $\mathbf{r}_{C_\alpha}(t_j)$ is the position of the $C_\alpha$ atom at jth time step, $\mathbf{r}_{C_\alpha}^{\text{ref}}$ is the position of the $C_\alpha$ atom in the reference structure, $N_\alpha$ is the total number of $C_\alpha$ atoms in the subunit. The RMSD values are calculated after removing the roto-translational body motions of the single subunits [109]. Upper row panels: unrestrained; lower row panels: restrained trajectory. Right panels: Root Mean Square Fluctuations (RMSFs) (in Å) of $C_\alpha$ carbon atoms in the individual subunits with respect to the average structure calculated along the final 100ns of the trajectories. The curves are colored according to the scheme in the inset. The bold black line in the RMSFs plot is the average over the five subunits. Upper row panels:
unrestrained; lower row panels: restrained trajectory.

(TIF)

**S4 Fig. Stability assessment.** Left panels: $C_α$ atoms Root Mean Square Deviation (in Å) calculated from the corresponding starting conformation after removing the roto-translational body motions [109], along the four free trajectories. Right panels: Root Mean Square Fluctuations (RMSFs) (in Å) of the protein $C_α$ atoms with respect to the average structure calculated along the final 100ns of the trajectory, for the four free trajectories. Both RMSD and RMSFs values are averaged over the five subunits.

(TIFF)

**S5 Fig. Crossed distances.** Time evolution of the crossed distances between Val13' (upper panel) and Leu16' (lower panel) $C_α$ atoms in facing subunits, for each pair of subunits, in the first 20ns of MD including the equilibration stage. Solid lines: along the restrained simulation; dotted line: along the unrestrained simulation. The curves are colored according to the scheme in the inset in S3 Fig.

(TIF)

**S6 Fig. Hydration behavior.** Top panel: time evolution of water count in the pore lumen delimited by the M2 helices. Bottom panel: time evolution of water count in a region of 10 Å centered at the constriction point 13' (Val246). Black curve: unrestrained; red curve: restrained trajectory. Note dewetting-rewetting transitions occurring along the trajectory Free₃ (violet curve) at ~ 115ns, ~ 135ns, ~ 145ns, ~ 180ns.

(TIF)

**S7 Fig. Hydration behavior.** Representative snapshots of water molecule arrangements inside the pore at two different times along the Free₄ trajectory (in the portion 145–150ns). Upper row: water molecules within 3Å from any atom of Ser10' residues. Lower row: water molecules within 3Å from any atom of Thr6' residues. All Ser10' and Thr6' residues are in vdW representation, colored in cyan.

(TIFF)

**S8 Fig. Pore profiles.** Pore radius profiles along the pore region in configurations averaged over the last 100ns of the four free MD trajectories. Also reported the profile for the collapsed channel. Pore radius computed with MolAxis [80].

(TIF)

**S9 Fig. Ion distributions.** Ion distributions in the free trajectories. Blue bars: cations, that can permeate the nicotinic channel; red bars: anions. Top panels show the channel pore profiles along the full-length protein. The cyan region highlights the channel between Glu20' and Leu9'.

(TIF)

**S1 Text. Epibatidine parameters.** Topology and epibatidine partial atomic charges.

(PDF)

## Acknowledgments

We acknowledge CINECA award under the ISCRA initiative and SFI/HEA Irish Center for High-End Computing (ICHEC), for the availability of high performance computing resources and support; Sergio Giuffrida for his valuable help in the parametrization of the ligand force field terms; Daniele Di Marino for useful discussions. This work was supported by SFI-Science Foundation of Ireland (grant 08-IN.1-I1869) and Istituto Italiano di Tecnologia (IIT, seed
project grant N. 259SIMBEDD). TM acknowledges CNRS and Institut Pasteur for financial support and Dr. P. J. Corringer for fruitful discussions.

## Author Contributions

Conceived and designed the experiments: TEM LM G. Cottone. Performed the experiments: LC TEM G. Cottone LM. Analyzed the data: LC TEM G. Cottone LM. Contributed reagents/ materials/analysis tools: LC TEM G. Cottone LM. Wrote the paper: G. Cottone LC TEM LM G. Ciccotti.

## References

1. Taly A, Corringer PJ, Guedin D, Lestage P, Changeux JP. Nicotinic receptors: allosteric transitions and therapeutic targets in the nervous system. Nat Rev Drug Discov. 2009; 8:733–750. doi: 10.1038/nrd2927 PMID: 19721446

2. Changeux JP, Edelstein SJ. Nicotinic Acetylcholine Receptors: From Molecular Biology to Cognition. New York: Odile Jacob, 1st ed.; 2005.

3. Changeux JP, Edelstein SJ. Allosteric mechanisms of signal transduction. Science. 2005; 308:1424–1428. doi: 10.1126/science.1108595 PMID: 15933191

4. Nys M, Kesters D, Ulens C. Structural insights into Cys-loop receptor function and ligand recognition. Biochem Pharm. 2013; 86:1042–1053. doi: 10.1016/j.bcp.2013.07.001 PMID: 23850718

5. Whatley VJ, Harris RA. The cytoskeleton and neurotransmitter receptors. Int Rev Neurobiol. 1996; 39:113–143. doi: 10.1016/S0074-7742(08)60665-0 PMID: 8894846

6. Larsson A, Engel JA. Neurochemical and behavioral studies on ethanol and nicotine interactions. Neurosci Biobehav Rev. 2004; 27:713–720. doi: 10.1016/j.neubiorev.2003.11.010 PMID: 15019421

7. Bordarenko V, Mowrey D, Tillman T, Seyoum E, Xe Y, Tang P. NMR structures of the human α7 nAChR transmembrane domain and associated anesthetic binding sites. Biochim Biophys Acta. 2014; 1838:1389–1395. doi: 10.1016/j.bbamem.2013.12.018

8. Auerbach A. The Energy and Work of a Ligand-Gated Ion Channel. J Mol Biol. 2013; 425:1461–1475. doi: 10.1016/j.jmb.2013.01.027 PMID: 23357172

9. Unwin N. Refined structure of the nicotinic acetylcholine receptor at 4 Å resolution. J Mol Biol. 2005; 346:967–989. doi: 10.1016/j.jmb.2004.12.031 PMID: 15701510

10. Unwin N, Fujiyoshi Y. Gating movement of acetylcholine receptor caught by plunge freezing. J Mol Biol. 2012; 422:617–634. doi: 10.1016/j.jmb.2012.07.010 PMID: 22841691

11. Nury H, Renterghem CV, Weng Y, Tran A, Baaden M, Dufresne V, et al. X-ray structures of general anaesthetics bound to a pentameric ligand-gated ion channel. Nature. 2011; 469:428–431. doi: 10.1038/nature09647 PMID: 21248852

12. Hilf RJ, Bertozzi C, Zimmermann I, Reiter A, Trauner D, Dutzler R. Structural basis of open channel block in a prokaryotic pentameric ligand-gated ion channel. Nat Struct Mol Biol. 2010; 17:1330–1336. doi: 10.1038/nsmb.1933 PMID: 21037567

13. Nury H, Poitevin F, Renterghem CV, Changeux JP, Corringer PJ, Delarue M, et al. One-microsecond molecular dynamics simulation of channel gating in a nicotinic receptor homologue. Proc Natl Acad Sci USA. 2010; 107:6275–6280. doi: 10.1073/pnas.1001832107 PMID: 20308576

14. Hilf RJ, Dutzler R. Structure of a potentially open state of a proton-activated pentameric ligand-gated ion channel. Nature. 2009; 457:115–118. doi: 10.1038/nature07461 PMID: 18987630

15. Bocquet N, Nury H, Baaden M, Poupon CL, Changeux JP, Delarue M, et al. X-ray structure of a pentameric ligand-gated ion channel in an apparently open conformation. Nature. 2009; 457:111–114. doi: 10.1038/nature07462 PMID: 18987633

16. Zimmermann I, Dutzler R. Ligand activation of the prokaryotic pentameric ligand-gated ion channel ELIC. PLoS Biol. 2011; 9:e1001101. doi: 10.1371/journal.pbio.1001101 PMID: 21713033

17. Hilf RJ, Dutzler R. X-ray structure of a prokaryotic pentameric ligand-gated ion channel. Nature. 2008; 452:375–379. doi: 10.1038/nature06717 PMID: 18322461

18. Pan J, Chen Q, Willenbring D, Yoshida K, Tillman T, Kashlan OB, et al. Structure of the Pentameric Ligand-Gated Ion Channel ELIC Co-Crystallized with its Competitive Antagonist Acetylcholine. Nat Commun. 2012; 3:714. doi: 10.1038/ncomms1703 PMID: 22395605
19. Gonzalez-Gutierrez G, Cuello L, Nair S, Grosman C. Gating of the proton-gated ion channel from Gloeobacter violaceus at pH 4 as revealed by X-ray crystallography. Proc Natl Acad Sci USA. 2013; 110:18716-18721. doi: 10.1073/pnas.1313156110 PMID: 24167270

20. Prevost MS, Sauguet L, Nury H, Renterghem CV, Huon C, Poitevin F, et al. A locally closed conformation of a bacterial pentameric proton-gated ion channel. Nat Struct Mol Biol. 2012; 19:642-649. doi: 10.1038/nsmb.2307 PMID: 22580559

21. Celie PH, van Rossum-Fikkert SE, van Dijk WJ, Brejc K, Smit AB, Sixma K. Nicotine and carbamyl-choline binding to nicotinic acetylcholine receptors as studied in AChBP crystal structures. Neuron. 2004; 41:907-914. doi: 10.1016/S0896-6273(04)00115-1 PMID: 15046723

22. Celie PH, Kasheverov IE, Mordvintsev DY, Hogg RC, van Nierop P, van Elk R, et al. Crystal structure of nicotinic acetylcholine receptor homolog AChBP in complex with an alpha-conotoxin PnIA variant. Nat Struct Biol. 2005; 12:582-588. doi: 10.1038/nsmb951

23. Celie PH, Klaassen RV, van Rossum-Fikkert SE, van Elk R, van Nierop P, Smit AB, et al. Crystal structure of acetylcholine-binding protein from *Bulinus truncatus* reveals the conserved structural scaffold and sites of variation in nicotinic acetylcholine receptors. J Biol Chem. 2005; 280:26457-26466. doi: 10.1074/jbc.M414476200 PMID: 15899893

24. Hibbs RE, Radic Z, Taylor P, Johnson DA. Influence of agonists and antagonists on the segmental motion of residues near the agonist binding pocket of the acetylcholine-binding protein. J Biol Chem. 2006; 281:39708-39718. doi: 10.1074/jbc.M604752200 PMID: 17068341

25. Hibbs RE, Sulzenbacher G, Shi J, Talley TT, Conrod S, Kem WR, et al. Structural determinants for interaction of partial agonists with acetylcholine binding protein and neuronal alpha7 nicotinic acetyl-choline receptor. EMBO J. 2009; 28:3040-3051. doi: 10.1038/emboj.2009.227 PMID: 19696737

26. Hansen SB, Sulzenbacher G, Huxford T, Marchot P, Taylor P, Bourne Y. Structures of Aplysia AChBP Complexes with Nicotinic Agonists and Antagonists Reveal Distinctive Binding Interfaces and Conformations. EMBO J. 2005; 24:3635-3646. doi: 10.1038/sj.emboj.7600828 PMID: 16193063

27. Ulens C, Akdemir A, Jongejan A, van Elk R, Bertrand S, Perrakis A, et al. Use of acetylcholine binding protein in the search for novel alpha7 nicotinic receptor ligands. In silico docking, pharmacological screening, and X-ray analysis. J Med Chem. 2009; 52:2372-2383. doi: 10.1021/jm801400g PMID: 19331415

28. Ulens C, Hogg RC, Celie PH, Bertrand D, Tsetlin V, Smit AB, et al. Structural determinants of selective alpha-conotoxin binding to a nicotinic acetylcholine receptor homolog AChBP. Proc Natl Acad Sci USA. 2006; 103:3615-3620. doi: 10.1073/pnas.0507889103 PMID: 16505382

29. Bourne Y, Hansen SB, Sulzenbacher G, Talley TT, Huxford T, Taylor P, et al. Structural comparison of three crystalline complexes of a peptidic toxin with a synaptic acetylcholine recognition protein. J Mol Neurosci. 2006; 30:103-104. doi: 10.1385/JMNN:30:1:103 PMID: 17192648

30. Bourne Y, Radic Z, Aráoz R, Talley TT, Benoit E, Servent D, et al. Structural determinants in phycotoxins and AChBP conferring high affinity binding and nicotinic AChR antagonism. Proc Natl Acad Sci USA. 2010; 107:6076-6081. doi: 10.1073/pnas.0912372107 PMID: 20224036

31. Bourne Y, Talley TT, Hansen SB, Taylor P, Marchot P. Crystal structure of a CBTX-AChBP complex reveals essential interactions between snake alpha-neurotoxins and nicotinic receptors. EMBO J. 2005; 24:112-115. doi: 10.1038/sj.emboj.7600620

32. Hansen SB, Taylor P. Galanthamine and non-competitive inhibitor binding to ACh-binding protein: evidence for a binding site on non-alpha-subunit interfaces of heteromeric neuronal nicotinic receptors. J Mol Biol. 2007; 369:895-890. doi: 10.1016/j.jmb.2007.03.067 PMID: 17481657

33. Dutertre S, Ulens C, Buttner R, Fish A, van Elk R, Kendel Y, et al. AChBP-targeted alpha-conotoxin correlates distinct binding orientations with nAChR subtype selectivity. EMBO J. 2007; 26:3858-3860. doi: 10.1038/sj.emboj.7601785 PMID: 17660751

34. Edink E, Rucktooa P, Retra K, Akdemir A, Nahar T, Zuiderveld O, et al. Fragment Growing Induces Conformational Changes in Acetylcholine-Binding Protein: A Structural and Thermodynamic Analysis. J Am Chem Soc. 2011; 133:5363-5371. doi: 10.1021/ja110571r PMID: 21322593

35. Rucktooa P, Smit AB, Sixma TK. Insight in nAChR subtype selectivity from AChBP crystal structures. Biochem Pharmacol. 2009; 78:777-787. doi: 10.1016/j.bcp.2009.06.098 PMID: 19576182

36. Tsetlin V, Hucho F. Nicotinic acetylcholine receptors at atomic resolution. Curr Opin Pharmacol. 2009; 9:306-310. doi: 10.1016/j.coph.2009.03.005 PMID: 19428299

37. Costa V, Nistri A, Cavalli A, Carloni P. A structural model of agonist binding to the α3β4 neuronal nicotinic receptor. Br J Pharmacol. 2003; 140:921-931. doi: 10.1038/sj.bjp.0705498 PMID: 14504134

38. Gao F, Bren N, Burghardt TP, Hansen S, Henchman RH, Taylor P, et al. Agonist-mediated conformational changes in Acetylcholine-binding Protein revealed by simulation and Intrinsic Tryptophan Fluorescence. J Biol Chem. 2005; 280:8443-8451. doi: 10.1074/jbc.M412389200 PMID: 15591050
39. Amiri S, Sansom MS, Biggin PC. Molecular dynamics studies of AChBP with nicotine and carbamyl-choline: the role of water in the binding pocket. Protein Eng Des Sel. 2007; 20:353-359. doi: 10.1093/protein/gzm029 PMID: 17595341

40. Mohammad Hosseini Naveh Z, Malliavin TE, Maragliano L, Cottone G, Ciccotti G. Conformational Changes in Acetylcholine Binding Protein Investigated by Temperature Accelerated Molecular Dynamics. PLoS ONE. 2014; 9:e88555. doi: 10.1371/journal.pone.0088555 PMID: 24551117

41. Hibbs R, Gouaux E. Principles of activation and permeation in an anion-selective Cys-loop receptor. Nature. 2011; 474:54-60. doi: 10.1038/nature10139 PMID: 21572436

42. Althoff T, Hibbs RE, Banerjee S, Gouaux E. X-ray structures of GluCl in apo states reveal a gating mechanism of Cys-loop receptors. Nature. 2014; 512:333-337. doi: 10.1038/nature13669 PMID: 25143115

43. Miller PS, Aricescu AR. Crystal structure of a human GABA(A) receptor. Nature. 2014; 512:270-275. doi: 10.1038/nature13293 PMID: 24909990

44. Hassaine G, Deluz C, Grasso L, Wyss R, Tol MB, Hovius R, et al. X-ray structure of the mouse serotinin 5-HT3 receptor. Nature. 2014; 512:276-281. doi: 10.1038/nature13552 PMID: 25119048

45. Law RJ, Henchman RH, McCammon JA. A gating mechanism proposed from a simulation of a human alpha7 nicotinic acetylcholine receptor. Proc Natl Acad Sci USA. 2005; 102:6813-6818. doi: 10.1073/pnas.0407739102 PMID: 15857954

46. Cheng X, Wang H, Grant B, Sine SM, McCammon JA. Targeted molecular dynamics study of C-loop closure and channel gating in nicotinic receptors. PLoS Comput Biol. 2006; 2:1173-1184. doi: 10.1371/journal.pcbi.0020134

47. Cheng X, Lu B, Grant B, Law R, McCammon J. Channel Opening Motion of alpha7 Nicotinic Acetylcholine Receptor as suggested by Normal Mode Analysis. J Mol Biol. 2006; 355:310-324. doi: 10.1016/j.jmb.2005.10.039 PMID: 16307758

48. Cheng X, Ivanov I, Wang H, Sine SM, McCammon JA. Nanosecond-timescale conformational dynamics of the human alpha7 nicotinic acetylcholine receptor. Biophys J. 2007; 93:2622-2634. doi: 10.1529/biophysj.107.109843 PMID: 17573436

49. Taly A, Delarue M, GRutter T, Nilges M, Novere NL, Corringer PJ, et al. Normal Mode Analysis Suggests a Quaternary Twist Model for the Nicotinic Receptor Gating Mechanism. Biophys J. 2005; 88:3954-3965. doi: 10.1529/biophysj.104.050229 PMID: 15805177

50. Haddadian EJ, Cheng MH, Coalson RD, Xu Y, Tang P. In Silico Models for the Human alpha 4 beta 2 Nicotinic Acetylcholine Receptor. J Phys Chem B. 2008; 112:13981-13990. doi: 10.1021/jp804868s PMID: 18847252

51. Bergmann R, Kongsbak K, SÃrensen P, Sander T, Balle T. NMR structures of the human α7 nAChR transmembrane domain and associated anesthetic binding sites. PLoS One. 2013; 8:e52323. doi: 10.1371/journal.pone.0052323 PMID: 23308109

52. YU R, Hurdiss E, Greiner T, Lape R, Sivilotti L, Biggin P. Agonist and Antagonist Binding in Human Glycine Receptors. Biochemistry. 2014; 53:6041-6051. doi: 10.1021/bi500815f PMID: 25184435

53. Calimet N, Simoes M, Changeux JP, Karplus M, Taly A, Cecchini M. A gating mechanism of penta-meric ligand-gated ion channels. Proc Natl Acad Sci USA. 2013; 110:3987-3996. doi: 10.1073/pnas.1313785110

54. Yoluk O, Bromstrup T, Bertaccini EJ, Trudell JR, Lindahl E. Stabilization of the GluCl Ligand-Gated Ion Channel in the Presence and Absence of Ivermectin. Biophys J. 2013; 105:640-647. doi: 10.1016/j.bpj.2013.06.037 PMID: 23931312

55. Cheng MH, Coalson RD. Energetics and Ion permeation Characteristics in a Glutamate-Gated Chloride (GluCl) Receptor Channel. J Phys Chem B. 2012; 116:13637-13643. doi: 10.1021/jp3074915 PMID: 23088363

56. Bouzat C, Gumilar F, Spitzmaul G, Wang HL, Rayes D, Hansen SB, et al. Coupling of agonist binding to channel gating in an ACh-binding protein linked to an ion channel. Nature. 2004; 430(7002):896-900. doi: 10.1038/nature02753 PMID: 15318223

57. Lee WY, Sine SM. Principal pathway coupling agonist binding to channel gating in nicotinic receptors. Nature. 2005; 438(7065):243-247. doi: 10.1038/nature04156 PMID: 16281039

58. Mukhtasimova N, Free C, Sine SM. Initial Coupling of Binding to Gating Mediated by Conserved Residues in the Muscle Nicotinic Receptor. J Gen Physiol. 2005; 126:23-39. doi: 10.1085/jgp.200509283 PMID: 15955875

59. Reeves DC, Jansen M, Bali M, Lemster T, Akabas MH. A Role for the β1-β2 Loop in the Gating of 5-HT3 Receptors. J Neurosci. 2005; 25:9358-9366. doi: 10.1523/JNEUROSCI.1045-05.2005 PMID: 16221844
60. NovÃ˜re NL, Grutter T, Changeux J. Models of the extracellular domain of the nicotinic receptors and of agonist- and Ca²⁺ binding sites. Proc Natl Acad Sci USA. 2002; 98:3210–3215.

61. Grazioso G, Cavalli A, Amici MD, Recanatini M, Micheli CD. Alpha7 nicotinic acetylcholine receptor agonists: prediction of their binding affinity through a molecular mechanics Poisson Boltzmann surface area approach. J Comput Chem. 2008; 29:2593–2602. doi: 10.1002/jcc.21019 PMID: 18478580

62. MartÃ-Renom MA, Stuart AC, Fiser A, SÃ¡nchez R, Melo F, Sali A. Comparative Protein Structure Modeling of Genes and Genomes. Annu Rev Bioph Biom. 2000; 29:291–325. doi: 10.1146/annurev.biophys.29.1.291

63. Forrest L, Tang C, Honig B. On the Accuracy of Homology Modeling and Sequence Alignment Methods Applied to Membrane Proteins. Biophys J. 2006; 91:508–517. doi: 10.1529/biophysj.106.082313 PMID: 16648166

64. Notredame C, Higgins D, Heringa J. T-Coffee: A novel method for fast and accurate multiple sequence alignment. J Mol Biol. 2000; 302:205–217. doi: 10.1006/jmbi.2000.4042 PMID: 10964570

65. Eswar N, Webb B, Marti-Renom MA, Madhusudhan MS, Eramian D, Shen MY, et al. Comparative Protein Structure Modeling using Modeller. Current Protoc Bioinformatics. 2006; 5:1–30.

66. Laskowski R, MacArthur M, Moss D, Thornton J. PROCHECK: a program to check the stereochemical quality of protein structure. J Appl Cryst. 1993; 26:283–291. doi: 10.1107/S0021889892009944

67. Chen V, Arendall W, Headd J, Keedy D, Immormino R, Kapral G, et al. MolProbity: all-atom structure validation for macromolecular crystallography. Acta Crystallographica D. 2010; 66:12–21. doi: 10.1107/S0907444909042073

68. Lovell S, Davis I, Arendall W, de Bakker P, Word J, Prisant M, et al. Structure validation by Calpha geometry: phi,psi and Cbeta deviation. Proteins. 2003; 50:437–450. doi: 10.1002/prot.10286 PMID: 12557186

69. Benkert P, Künzli M, Schwede T. QMEAN Server for Protein Model Quality Estimation. Nucleic Acids Res. 2009; 37:510–514. doi: 10.1093/nar/gkp322

70. Delano WL. The PyMOL Molecular Graphics Systemml: http://www.pymol.org. DeLano Scientific, Palo Alto, CA, USA; 2002. Available fromml: http://www.pymol.org.

71. Bahsford D, Karplus M. pKa of Ionizable Groups in Proteins: Atomic Detail from a Continuum Electrostatic Model. Biochemistry. 1990; 29:10219–10225. doi: 10.1021/bi00496a010

72. Anandakrishnan R, Aguilar B, Onufriev A. H++ 3.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulation. Nucleic Acids Res. 2012; 40:537–541. doi: 10.1093/nar/gks375

73. Phillips J, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. J Comput Chem. 2005; 26:1781–1802. doi: 10.1002/jcc.20289 PMID: 16222654

74. Quigley D, Probert M. Langevin dynamics in constant pressure extended systems. J Chem Phys. 2004; 120:11432–11441. doi: 10.1063/1.1755657 PMID: 15268177

75. Wang J, Wang W, Kollmann P, Case D. Antechamber, An Accessory Software Package For Molecular Mechanical Calculation. J Comput Chem. 2005; 25:1157–1174. doi: 10.1002/jcc.20035

76. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al. Gaussian 03, Revision B.02. Pittsburgh PA: Gaussian, Inc.; 2003.

77. Essmann U, Perera L, Berkowitz M, Darden T, Lee H, Pedersen L. A smooth particle mesh Ewald method. J Chem Phys. 1995; 103:8577–8593. doi: 10.1063/1.470117

78. Ryckaert J, Ciccotti G, Berendsen H. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys. 1977; 23:327–341. doi: 10.1016/0021-9991(77)90098-5

79. LeBard D, Henin J, Eckenhoff R, Klein M, Brannigan G. General Anesthetics Predicted to Block the GLIC Pore with Micromolar Affinity. PLoS Comput Biol. 2012; 8:e1002532. doi: 10.1371/journal.pcbi.1002532 PMID: 22693438

80. Yaffe E, Fishelovitch D, Wolfson H, Halperin D, Nussinov R. MolAxis: A server for Identification of Channels in Macromolecule. Nucleic Acids Research. 2008; 36:210–215. doi: 10.1093/nar/gkn223

81. Baptista-Hon DT, Deeb TZ, Lambert JJ, Peters JA, Hales TG. The minimum M3-M4 loop length of neurotransmitter-activated pentameric receptors is critical for the structural integrity of cytoplasmic portals. J Biol Chem. 2013; 288:21558–21568. doi: 10.1074/jbc.M113.481689 PMID: 23740249

82. Henchman RH, Wang HL, Sine SM, Taylor P, McCammon JA. Asymmetric Structural Motions of the Homomeric alpha7 Nicotinic Receptor Ligand Binding Domain Revealed by Molecular Dynamics Simulations. Biophys J. 2003; 85:3007–3018. doi: 10.1016/S0006-3495(03)74720-1 PMID: 14581202
83. Beckstein O, Sansom MSP. The influence of geometry, surface character, and flexibility on the permeation of ions and water through biological pores. Phys Biol. 2004; 1:42–52. doi: 10.1088/1478-3967/1/1/005 PMID: 16204821

84. Murail S, Wallner B, Trudell JR, Bertaccini E, Lindahl E. Microsecond Simulations Indicate that Ethanol Binds between Subunits and Could Stabilize an Open-State Model of a Glycine Receptor. Biophys J. 2011; 100:1642–1650. doi: 10.1016/j.bpj.2011.02.032 PMID: 21463577

85. Taly A, Corringer PJ, Grutter T, Prado L, Karplus M, Changeux JP. Implications of the quaternary twist allosteric model for the physiology and pathology of nicotinic acetylcholine receptors. Proc Natl Acad Sci USA. 2006; 103:16965–16970. doi: 10.1073/pnas.0607477103 PMID: 17077146

86. Samson AO, Levitt M. Inhibition mechanism of the acetylcholine receptor by alpha-neurotoxins as revealed by normal-mode dynamics. Biochemistry. 2008; 47:4065–4070. doi: 10.1021/bi702272j PMID: 18327915

87. Hung A, Tai K, Sansom MS. Molecular dynamics simulation of the M2 helices within the nicotinic acetylcholine receptor transmembrane domain: structure and collective motions. Biophys J. 2005; 88 (5):3321–3333. doi: 10.1529/biophysj.104.052878 PMID: 15722430

88. Mowrey D, Cheng MH, Liu LT, Willenbring D, Lu X, Wymore T, et al. Asymmetric Ligand Binding Facilities Conformational Transitions in Pentameric Ligand-Gated Ion Channels. J Am Chem Soc. 2013; 135:2172–2180. doi: 10.1021/ja307275v PMID: 23339564

89. Song C, Corry B. Ion Conduction in Ligand-Gated Ion Channels: Brownian Dynamics Studies of Four Recent Crystal Structures. Biophysical Journal. 2010; 98(3):404–411. doi: 10.1016/j.bpj.2009.10.032 PMID: 20141753

90. Wang H, Cheng X, Taylor P, McCammon J, Sine S. Control of Cation Permeation through the Nicotinic Receptor Channel. PLoS Comput Biol. 2008; 4:e41. doi: 10.1371/journal.pcbi.0040041 PMID: 18282090

91. Sauguet L, Poitevin F, Murail S, Renterghem CV, Moraga-Cid G, Malherbe L, et al. Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels. EMBO J. 2013; 32:728–741. doi: 10.1038/emboj.2013.17 PMID: 23403925

92. Xiu X, Hanek A, Wang J. A unified view of the role of electrostatic interactions in modulating the gating of cys loop receptors. J Biol Chem. 2005; 280:41655–41666. doi: 10.1074/jbc.M508635200 PMID: 16216879

93. Li SX, Huang S, Bren N, Noridomi K, Dellisanti CD, Sine SM, et al. Ligand-binding domain of an α7-nicotinic receptor chimera and its complex with agonist. Nat Neurosci. 2011; 14:1253–1259. doi: 10.1038/nn.2908 PMID: 21909087

94. Zouridakis M, Giastas P, Zarkadas E, Chroni-Tzartou D, Bregestovski P, Tzartos SJ. Crystal structures of free and antagonist-bound states of human α9 nicotinic receptor extracellular domain. Nat Struct Mol Biol. 2014; 21:976–980. doi: 10.1038/nsmb.2900 PMID: 25282151

95. Gao F, Bern N, Little A, Wang H, Hansen S, Talley T, et al. Curariform antagonists bind in different orientations to acetylcholine-binding protein. J Biol Chem. 2003; 278(25):23020–23026. doi: 10.1074/jbc.M301151200 PMID: 12682067

96. Yan D, White M. Spatial orientation of the antagonist granisetron in the ligand-binding site of the 5-HT3 receptor. Mol Pharmacol. 2005; 68(2):365–71. PMID: 15914697

97. Zhu F, Hummer G. Pore opening and closing of a pentameric ligand-gated ion channel. Proc Natl Acad Sci USA. 2010; 107:19814–19819. doi: 10.1073/pnas.1009313107 PMID: 21041674

98. Howard R, Murail S, Ondricek K, Corringer P, Lindahl E. Structural basis for alcohol modulation of a pentameric ligand-gated ion channel. Proc Natl Acad Sci USA. 2011; 108:12149–12154. doi: 10.1073/pnas.1104480108 PMID: 21730162

99. Schueler-Furman O, Wang C, Bradley P, Misura K, Baker D. Progress in modeling of protein structures and interactions. Science. 2005; 28:638–642. doi: 10.1126/science.1112160

100. Khalili-Araghi F, Jogini V, Yarov-Yarovoy V, Tajkhorshid E, Roux B, Schulten K. Calculation of the gating charge for the Kv1.2 voltage-activated potassium channel. Biophys J. 2010; 98:2189–2198. doi: 10.1016/j.bpj.2010.02.056 PMID: 20483327

101. Lindert S, Meiler J, McCammon J. Iterative Molecular Dynamics-Rosetta Protein Structure Refinement Protocol to Improve Model Quality. Chem Theory Comput. 2013; 9:e3843–3847. doi: 10.1021/ct400260c

102. Raval A, Piana S, Eastwood M, Dror R, Shaw D. Refinement of protein structure homology models via long, all-atom molecular dynamics simulations. Proteins. 2012; 80:2071–2079. PMID: 22513870

103. Shahsavar A, Kastrup JS, Nielsen E, Kristensen JL, Gajhede M, Balle T. Crystal structure of Lymnaea stagnalis AChBP complexed with the potent nAChR antagonist DHβE suggests a unique mode of antagonism. PLoS One. 2012; 7:e40757. doi: 10.1371/journal.pone.0040757 PMID: 22927902
104. Sauguet L, Shahsavar A, Poitevin F, Huon C, Menny A, Nemecz A, et al. Crystal structures of a penta-meric ligand-gated ion channel provide a mechanism for activation. Proc Natl Acad Sci USA. 2014; 111:966-971. doi: 10.1073/pnas.1314997111 PMID: 24367074

105. Kinde M, Chen Q, Lawless M, Mowrey D, Xu J, Saxena S, et al. Conformational Changes Underlying Desensitization of the Pentameric Ligand-Gated Ion Channel ELIC. Structure. 2015; 23:995-1004. doi: 10.1016/j.str.2015.03.017 PMID: 25960405

106. Maragliano L, Fischer A, Vanden-Eijnden E, Ciccotti G. String method in collective variables: Minimum free energy paths and isocommittor surfaces. J Chem Phys. 2006; 125:024106-14. doi: 10.1063/1.2212942

107. Bolhuis P, Chandler D, Dellago C, Geissler P. TRANSITION PATH SAMPLING: Throwing Ropes Over Rough Mountain Passes, in the Dark. Ann Rev Phys Chem. 2002; 53:291-318. doi: 10.1146/annurev.physchem.53.082301.113146

108. Czerminski R, Elber R. Self avoiding walk between 2 fixed end points as a tool to calculate reaction paths in large molecular systems. Int J Quantum Chemistry. 1990; 38:167-186. doi: 10.1002/qua.560382419

109. Kabsch W. A solution for the best rotation to relate two sets of vectors. Acta Crystallographica. 1976; 32:922-923. doi: 10.1107/S0567739476001873