The α7 nicotinic acetylcholine receptor: Molecular modelling, electrostatics, and energetics

Abstract

The structure of a homopentameric α7 nicotinic acetylcholine receptor is modelled by combining structural information from two sources: the X-ray structure of a water soluble acetylcholine binding protein from Lymnea stagnalis, and the electron microscopy derived structure of the transmembrane domain of the Torpedo nicotinic receptor. The α7 nicotinic receptor model is generated by simultaneously optimising: (i) chain connectivity, (ii) avoidance of stereochemically unfavourable contacts, and (iii) contact between the β1–β2 and M2–M3 loops that have been suggested to be involved in transmission of conformational change between the extracellular and transmembrane domains. A Gaussian network model was used to predict patterns of residue mobility in the α7 model. The results of these calculations suggested a flexibility gradient along the transmembrane domain, with the extracellular end of the domain more flexible that the intracellular end. Poisson-Boltzmann (PB) energy calculations and atomistic (molecular dynamics) simulations were used to estimate the free energy profile of a Na⁺ ion as a function of position along the axis of the pore-lining M2 helix bundle of the transmembrane domain. Both types of calculation suggested a significant energy barrier to exist in the centre of the (closed) pore, consistent with a ‘hydrophobic gating’ model. Estimations of the PB energy profile as a function of ionic strength suggest a role of the extracellular domain in determining the cation selectivity of the α7 nicotinic receptor. These studies illustrate how molecular models of members of the nicotinic receptor superfamily of channels may be used to study structure-function relationships.

• Nicotinic receptor
• molecular modelling
• electrostatics
• Gaussian network model

Introduction

Ligand gated ion channels (LGICs) are integral membrane proteins responsible for fast neurotransmission at both neuromuscular and neuronal synapses. They also play important roles at presynaptic or extrasynaptic loci, where they may modulate transmitter release or neuronal excitability. Underlying all of the roles, the binding of a ligand (i.e., a neurotransmitter) to the extracellular receptor a fast (∼1 ms timescale) conformational change occurs that opens a transmembrane ion pore, thus allowing the passage of ions. There are a number of major classes of LGICs, including the nicotinic acetylcholine receptor (nAChR) and its homologues (e.g. glycine receptors, GABA_A and GABA_C receptors, and 5-HT₃ receptors (Ashcroft [2000]; Hille [2001]; Karlin [2002]; Lester et al. [2004]; Reeves & Lummis [2002]), and glutamate receptors (Mayer & Armstrong [2004]). The nAChR and its homologues are sometimes known as “Cys-loop” receptors due to a pair of conserved cysteines in the ligand-binding domain (LBD) (Kash et al. [2003]; Schofield et al. [2003]). The nAChR and 5-HT₃ receptors have cation-selective pores, while the GABA and glycine receptor channels are anion-selective.

The nAChR is made up of five subunits, arranged with fivefold or pseudo-fivefold rotational symmetry about a central pore. The subunits can either be homologous (in heteropentameric nAChR, e.g. α₂β?δ in the Torpedo nAChR, or α₂β₃ in many neuronal nAChR) or identical (in homopentameric nAChR such as α7–see below). Each subunit is made up of a large extracellular (EC) domain, a transmembrane (TM) domain, and an intracellular (IC) domain. Subunits of the α type are essential for ligand-binding and are identified by a conserved pair of vicinal cysteines in loop C of the EC domain (Karlin [2002]). The TM domain of each subunit contains four membrane-spanning α-helices (M1 to M4). The pore is lined by a bundle of five (i.e., one from each subunit) M2 helices. The M2 helices play an important role in governing the ion permeation and gating properties of the pore (Changeux et al. [1992]; Corringer et al. [2000]; Karlin & Akabas [1995]; Lester [1992]; Oiki et al. [1988]; Opella et al. [1999]).

A high resolution X-ray structure is not available for a nAChR. However, there are two structures which enable us to ‘reconstruct’ much of the structure of a nAChR. The work of Unwin and colleagues (Miyazawa et al. [2003]; Unwin [1993], [1995], [2000]; Unwin et al. [2002]) using cryo-electron microscopy (EM) has resulted in a 4 Å resolution image of the TM domain of the nAChR from the electric ray Torpedo marmorata, and a model has been built into this image (Protein Data Bank (PDB) code 1OED). The EC domain has also been imaged by EM (Unwin et al. [2002]) but model coordinates for this are not as yet deposited. However, high resolution X-ray structures are available for a water-soluble homologue of the EC domain, namely the acetylcholine binding protein (AChBP) from the snail Lymnea stagnalis (Brejc et al. [2001]; Celie et al. [2004]). The AChBP (PDB code 1I9B for the 2.7 Å resolution structure) has 20% to 26% sequence identity to nAChRs. Furthermore, a chimeric protein formed by fusing the AChBP with the TM domain of a nAChR homologue (the 5-HT_3A receptor) has shown to form ligand-gated channels (Bouzat et al. [2004]). Thus we have molecular models of the TM domain of a nAChR based on EM data, and of a homologue of the EC domain. The problem remains of how to reassemble these fragments into a model of a (more or less) complete nAChR (see Figure 1).

Figure 1. Modelling the intact α7 nAChR. The extracellular (A: EC; light grey) domain was homology modelled using the snail AChBP (PDB code 1I9B) as a template; the transmembrane (B: TM; dark grey) domain was homology modelled using the corresponding Torpedo cryo-electron microscopy structure (PDB code 1OED) as a template. These two domains were then combined (see Methods for details) to form a model of the intact receptor (C). The intracellular domain (IC), for which template coordinates are unavailable, is shown schematically as an ellipse, with arrows indicating the sideways entrance/exit of ions via the windows (dotted line ellipses) in the wall of this domain.

Computational methods provide an opportunity to extend the experimental data, and facilitate extrapolation from molecular models of nAChR structure to interpretation of physiological function. Much of the early work in this area focussed on the transmembrane domain and in particular on the pore-lining M2 helix (Law et al. [2000]; Law et al. [2003]; Saiz et al. [2004]; Saiz & Klein [2002]), which had been shown experimentally as an isolated synthetic peptide to form ion channels in lipid bilayers (Montal [1995]; Oiki et al. [1988]; Opella et al. [1999]). Since the determination of the structure of the AChBP, a number of investigators have used this as the basis of modelling and simulation studies of the EC domain of the nAChR and its interactions with various ligands (Le Novere et al. [2002]). From a modelling and simulation perspective, the fivefold symmetry of the homopentameric nAChRs provides considerable advantages. Homopentameric nAChRs may be formed of either α7, α8, or α9 subunits (Ashcroft [2000]). Of these, the chick α7 receptor has been the most intensively studied, both experimentally (Bertrand et al. [1992]; Corringer et al. [2000]; Galzi et al. [1992]) and computationally (Henchman et al. [2003]; Le Novere et al. [2002]). The α7 nAChR is a neuronal receptor and has been implicated in cognitive function, neurological disorders, and nicotine addiction (Karlin [2002]). It is relatively abundant in the central and peripheral nervous systems (Couturier et al. [1990]; Itier & Bertrand [2001]; Marinou & Tzartos [2003]). It shows some functional differences from muscle nAChR, in particular a relatively high permeability to Ca^2 + ions (Fucile [2004]).

Molecular modelling and simulations may complement experimental structural studies by providing further clues as to mechanisms of channel function. There have been a number of simulation studies of the TM domain components of the nAChR. Individual M2 helices, as well as M2 helix bundles, have been simulated in a variety of environments (Law et al. [2000]; Law et al. [2003]; Saiz et al. [2004]; Saiz & Klein [2002]) to explore the dynamics of helix conformation. A recent model of a pentameric M2 helix bundle has been proposed (Kim et al. [2004]). These studies have provided a detailed picture of the nature of the M2 helix bundle in isolation. There have also been a number of modelling studies of the EC domain of nAChR (Le Novere et al. [2002]) and of related receptors (Reeves & Lummis [2002]) based on the structure of the snail AChBP, which have largely focussed on receptor/ligand interactions. An MD simulation of the α7 EC domain homopentamer provided insights into symmetry breaking patterns in subunit dynamics that may be related to channel gating (Henchman et al. [2003]). There have been rather fewer studies of models of the intact nAChR. A preliminary study of the electrostatics contribution to the energetics of cation permeation through open vs. closed models of the acetylcholine receptor has been reported (Corry [2004]). Here we provide a detailed description of a method to generate an optimal model of an intact α7 nAChR, and examine aspects of its predicted dynamic behaviour in addition to exploring the energetics of ion permeation through this closed state model.

Methods

Modelling

The sequence alignment of the Torpedo nAChR α subunit sequence, of the chick α7 nAChR sequence, and of the snail AChBP sequence is shown in Figure 2. To avoid ambiguity, the sequence numbering scheme used in the remainder of the paper corresponds to that of the CHICK sequence in the alignment in Figure 2. The M2 helix corresponds to residues 263–290, which in turn corresponds to 1′ to 28′ using the ‘prime’ numbering scheme for M2 of Lester and others (Lester [1992]; Lester et al. [2004]). Note that residue D289 (alignment) = D27′, at the C-terminus of M2, corresponds to D5* in the ‘star’ numbering scheme used by (Lester et al. [2004]) for the M2–M3 linker. Residue 5* is the one proposed by Unwin and colleagues (Miyazawa et al. [2003]) to interact with a residue at the tip of the β1–β2 loop of the EC domain. Thus in the chick α7 model we would propose an equivalent β1–β2 loop to M2-M3 loop interaction between residues K68 and D289 (alignment numbering).

Figure 2. Alignment of (partial) sequences of the chick α7 nAChR (CHICK) and Torpedo marmorata α subunit nAChR (TORP) with the Lymnea stagnalis acetylcholine binding protein (ACHBP) sequence. The small vertical arrow indicates the boundary between the AChBP EC and nAChR TM templates, and the horizontal bar indicates the extent of the pore-lining M2 helices. The dotted line connects the two residues (K44 and D264 of the model) that formed the basis of the “other criteria” restraint (see Methods for details).

An overview of our modelling procedure is given in Figure 3A. Briefly, homology models of the α7 TM and EC domains were constructed independently of one another, using the Torpedo TM domain (pdb code 1OED) and snail AChBP (pdb code 1I9B) structures respectively as templates. The resultant models of the α7 TM and α7 EC domains were then ‘docked and joined’ by the procedure described in more detail below. The stereochemistry of the resultant model was then adjusted by energy minimization.

Figure 3. (A) Flow diagram of the procedure used to generate the α7 model. (B) Contour plot of a linear combination of three scoring functions (“chain termini”, “bad contacts”, and “other criteria”) as a function of translation and rotation of the EC domain relative to the TM domain. The scores are normalized from 0 (dark grey, corresponding to a “good” model) to 3.0 (white, corresponding to a “poor” model). The model used in subsequent calculations, i.e., with the optimal rotation and translation, is indicated by the white square.

Homology modelling employed Modeller version 6v2 (Fiser et al. [2000]; Sali & Blundell [1993]). Both the EC domain and the TM domain were modelled as homo-pentamers. For each domain, 100 structures were generated by Modeller and ranked on the basis of energies. The ten best structures from Modeller were analysed in terms of protein stereochemistry using Procheck (Laskowski et al. [1993]) and the structure with the best Procheck score was used in subsequent modelling and calculations.

The first stage of the procedure to model the intact α7 nAChR from the models of the two component domains is to align the fivefold rotational axes of the EC and TM domain models. The parameter space to be scanned for ‘docking’ the two domains is then calculated, in terms of the range of translations along the pore (z) axis and rotations about that axis. A grid search of translations (stepsize 1 Å) and rotations (stepsize 1°) is then undertaken, scoring each resultant model on the basis of three criteria: (i) B, the number of bad contacts between atoms of the EC and TM domains, where a ‘bad contact’ between two atoms is defined as one where the interatomic distance is less than 1.5 Å; (ii) T, the distance between the C-terminus of the EC domain and the N-terminus of the TM domain; and (iii) A, additional criteria, in this case the distance between the Cα atoms of residues K68 and D289 (numbering as in Figure 2) corresponding to the interaction between a Val at the tip of the EC domain β1-β2 loop and a Ser in the M2-M3 linker of Torpedo nAChR that have been suggested to form a contact playing a crucial role in linking conformational transitions in the EC and TM domains (Miyazawa et al. [2003]). The overall scoring function S=w_BB + w_TT + w_AA is given by a linear combination of these individual scores, the individual weights w_B, w_T, w_A having been adjusted to give equal weight to the three individual criteria. A contour plot of the overall score S as a function of translation along and rotation about the z axis is shown in Figure 3B. Note that the ‘best’ model corresponds to that with the lowest value of S. This model was selected and energy minimized using Gromacs (see www.gromacs.org) before being used in subsequent calculations and analyses.

Gaussian network models

Fortran code for Gaussian network model calculations was downloaded from the website of Dr Robert L. Jernigan's laboratory (Iowa State University, http://ribosome.bb.iastate.edu/software.html) and modified locally to allow for larger protein molecules (Atilgan et al. [2001]; Bahar et al. [1997]; Erkip & Erman [2003]; Ming et al. [2002]).

Poisson-Boltzmann calculations

Poisson-Boltzmann (PB) calculations were used to estimate the Born energy of a Na⁺ ion placed at successive points along the pore axis. These calculations were performed as described in Beckstein et al. ([2004]). The program, HOLE (Smart et al. [1996]) yielded the pore radius profile of the nAChR model, and also provided sample points along the pore axis at which to place the ion. The energy minimized model was prepared for the PB calculations: WHATIF (Vriend [1990]) was used to add hydrogen atoms to apolar carbons, and PDB2PQR (Dolinsky et al. [2004]) was then used to assign partial charges and radii for all of the atoms. The net charge on the α7 nAChR thus generated was −15 e. The nAChR structure was then embedded in a low-dielectric slab with thickness 30 Å to simulate the effects of the lipid bilayer.

The PB calculation box contained the slab-embedded nAChR structure (the pore being oriented along the z-axis) and one cation, with dimensions 96×96×210 ų. For each data point, the cation was placed at a sample point on the pore centre line as identified by HOLE. This cation had charge +1 e and radius 1.68 Å, equivalent to the Born radius of sodium (Rashin & Honig [1985]). The program APBS (Baker et al. [2001]) was used to perform the PB profiles. The sample points were each 1 Å apart along the z-axis. The number of grid points was 97×97×193; using more grid points does not significantly alter the shape of the resultant PB energy profile. Fine-grid focusing around the ion was used, in a cube with dimensions 10×10×10 ų. The dielectric constant for the nAChR and the membrane-mimetic slab was 2.0; that for the solvent (water) was 78.5. The radius of a solvent probe sphere was 1.4 Å, the effective radius for a water molecule. The PB potential of mean force (PMF) for inserting the cation at a sample point with coordinate z was calculated as:

Umbrella sampling by molecular dynamics

The free energy profile (or potential of mean force, PMF) of a Na⁺ ion along the pore (z) axis of the α7 M2 helix bundle was calculated by umbrella sampling along the axis combined with the results of an ∼60 ns equilibrium molecular dynamics (MD) simulation of the M2 bundle embedded in a simple model membrane, using methods described in more detail elsewhere (Beckstein & Sansom [2004]; Beckstein et al. [2004]). Ions equivalent to approximately 1 M NaCl were present during these simulations. The MD simulations were carried out with Gromacs 3.2.1 (Lindahl et al. [2001]) with the GROMOS96 force field and the SPC water model (Hermans et al. [1984]).

The simulation system (M2 helix bundle, membrane model, water, ions) consisted of 9626 atoms (1275 protein atoms, 478 membrane atoms, 2589 water molecules, 58 Na⁺ ions and 48 Cl⁻ ions) in a simulation box of dimensions 46×46×80 ų. The simulation was performed using a 2 fs timestep, with Berendsen (Berendsen et al. [1984]) temperature (300 K) and pressure coupling (1 bar) in the z direction only. Long-range electrostatic interactions were treated using PME (Darden et al. [1993]; Essmann et al. [1995]).

Umbrella sampling (Valleau & Torrie [1977]) was performed for 101 windows (spacing 0.406 Å), sampling for 0.6 ns per window. For each window, a frame with positions and velocities from the equilibrium simulation trajectory was selected with a Na⁺ ion close to the window centre. The ion was restrained with a force constant of 14.83 kJ mol⁻¹ Å⁻²; this allows diffusion of the ion into neighbouring windows (the restraint energy is ∼1.5 kT at the centre of the neighbouring windows). Each window simulation was started from the positions and velocities of the corresponding equilibrium trajectory frame. The z-coordinate of the ion was recorded at every time step (2 fs). The WHAM method (Kumar et al. [1992]) was used to unbias the ion distribution. The tolerance in the self-consistency cycle was 10⁻⁶kT. The initial 0.1 ns of each 0.6 ns sample were discarded as equilibration, thus using 0.5 ns per window. The umbrella potential of mean force (PMF) and equilibrium MD derived PMF (the latter derived using W(z) = − kT ln n(z)/n₀, where n(z) is the density at z and n₀ the density in the bulk) were combined by matching the graphs in the channel entrance regions where both methods are accurate and thus provide overlapping results.

Results

Modelling

The modelling procedure allowed us to generate a model consistent with the input restraints, namely chain connectivity and minimal stereochemically unfavourable contacts between the EC and TM domains, plus placing the EC domain β1–β2 loop close to the M2–M3 loop in the TM domain. As can be seen from Figure 3B these criteria yield a single, well defined minimum on the scoring function contour plot, and thus a model can be unambiguously assigned.

An overview of the model is provided in Figure 1C. The EC domain interacts closely with the TM domain. The M4 helices are packed rather more loosely on the surface of the M1–M3 TM helix bundle. Note the absence of a model for the intracellular domain of the protein (formed by the region of polypeptide chain between M3 and M4) which is due to the absence of a template structure for this region. This may have important consequence with respect to future calculations of ion conductance and selectivity (Kelley et al. [2003]).

Having generated an energy-minimized model of the combined EC and TM domains, it is possible to analyse predicted dynamic and energetic properties of this model in relation to nAChR function.

Gaussian network model analysis

The model we have generated is assumed to correspond to a closed conformation of the nAChR (see below). It is therefore of interest to attempt to predict the large scale motions of the protein in order to understand their possible relationship to the conformational transition between a closed and an open state of the receptor. Although in principle this might be achieved by long time scale atomistic molecular dynamics simulations (such methods have been successfully employed for e.g., the EC domain of the α7 nAChR (Henchman et al. [2003])), the absence of the IC domain from the model may result in conformational instability in such simulations in the absence of restraints (Hung & Sansom, unpublished results), and also extended timescales (≫20 ns) would be required to see any significant motions. We have therefore elected to use a more ‘coarse-grained’ approach, which omits atomistic detail but has been shown to allow prediction of overall patterns of mobility in a number of proteins. The method used is a Gaussian network model (GNM) in which each residue of a protein is modelled as a particle (“ball”) connected to spatially adjacent residues by harmonic functions (“springs”). Thus, the greater the number of neighbours a residue has, the less freedom of motion a residue is allowed. This model may be used to predict the equivalent of crystallographic B-values (i.e., residue-by-residue fluctuations). Such models have been shown to give good predictions of experimental B-values for a number of proteins (Atilgan et al. [2001]; Bahar et al. [1997]; Erkip & Erman [2003]; Ming et al. [2002]).

We performed a GNM model calculation for the homo-pentameric α7 model. From visual inspection of the calculated relative mobility mapped onto the structure of the model (Figure 4A) it can be seen that there are number of “hotspots” i.e., regions for which the predicted mobility is higher than for the remainder of the structure. These correspond to two loops in the EC domain, those around residues E184 and C212 (alignment numbering), and the upper (i.e., extracellular) ends of TM helices M1, M3 and especially M4 (see Figure 4B).

Figure 4. (A) Gaussian network model (GNM) results mapped onto the Cα trace of a single subunit, where each residue is coloured according to the calculated mobility (blue = low mobility, red = high mobility). Note that the calculations were on the pentameric assembly, but the results for a single subunit are shown for clarity. The two key residues (K68 of the EC domain and D289 of the TM domain) at the EC/TM domain interface are shown in spacefill format. The broken grey line is indicative of the unmodelled sequence between TM helices M3 and M4. (B) Predicted relative mobility (from the GNM) of residues as a function of residue number in the model; to convert to the alignment numbering scheme subtract 24).

The two high mobility loops in the EC domain correspond to two binding sites, one (residues 213–220) for the neurotransmitter acetylcholine and its analogues (this site is also a high affinity bungarotoxin, BgTx, binding site), and one (residues 183–186) a low affinity BgTx binding site (Celie et al. [2004]; Marinou & Tzartos [2003]; Samson et al. [2001]). Interestingly, the region of the polypeptide linking the EC domain to the TM domain is predicted to be quite flexible. There is a gradient in the predicted mobility such that the EC ends of the TM helices are on average of higher mobility than the IC ends. Also, the outer helices, especially M4, have a higher predicted mobility than the inner ones. This is consistent with preliminary MD simulation studies of the TM domain in a DMPC bilayer (Hung and Sansom, unpublished data) but is somewhat different from the model of a relatively inflexible outer helix bundle (M1, M3, M4) with a flexible inner M2 bundle proposed on the basis of the EM images (Miyazawa et al. [2003]; Unwin [1995]).

One should be cautious in over-interpreting such results in terms of possible gating models for the nAChR. However, the gradient in mobilities of the TM helices is suggestive of a model in which both inter- and intra-helical motions contribute to channel gating. We note that both atomistic simulations (Hung et al. 2004) and F-value analysis of the kinetics of mutant nAChRs (Cymes et al. [2002]) are consistent with a model where the upper half of the M2 helix moves more than the lower (N-terminal) half. It is interesting that the EC-M1 linker is of high predicted mobility whereas the mobility in the vicinity of the suggested contact between the β1–β2 loop of the EC domain and the M2–M3 loop is of lower mobility. The β1–β2 to M2–M3 contact in the α7 model seems to be mediated by a salt bridge (K68 to D289; see Figure 4A). This may be part of the mechanism responsible for transmitting ligand-induced conformational change from the EC to the TM domain.

Pore properties

We have also explored the properties of the pore through the centre of the α7 nAChR model. It must be remembered that this model corresponds to a closed state of the channel: however, the pore properties are still quite informative. The pore radius profile was analysed using the program HOLE (Smart et al. [1996]). The resultant pore-lining surface (Figure 5A) and pore radius profile (Figure 5C) reveal a number of constrictions along the length of the pore, which spans from z ∼−25 Å (at the intracellular end) to z ∼+80 Å at the extracellular mouth. The main constriction of the pore is in the TM domain, where the radius drops to ∼2 Å in the vicinity of the ring of Thr residues at 6′ and to ∼3 Å in the hypothesized location of the gate of the nAChR at Leu 9′ and Val 13′. The latter region thus consists of two conserved rings of hydrophobic residues, on the M2 of each subunit. There is a less pronounced constriction (radius ∼5 Å) in the EC region at z ∼+45 Å. When considering these constrictions it is useful to recall that the radius of a solvated Na⁺ ion, allowing for just a single solvation shell, is ∼4 Å.

Figure 5. (A) Model of the α7 nAChR oriented such that the pore (z) axis coincides with that of the graphs in (B) and (C). The solid grey surface represents the pore lining surface, as calculated using HOLE (Smart et al. [1996]). Note that the graphs in (B) and (C) are matched to the dimensions of (A) such that the minima in the pore radius plot (C) correspond to the constrictions of the pore lining surface in (A), as indicated by the two vertical arrows. (B) Comparison between PB energy profiles (for the intact α7 nAChR model, black broken line; and for just the M2 helix bundle, solid grey line) and MD-derived free energy profiles (PMFs for the α7 M2 helix bundle, solid black lines), calculated in the presence of 1 M NaCl. The PMFs calculated by equilibrium simulations (thin solid black line) and by umbrella sampling (thick solid black line) are shown (see Methods for details). The profiles are aligned with the pore lining residues from the M2 helices (indicated by the vertical dotted lines). Note that a negative z value corresponds to the intracellular mouth of the pore, and a positive value to the extracellular mouth of the pore. The M2 helices that line the pore extend from z ∼ −20 Å to z ∼ +20 Å. (C) Pore radius profile calculated for the intact α7 nAChR model using HOLE.

In addition to the pore radius profile, it is possible to estimate the energetics of a Na⁺ ion as it is translated along the pore axis. For the TM region of the channel, we have estimated the free energy profile vs. z (i.e., the potential of mean force, or PMF) via atomistic MD simulations of the M2 bundle embedded in a membrane-like slab and with explicit water molecules and ions present (see Methods for details). For the wider regions of the channel, this PMF can be estimated with suitable accuracy via an extended (∼60 ns) equilibrium simulation. For the constricted centre of the M2 bundle, a more accurate estimate is provided via umbrella sampling, whereby the ion is constrained to lie at successive positions along the pore axis. It is important to note that in both components of the PMF calculation, whilst the M2 helix backbone atoms are constrained the sidechain atoms and water molecules are free to move. The resultant free energy profile (solid lines in Figure 5B) clearly reveals a barrier in the centre of the TM pore formed by the M2 helix bundle. The barrier extends from the T6′ through the L9′ to the V13′ sidechain rings. The height of the barrier is ∼10 kT, i.e., one order of magnitude greater than the mean thermal energy of an ion, thus providing a significant barrier to permeation and demonstrating that the model of the channel indeed corresponds to a closed state. It is interesting that both the equilibrium and umbrella sampling estimates of the PMF indicate a free energy well in the vicinity of the E20′ sidechain ring, a region that mutation studies have implicated in ion selectivity (Galzi et al. [1992]).

A more coarse-grained calculation of the free energy profile allows us to extend the analysis to the influence of the extended lumen formed by the EC domain. This is by using the Poisson-Boltzmann (PB) equation to estimate the Born energy of a Na⁺ ion along the length of the whole channel (broken line in Figure 5B). Note that in this calculation the protein atoms are not free to move, and water molecules are treated as a continuum dielectric rather than as individual molecules. In order to enable a valid comparison with the atomistic PMF calculations, which were based on simulations in the presence of ∼1 M NaCl, the initial PB energy calculations were also for an ionic strength corresponding to 1 M. Once again there is a pronounced barrier to ion permeation in the constricted region of the M2 helix bundle. The height of the barrier (>25 kT) is considerably greater than that for the atomistic simulations. This difference is less pronounced if the PB energy profile calculation is repeated for just the M2 helix bundle of the α7 model, which results in a barrier height of ∼15 kT. Even for just the M2 helix bundle the PB energy barrier (∼15 kT) is greater than that from the atomistic simulations (∼10 kT). This latter difference may be a consequence of intrinsic overestimation of barrier heights for narrow (<3 Å) channels by PB calculations vs. full atomistic PMFs (Beckstein et al. [2004]). However, the same conclusion may be drawn as for the atomistic PMF calculations, namely that the model corresponds to a closed state of the channel.

Based on the free energy profile for the M2 bundle one may obtain an upper limit (g_MAX) on the conductance of the closed state of the α7 channel, using a method developed for the gramicidin channel (Allen et al. [2004]). This yields an estimate of g_MAX=0.02 pS; this is well below the detection limit of patch clamp measurements of single channel currents, and thus would evidently correspond to a closed state of the channel.

The PB calculation also allows one to probe the energetics of the ion within the vestibule. From Figure 5B it can be seen that there is a broad energy well (depth ca. −2 kT) for Na⁺ in the lower part of the vestibule (from z ∼+10 Å to +40 Å). However, it must be recalled that these calculations are at 1 M ionic strength, and so to examine the possible role of the vestibule in more detail, a range of conditions must be considered.

In Figure 6 we compare PB energy profiles for a Na⁺ ion at three different ionic strengths, corresponding to 0, 0.15 and 1.0 M NaCl, and for a Cl⁻ ion at 0.15 M NaCl. At zero ionic strength, the extracellular vestibule generates a very deep (∼−40 kT) well for a Na⁺ ion. However, the depth of this well is very dependent on the ionic strength (see e.g., Corry [2004]). At a physiological ionic strength (0.15 M) the vestibule generates a significant well for Na⁺ ions and so may be expected to play a role in the overall ion selectivity of the α7 nAChR. To confirm this, we calculated the PB energy profile at 0.15 M of a Cl^−. ion (grey broken line in Figure 6). This reveals a considerable (∼20 kT) barrier to anion permeation presented by the vestibule, in addition to the desolvation barrier in the M2 bundle. Thus, the calculations suggest that the extracellular vestibule may contribute to the ion selectivity of the channel (Unwin [1989]).

Figure 6. PB profiles of a monovalent cation (Na⁺) through the intact α7 nAChR channel at various ionic strengths (zero ionic strength = solid black line; 0.15 M NaCl = broken black line; 1.0 M = solid grey line). The PB profile is also shown for a monovalent anion (Cl⁻) in the presence of 0.15 M NaCl (broken grey line).

Discussion

Although the model of the α7 nAChR is incomplete, lacking the intracellular domain, calculations based on this model provide some valuable insights into the structural basis of nAChR function. The GNM predictions of trends in intrinsic mobility/flexibility of the α7 nAChR are of some interest in the context of trying to understand channel gating. Note that it is uncertain whether the AChBP structure corresponds to the closed, open or desensitized state of the nAChR EC domain (Lester et al. [2004]). As the EM structure of the TM domain was obtained in the absence of agonist (acetylcholine), this corresponds to a closed state structure. Thus, the α7 model corresponds most closely to the receptor in a closed conformation. The GNM results reveal several regions of elevated flexibility. The more flexible EC loops (183–186 and 213–220) correspond to those which interact with ligands. This is suggestive of a model in which ligand binding may lead to a modulation of a conformational equilibrium (Goh et al. [2004]) as a component of the conformational transition leading to channel opening. The second indication from the GNM results is that there is a flexibility gradient along the length of the TM domain, down the pore axis. This suggests that the extracellular half of the transmembrane domain may be move relative to the intracellular half. Subtle analysis of the effects of point mutations on the gating kinetics of the nAChR have suggested a conformational wave that is propagated along the length of the molecule from the EC domain to the TM domain (Cymes et al. [2002]; Grosman et al. [2000]). This proposal is consistent with the TM mobility gradient seen in the GNM results, if the wave is propagated along the TM helices from the EC to the IC end of the domain.

The two approaches to permeation energetics that we have adopted, PB energy calculations and fully atomistic MD based PMF calculations, both reveal a significant energetic barrier in the centre of the M2 helix bundle. This is consistent with a ‘hydrophobic gate’ model, in which a relatively narrow (but not fully occluded) hydrophobic region in the centre of the pore forms the barrier to ion permeation in the closed channel. As shown by model calculations, such a gate could be opened by a modest increase in radius and/or polarity of this region (Beckstein et al. [2003]; Beckstein et al. [2001]; Beckstein & Sansom [2004]). Free energy profile calculations on the TM domain of the Torpedo nAChR suggest a similar hydrophobic gate region in the centre of the pore (Beckstein & Sansom, ms. in preparation). A comparable gating model has been proposed based on simulations of the bacterial mechanosensitive channel MscS (Anishkin & Sukharev [2004]). A hydrophobic gate for the nAChR is consistent with the structural (EM) data (Miyazawa et al. [2003]; Unwin [1993], [1995]) and with the single channel kinetic analysis of mutational data (Cymes et al. [2002]). It has been argued on the basis of accessibility studies of cysteine mutations that the gate is at the intracellular end of the transmembrane pore rather than in the centre (Karlin [2002]; Karlin & Akabas [1995]; Wilson & Karlin [1998]; Wilson et al. [2000]). This does not readily agree with the current calculations. However, we note that application of a similar cysteine accessibility procedure to the 5HT₃ receptor (Panicker et al. [2002]) has yielded results consistent with a gate located midway along M2, consistent with the current study.

A number of other models of intact nAChR have been presented (e.g., Corry [2004]). The current study is the first to: (a) present an explicit and objective protocol for combining the EC and TM domains in one model; and (b) provide a comparison of continuum electrostatic and atomistic (i.e., MD based) free energy analyses of the barrier(s) to cation permeation in a closed state nAChR. As seen in calculations on model nanopores (Beckstein et al. [2004]) it is important to compare the two energetic approaches in order to obtain realistic bounds on the height of the energetic barrier, and hence on the predicted magnitude of the (closed state) single channel conductance.

It is important to assess critically possible limitations of the approaches adopted in this study. Perhaps the main limitation is that the TM domain model is based upon relatively low resolution structural data. However, in the absence of higher resolution structural data there is little one may do to remedy this: it is simply important to be aware of this limitation. In terms of the GNM calculations of residue mobility, it should be remembered that this is a rather ‘broad brush’ approach. It therefore yields results that are suggestive of possible (closed state) motions that may be related to gating, but does not provide a detailed atomic resolution mechanism for channel gating.

The permeation energetics calculations are likely to exhibit some sensitivity of the exact height of the energy barrier to relatively small changes in the (TM domain) model. However, calculations of the sensitivity of the PMF calculations to a fivefold reduction in the polypeptide backbone restraints suggest one should not overestimate the sensitivity of the energetic calculations to small changes in model (Beckstein & Sansom, in preparation). Such changes in TM model coordinates are unlikely to qualitatively change the conclusion concerning hydrophobic gating, especially given the comparable results of PMF calculations on the chick α7 and Torpedo nAChRs.

Future computational studies are likely to benefit from improved structural data on which to base models, especially for the IC domain with its potential effects on ion selectivity (Kelley et al. [2003]). A model including the IC domain would also allow unrestrained MD simulations of the full structure of a nAChR in a lipid bilayer, which at present are rather challenging due to the likely effect of the absence of an IC domain model on structural stability during such simulations. Furthermore, once such a model is obtained, it will be timely to develop an open state model by using computational studies to integrate (low resolution) structural data and indirect (e.g., mutational) data in an objective fashion. Such a model could be used to probe, for example, the high permeability of α7 nAChR to Ca^2 + ions, relative to other nAChR species.

Our thanks to all of our colleagues, especially Andrew Hung, Jeff Campbell, David Sattelle, and Nigel Unwin, for their advice and encouragement. Our thanks to the BBSRC, the EPSRC, and the Wellcome Trust for financial support.

Note added at proof stage

A more complete EM structure of the Torpedo nAChR has appeared subsequent to acceptance of this paper (Unwin, N. (2006) J Mol Biol 346: 967–/ 989). This structure is consistent with the model proposed above.

This paper was first published online on prEview on 21 April 2005.