| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Department of Biochemistry, University of Oxford, Oxford OX1 3QU, United Kingdom
Correspondence: Address reprint requests to Mark S. P. Sansom, Tel.: 44-1865-275371; Fax: 44-1865-275182; E-mail: mark.sansom{at}biop.ox.ac.uk.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
L251 and V255). The bending motions of the M2 helices lead to a degree of dynamic narrowing of the pore in the region of the proposed hydrophobic gate. Calculations of Born energy profiles for various structures along the simulation trajectory suggest that the conformations of the M2 bundle sampled correspond to a closed conformation of the channel. Principal components analyses of each of the M2 helices, and of the five-helix M2 bundle, reveal concerted motions that may be relevant to channel function. Normal mode analyses using the anisotropic network model reveal collective motions similar to those identified by principal components analyses. |
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
; Karlin, 2002; Lummis, 2004). It has a pentameric structure, comprising five subunits arranged around an approximate five-fold axis. Subtypes of nAChR are characterized by subunit composition and are broadly divided into neuronal types (composed of homopentamers or heteropentamers of subunits 1–10 and ß1–9) and muscle types (heteropentamers of subunits , ß, 
or 
, and 
). Numerous mutagenesis and labeling studies of the receptor (Changeux et al., 1992; Corringer et al., 2000; Karlin and Akabas, 1995; Lester, 1992) have established the overall topology of the transmembrane (TM) domain of the protein, identifying M2 as the pore-lining helix, shielded from the surrounding lipid environment by the M1, M3, and M4 lipid-exposed helices. More recent mutagenesis studies have explored the nature of the conformational transitions linking acetylcholine binding to the extracellular ligand binding domain (LBD) to opening of the channel formed by the TM domain (Cymes et al., 2002; Grosman et al., 2000; Lester et al., 2004).
Cryolectron microscopy (EM) of the closed (Unwin, 1993) and open (Unwin, 1995) states of the Torpedo nAChR yielded 9 Å resolution images of the inner M2 helix bundle, enabling modeling and simulation studies of the M2-lined pore in relationship to channel function (Adcock et al., 1998, 2000; Sankararamakrishnan et al., 1996). However, such studies were limited by the resolution of the available structural data. More recently, a higher resolution (4 Å) structure of the TM domain of Torpedo nAChR in the closed state has been obtained (Miyazawa et al., 2003). On the basis of this structure, it has been suggested that short loops between the ß-strands of the LBD interact with the M2-M3 loops of the TM domain to form the main direct contact between the LBD and the pore-lining M2 helices. It is also suggested that rotation of the inner sheets of the LBD after ligand binding is communicated through the -subunit M2-M3 loops, resulting subsequently in same-sense rotation of the M2 helices.
The gate of the nAChR is thought to be formed by a ring or rings of hydrophobic residues close to the center of the M2 helix bundle. A variety of experiments suggest the gate may be close to the L9' (i.e., L251; see Fig. 1 for M2 sequences and numbering) ring (Bertrand et al., 1993; Corringer et al., 2000; Lester, 1992; Lester et al., 2004). Model calculations suggest that such rings of hydrophobic residues can form an effective barrier to ion permeation (Anishkin and Sukharev, 2004; Beckstein et al., 2001, 2003; Beckstein and Sansom, 2003, 2004; Corry, 2004). Rotation of the M2 helices is thought to open the channel by increasing the radius of the pore and/or by increasing its polarity by moving polar side chains into the pore-lining region. Rotation of only two M2 -helices is suggested to be sufficient to disrupt the hydrophobic interactions between the five-helix M2 bundle, causing the concerted collapse of the helices against the outer wall, which results in pore widening by 
0.3 nm (Miyazawa et al., 2003).
|
; Grosman et al., 2000). Computer simulations complement structural studies by elucidating some of the dynamical properties of the protein and providing further clues to the dynamic mechanisms of channel activation. 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, 2003; Saiz and Klein, 2002, 2004) to explore the dynamics of helix conformation, thus providing a complement to functional (Montal, 1995; Oiki et al., 1988) and NMR studies of the M2 helix (Opella et al., 1999). More recently, another model of a pentameric M2 helix bundle has been proposed which appears to be consistent with NMR data on the M2 helix peptide (Kim et al., 2004). These studies have provided a detailed picture of the nature of the M2 helix bundle in isolation.
In the intact receptor, the M2 helices are shielded from direct contact with lipids, and instead are surrounded by the outer (i.e., M1, M3, and M4) helices. The interior of the pore is believed to be largely water filled (Hille, 2001). In this work, we seek to elucidate the influence of the environment provided by the outer helices, especially of M1 and M3, on the structure and dynamics of the pore-lining M2s. The simulation in this work therefore serves as a first approximation to the study of M2 dynamics within the framework of a current proposed gating mechanism (Miyazawa et al., 2003), in which the outer helices act as a relatively stationary scaffold within which the M2 helices move. Additionally, we have performed normal mode analyses (NMA) on the TM domain using a coarse-grained anisotropic network model (ANM) to identify possible low frequency collective motions. Comparison of these motions with those extracted from the molecular dynamics (MD) trajectory allows an estimate of the extent to which the atomistic simulation (which describes protein dynamics on a nanosecond timescale) was able to capture larger scale motions.
|
SIMULATION METHODOLOGY |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
). These calculations suggested that ionizable side chains with nonstandard charge states were restricted to Glu and Asp residues, all located in the outer scaffold helices (i.e., D238, E432, and E436 for chain (A); D244 for chain ß(B); E252 for chain (C); D238 and E432 for chain (D); and E477 for chain (E)). All titratable residues in the M2 helices were predicted to be in their standard charge states at pH = 7. In addition to the TM domain, the simulation cell included 1486 octanes, 15,383 simple point charge waters, 259 Na+ and 266 Cl– ions.
|
). The GROMOS96 (van Gunsteren et al., 1996) force field parameters were employed for the simulations. Temperature and pressure coupling were performed using the scheme described in (Berendsen et al., 1984). A constant pressure of 1 bar and pressure coupling constant = 1.0 ps was applied independently in all directions. Water, octane, and protein were coupled separately to a temperature bath at 300 K with temperature coupling constant = 0.1 ps. Electrostatic interactions were evaluated using the particle-mesh Ewald method (Darden et al., 1993; Essmann et al., 1995) with van der Waals interactions truncated at 1.0 nm. Integration time steps of 2 fs were used. Bond lengths were constrained via the LINCS algorithm (Hess et al., 1997).
Analyses of MD trajectories were performed using the GROMACS suite of programs. Secondary structure analysis employed DSSP (Kabsch and Sander, 1983). Helix-bending motions were analyzed using SWINK (Cordes et al., 2002). Pore radius profiles were determined using HOLE (Smart et al., 1996). Porcupine plots (Tai et al., 2001, 2002) of protein concerted motions were acquired using the DYNAMITE web server (Barrett et al., 2004). Visualization of system geometries and evaluation of protein secondary structure were performed using VMD (Humphrey et al., 1996).
Born energy calculations (see below) were performed using the Adaptive Poisson-Boltzmann Solver program (Baker et al., 2001), for which computational parameters were defined as follows: partial charges and radii for the protein atoms were assigned using the PDB2PQR web server (Dolinsky et al., 2004) based on the parameters of the AMBER force field (Cornell et al., 1995); dielectric constants for water and protein were defined to be 78.5 and 2, respectively, with the system temperature set at 300 K. The ionic strength was that of a 150 mM 1:1 electrolyte, NaCl. The Born radius of the Na+ test ion placed at various points within the pore was defined to be 0.169 nm (Rashin and Honig, 1985). The cell size dimensions around the protein were 10 x 10 x 10 nm3.
NMA were performed within the approximation of the ANM. In this approach (Atilgan et al., 2001; Bahar et al., 1997), each residue of the protein is represented by the corresponding C atom, and interacts only with those other residues residing within a specified cut-off radius, rC. The potential between each interacting residue pair is described by a Hookean function, with the sole parameter being the force constant , which is taken to be identical for all interresidue interactions. In this work, rC was defined to be 1.4 nm, whereas was defined as 4.18 kJ mol–1 nm–2. In all of the NMA performed, the M4 helices (each of which is unattached to the rest of its corresponding subunit) were removed, as their inclusion resulted in large and unphysical fluctuations which are unlikely to be present in the complete, intact receptor.
|
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
) water molecules. Water and octane molecules were equilibrated during a 1-ns MD run with positional restraints (force constant 1000 = kJ mol–1 nm–2) on all non-H protein atoms. Preliminary studies showed that in the absence of any restraints, there was large scale structural drift of the outer helices that disrupted the integrity of the TM domain. In particular, the C root mean-square deviations (RMSDs) of the M4 helices rose to 5.5 nm and of the M1 and M3 helices to 3 nm within 1 ns. This is perhaps not surprising given the absence of the two extramembraneous domains (i.e., the LBD N-terminal to helix M1, and the intracellular domain located between helices M3 and M4) from the current model. Thus, a 10-ns "production" trajectory was generated with positional restraints on the backbone atoms of the M1, M3, and M4 helices. The M2 helices, as well as the M2-M1 and M2-M3 linker loops were free to undergo unrestrained motions.
Conformational drift and residue flexibility
The M2 helices remain reasonably stable during the course of the simulation, undergoing little structural drift as evidenced by the formation of a plateau after 4 ns in the RMSD plot of the C atoms of all unrestrained segments (Fig. 3 A, black line). Decomposition of the total RMSD into contributions from the individual subunits (shaded lines, Fig. 3 A) reveals the approximate order in which the subunits attain structural stability. The RMSDs of chains B to E each exhibit a single stepwise increase between 1 and 3 ns, whereas chain A exhibits a relatively large stepwise increase at 4 ns, which is responsible for the minor "bump" in the total RMSD at the same time. The behavior of the total and subunit-specific RMSD, and in particular their tendency to undergo single, discrete increases during the run, may be explained by analysis of their secondary structure using DSSP with respect to simulation time (data not shown). The H-bonding patterns of the M2 helices are largely retained during the trajectory, apart from some disruption at the C-termini ends of the M2 from the B (ß-subunit) and C (-subunit) chains, with more significant uncoiling at the C-termini of M2 A (), D (), and E (). Additionally, there is an uncoiling of a single helical turn at the C-terminus of the M2 of chain A () at 4 ns, and collapse of this M2-M3 loop against the M1. Thus, the RMSD plots and DSSP analyses both suggest that the M2 bundle has apparently settled into a local energy minimum in conformational space after 4 ns of simulation, and for the purposes of current analyses the trajectory is assumed to be stable between 4 ns and 10 ns.
|
-subunit M2s) exhibit the highest fluctuations of the entire bundle. This may reflect a requirement for their flexibility in transmitting structural changes at the ligand binding sites to the M2s after ligand binding. Furthermore, fluctuations in the M2-M3 loops are shown to be correlated with rotations and bending motions of the M2 helices, which in turn have an impact on the pore dimensions, as discussed below.
We have determined the relative flexibility of the TM and extramembranous domains of the protein via calculation of block averaged mean square-fluctuation (MSF) plots (Faraldo-Gómez et al., 2003, 2004) taken for the last 5 ns of the simulation for the M1-M2, M2-M3, and M2 helical regions of the combined bundle (Fig. 3 C). This analysis demonstrates that the M2-M3 loops exhibit the highest overall flexibility, followed by the M1-M2 loops, and finally the relatively rigid M2 TM helix regions. The plateau of the MSF plots after time windows of 2.5 ns suggests that the motions of the unrestrained regions of the protein are reasonably well converged, and thus the structures sampled within the trajectory during this simulation time lie within a local minimum on the energy surface. Examination of the MSF plots for the individual subunits (data not shown) confirms that the highest fluctuations are for the M2-M3 loops of the -subunits.
M2 helix kink and swivel
Inspection of the conformations of the M2 helices during the trajectory revealed several structural characteristics which differ markedly from that of the initial crystal structure (Fig. 4 A). The overall -helical structure is maintained within the timescale of the trajectory, as evidenced by analysis using DSSP (discussed above). However, all of the helices exhibit a more significant kink (bending) compared to that of the initial EM-acquired structure, with a hinge point near the pore center. This is in qualitative agreement with the results of earlier MD simulations of isolated M2 helices (Law et al., 2000), which revealed a propensity of M2 for helix bending in the same region. It also correlates nicely with the results of 
-value analysis (Cymes et al., 2002), which are indicative of bending and/or swiveling around a residue in the center of M2. Kinked M2 helices have been included in a number of earlier models of the nAChR pore (e.g., Tikhonov and Zhorov, 1998).
|
; Cordes et al., 2002). For each helix, an initial analysis was performed on the last 5 ns of the trajectory to identify the average hinge point. Subsequent analyses on each helix were then performed with prespecified hinge points, taken as the trajectory-averaged hinge residue. The hinge points for all of the helices, with the exception of M2 (chain C), lie near the center of the pore in the vicinity of the hydrophobic residues L251 (L9') and V255 (L13'), identified as the sites of the proposed hydrophobic gate by (Lester et al., 2004) and by O. Beckstein and M. S. P. Sansom (unpublished) respectively. For M2 (chain C) the automatic location of the hinge point, was complicated by the presence of P279 near the extracellular side of the helix bundle. For the purposes of comparison, however, we have examined the distortion characteristics of this helix at the same predefined hinge point as those of the others. Having defined the hinge points, we analyzed the kink and swivel angle conformational sampling of all five helices for the final 5 ns of the simulation. The polar plot (Fig. 4 B) shows that the helices remain in a bent geometry throughout the simulation, but there are differences in the nature of the bending motion for each helix.
M2(A) exhibits the highest kink angle, varying from 15° to 45°, with an average of 30°. Helices M2ß(B), M2(D), and M2(E) have somewhat lower kink angles, with averages of 20°. Disregarding the proline-induced kink, chain M2(C) exhibits the lowest bending, exceeding no more than 15° in general. Examination of the swivel angles reveals that helices M2(A), M2ß(B), and M2(E) bend anisotropically, favoring a swivel angle range of 60°. M2(D) samples a wider range of swivel angles (90°), whereas M2(C) appears to bend isotropically, with no apparent preference for particular regions of swivel space. All of the helices exhibit swivel angle sampling in the same region of the polar plot. The anisotropic sampling of kink-swivel space may be understood by visual inspection of the structure, which reveals that all of the helices except M2(C) bend toward the outer helices, with the C-terminal end toward M3 and the N-terminal end toward M1, whereas the apex of the kink points face toward the pore center. Helix M2(C) apparently remains nearly unkinked throughout the simulation, deviating little from its initial structure. Although the differences in the average kink/swivel conformations of the five helices in the current simulation might be expected to be due to differences in the primary structures of the subunits, it is interesting to note that MD simulations of homopentameric M2 bundles (Law et al., 2003) also revealed nonequal average kink angle distributions for each of the five helices. It is possible that, for both the current and previous simulations, such an asymmetry in helical conformations may be due to insufficient simulation time, and that on (for example) millisecond timescales the kink/swivel conformations for each helix may converge to a single state. Possible implications of asymmetric motions are discussed below. Additionally, for the heteropentameric M2 bundle studied currently, differences in primary structure, especially in the M1 helix (see below), may well impart genuine differences on the M2 helical conformations and motions independent of timescale. In addition to structure, asymmetries also arise for the dynamical properties of each subunit's M2, as discussed later.
We note that the region of helix kinking does not correspond to any known helix-distorting sequence motif. Thus, for the nAChR, there may be several external factors which influence the kink-swivel behavior of the helices. Inspection of the motions of the M2s show that they shift closer toward the M1 and M3 helices in the early stages of the trajectory, possibly due to hydrophobic interactions between the M2s and the outer helices. However, for all subunits, a relatively bulky hydrophobic side chain near the center of M1 (F225 for chains (A) and (D), I for ß231(B), 239(C), and 233(E)) moves and points directly toward the M2 hinge point and may therefore cause, or enhance, the degree of M2 helix bending due to steric repulsion. It is noteworthy that reverse mutagenesis studies of the M1 F and I residues (in which the F and I residues were swapped between the subunits; Spitzmaul et al., 2004) revealed alterations in channel gating kinetics, but had no impact on ligand binding kinetics. Thus, given that the current simulation suggests a role for their side chains in causing M2 helical bending, it is possible that the extent and nature of M2 helix bending may have an impact on gating kinetics, although the mechanistic rationale for this requires further study. Additionally, for chains (A) and (D), there are persistent hydrogen bonds between the hydroxyl group of the S249 (M2, S6') side chain and the backbone carbonyl of L235 (M1) near the intracellular side, which may further help stabilize the M2 in a specific bent geometry. Finally, the motion of the M2-M3 loop is correlated with the degree of helix kinking (see below). Overall, then, there are forces driving the M2s toward the outer helices at the helix ends, and steric repulsion driving them toward the pore near the center. These combined forces create a tension which enhances the M2 bending angles from their initial values. For chain C, however, the tension is alleviated by pronounced bending at the proline residue near the extracellular side, preventing any more significant bending near the hydrophobic gate region. The current simulation therefore suggests an important role for the outer helix M1, and especially of the two residues discussed above, in imparting specific, anisotropic bending conformations upon the M2s, which may in turn influence their role in possible gating mechanisms.
Water and H-bonding to the M2 helices
In addition to hydrophobic interactions and specific residue-to-residue contacts with the M1 and M3 helices, the kinked conformations of M2 may also be stabilized by water-to-helix hydrogen-bonding interactions at or near the hinge point which compete with intrahelical H-bonds. To determine the existence of such hinge-stabilizing H-bonds between water and the protein backbone, we have estimated the water-to-backbone H-bond persistence ratio (RH) for the final 5 ns of the simulation for each residue of all five M2 helices in the vicinity of the hinge point. The persistence ratio is defined as
 |
(A), two H-bond persistence "hotspots" were identified, at L251 (L9') and V259 (V13'). These are both pore-facing residues 1 helical turn on either side of the hinge point. The profile for chain M2ß(B) shows H-bond longevity at ßF262, also 1 helix turn displaced from the hinge point. As previously discussed, chain M2(C) exhibits the lowest average degree of helix bending if the hinge point is defined to lie near the pore center. The comparatively low H-bond longevity is consistent with the small relative lack of helix distortion in this region. However, significant distortion near the extracellular end due to a Pro results in exposure of the backbone carbonyl to the pore waters, exhibiting a hotspot near this region. For chain M2(D), a single hotspot is located at L258 (L16'), 1 turn of the helix after the hinge point. Chain M2(E) exhibits high H-bond persistence at L260 (L9') and A261 (A10'), i.e., one turn of the helix before the hinge point in the M2 sequence. Taken together, these results suggest that there are somewhat stronger, more long-lived H-bonds between water and the M2 backbones often one helical turn before the actual hinge point. Pore waters were found to form transient H-bonds mainly to backbone carbonyl O, with the highest H-bond persistence located at residues most exposed to the solvent environment due to deviations from helix linearity, and which are therefore most susceptible to competitive H-bonding from water. This suggest that water-backbone H-bonds may play a role in stabilizing the kinked conformation of the M2 helices (as has been seen in water soluble proteins; Barlow and Thornton, 1988) although other factors likely have a greater influence on the conformation of the M2s.
|
M2s at 10 ns are shown in Fig. 6 A, alongside a diagram of the pore-lining surface, and it may be seen that the hinge points lie adjacent to the major constriction in the center of the pore (at V255, i.e., V13'). The radius profiles of the channel at equidistant points along its axis for both the initial cryo-EM structure and averaged over the final 5 ns of the MD simulation (dashed line) are compared in Fig. 6 B. For the cryoEM structure, two constriction points are situated at S248 (S6') and L251 (L9'), with a third constriction located at V255 (V13'). The constrictions at L251 and V255 correspond to the proposed hydrophobic gate, and free energy profile calculations for a Na+ ion within the M2 pore (O. Beckstein and M. S. P. Sansom, unpublished) indicate a significant energy barrier (8 kT) in this region. For the simulated structure, the constriction at S248 is lost, being replaced by two constriction points at L251 and V255, with the latter showing the lowest average pore radius (0.26 nm). However, as seen in Fig. 6 B, the pore dimensions at these points fluctuate between ±0.05 nm, so that at certain times during the trajectory it is possible for the relative radii of the two constriction points to be equal or even reversed. Thus, the pore radius calculations identify two regions of possible channel gating during the trajectory. The presence of multiple constriction regions may be responsible for possible ambiguities in identifying a unique residue of each subunit which forms the channel gate.
|
; Beckstein and Sansom, 2003, 2004). Thus, to obtain a preliminary quantification of the energetics of trying to move a cation through the (closed) Torpedo nAChR channel, we have computed the free energy of solvation (
Gsolv, or Born energy) of a Na+ ion at predefined positions along the pore axis. This calculation was performed within the framework of the continuum electrostatics model described by the Poisson-Boltzmann equation and implemented in the Adaptive Poisson-Boltzmann Solver code (Baker et al., 2001).
Born energy profiles were compared for four structures: the experimental EM structure (after energy minimization), as well as three structures acquired from the final 5 ns of the MD trajectory. The latter were chosen on the basis of constriction point radius, RCONSTRICT, i.e., the pore radius in the vicinity of the V255 (V13') constriction. Three structures were selected, with RCONSTRICT = 0.26 nm, 0.19 nm, and 0.29 nm. These structures had RCONSTRICT values close to, below, and above the mean for the simulation, respectively.
All of the Born energy profiles (Fig. 7) exhibit significant (i.e., >> kT) peaks in the region of the L251 (9') and V255 (13') rings, effectively constituting energetic barriers to cation permeation, and providing confirmation that the initial experimental as well as MD snapshot structures describe the channel in a closed state. Comparison of the Born energy profiles with the corresponding pore radius profiles (not shown) shows a close correspondence between Born energy barriers and pore constrictions (as might be anticipated). For the structure whose pore profile approximates that of the trajectory average (i.e., the RCONSTRICT = 0.26 nm snapshot), the permeation barrier height is 15 kT, comparable to that of the experimental structure with 18 kT. For the RCONSTRICT = 0.19 nm structure, the barrier height is 28 kT. This suggests that although on average the permeation barrier during an MD simulation is not significantly different from that of the EM structure, the barrier may fluctuate significantly. At the positions of the pore lined by the hydrophobic side chains of residues L251 (9') and V255 (13'), the cation solvation energy increases with decreasing pore radius. These observations are consistent with the expectation that narrower regions of the pore, if lined by apolar side chains, present a more hydrophobic environment in the neighborhood of the cation, increasing the energetic cost of placing it there. Thus this analysis illustrates the principle of hydrophobic gating adopted by this and other ion channels (Beckstein et al., 2001; Beckstein and Sansom, 2004).
|
E262 (20'), likely due to the stabilization of the Na+ by the ring of negatively charged side chains contributed from the - and ß-subunits. This stabilization is more significant for the initial experimental structure, which has a well depth of –5 kT, and is less so for the MD derived structures, with well depths of from –1 to –2 kT. We note that the pore radii at the 20' position for the MD structures are all greater than the experimental structure, thus placing the side chains further away from the Na+ pathway and reducing the favorable cation-anion interactions. Although in this work the electrostatic contributions from the ligand-binding domain are absent, it has been shown that for the 7 nAChR receptor the stabilization described above for the Torpedo nAChR appears to be present as suggested by an energy well at 20', although the presence of the large extracellular (ligand binding) domain introduces a plateau of negative Born energy on the extracellular side above this position (Amiri et al., 2005). The results discussed above therefore reveal that there are three points along the channel that may influence ion permeation: the L and V residues at the hydrophobic gate region forming an energetic barrier, and (perhaps less important) the negatively charged region forming a somewhat favorable electrostatic environment for cations. This suggests specific requirements for conformational changes that must occur during channel opening that will favor cation permeation, namely, enhancement of the pore width at 9' and 13' (and/or rotation of the M2 helices so that more polar backbone atoms line the pore) and reduction in the pore width at 20' to optimize possible Na+ side-chain interactions. Helical-bending motions of the M2 helices is one way by which the above may be achieved; a lesser bending angle may increase the pore radius near the hydrophobic center, and simultaneously reduce the width at the charged 20' position. It is interesting therefore to note that helix bending is one of the concerted motions identified during the MD simulation (see discussion on principal components analysis below).
We note that the influence of atomic-level interactions between the cation, water, and pore atoms, as well as entropic effects, have been neglected in the current model, averaged out by a continuum representation. These will certainly have an impact on the value of the solvation energy, especially for relatively narrow, hydrophobic regions of a pore. Atomic-level free energy calculations (such as umbrella sampling methods) are needed to obtain more quantitative results. We have performed investigations of the relationship between continuum models and atomic-level free energy calculations for simple models of channels (Beckstein et al., 2004), the results of which suggest that continuum models may be adequate for illuminating the trends in the Born energy as a function of pore width and polarity.
Principal components analysis
We have investigated the possibility of concerted motions within the individual M2 helices using principal components analysis (PCA) (or essential dynamics; Amadei et al., 1993), enabling the visualization of the directions and extents of the principal motions of C atoms within the protein which move in a correlated fashion for a particular eigenvector projected along the trajectory. The first three eigenvectors account for 45% of the motion observed in the last 5 ns of the simulation trajectory, with cosine contents (Hess, 2000) of 8%, 39%, and 3%, respectively. These eigenvectors for the M2(A) helix (along with its attached M1-M2 and M2-M3 loops) are shown in Fig. 8, A–C. All three eigenvectors indicate concerted motions between the extramembranous loops (and the M2-M3 loop in particular) with movements of the helical region atoms. However, the particular helical motions differ between the eigenvectors. The first eigenvector (EV1) shows unidirectional motion of the helical segments, all of which move toward the direction of the outer helices, with some bias toward the extracellular side, corresponding to collapse of M2 against the rigid scaffold. EV2 exhibits intrahelical movements, with downwards motions for the upper half of the helix, and upwards motions for the lower half, resembling a helix-bending motion, as identified by SWINK analysis (discussed above). Finally, EV3 suggests rotational motion of the entire helix. Similar analyses for the other M2 helices reveal qualitatively similar results (not shown). In sum, the first three EVs indicate that the motions of the M2-M3 linker loop are highly correlated with bending, translational, and rotational motions of the M2 helix relative to the scaffold helices. Since the mechanism of gating is currently proposed to involve such motions, this analysis reinforces the notion that conformational changes at the LBD may be transmitted to the pore via the M2-M3 linker and furthermore suggests that helix bending may also play a role in gating. In particular, helix-bending motions may play an important role in channel gating by altering the height of the permeation energy barrier at the hydrophobic girdle as well as the depth of the energy well near the extracellular end of the TM domain, thus mediating the degree of interaction between cations and the pore.
|
-subunits and the -subunit appear to rotate in the same sense. The ß M2 also undergoes some degree of rotation, to a lesser extent. The -subunit, however, mainly undergoes to-and-fro translational motions relative to the outer helices. The helices therefore rotate asymmetrically, as might be expected from the heteropentameric nature of the bundle. Dynamical asymmetry has also been observed in simulations of a homology model of the homopentameric 7 LBD (Henchman et al., 2003) and in simulations of the intracellular ligand-binding domain of inward rectifier (Kir) channels (Haider et al., 2005). These results suggest that asymmetric motions may be inherent in multimeric ion channels regardless of subunit composition. If so, this might suggest that a fully concerted model of conformational change (Galzi and Changeux, 1994) may not apply to ion channel gating, and that instead sequential propagation of conformational change through the constituent subunits (Cymes et al., 2002; Grosman et al., 2000) may be a more appropriate model, as suggested by inter alia (Lester et al., 2004).
The rotational motion apparent in EV1 is reminiscent of the rotation of the M2 helices in the proposed gating model; inspection of the positions of the constriction point V side chains between the two extreme projections of EV1 shows their motion away from the center of the pore, and HOLE calculations of the pore radius shows a change of 0.2 Å between the two structures. However, caution should be exercised in interpreting the results of the PCA. The short timescale of the current simulation is such that an actual and complete channel opening event is unlikely to be observed. Indeed, the pore profile calculation shows that, although there is some fluctuation of the minimum pore radius during the trajectory, the channel is essentially in a closed state throughout. Thus, the PCA results are simply indicative of the intrinsic flexibility of the gate region of the TM domain. Nonetheless, it is in principle feasible that, without the presence of the LBD to hold the M2-M3 loops in place, there would be greater freedom for the M2 helices to fluctuate between closed and open states, and thus it is possible that even on a 10 ns timescale motions which are important for channel gating may be observed.
Normal mode analysis of a coarse-grained model
A limitation of the preceding analysis is the relatively short timescale (10 ns) of atomistic MD simulations compared to that of gating transitions within the (intact) nAChR. We have attempted to estimate the extent to which the atomistic simulations are able to describe low frequency collective motions (which may be relevant to our understanding of channel gating) by using a more coarse-grained model of the protein motions. Accordingly, we have applied NMA within an ANM to identify possible global modes for the motions of both the initial (i.e., from cryo-EM) and final MD (i.e., t = 10 ns) TM structures to obtain qualitative insights into their collective motions as well as to serve as a comparison with results of PCA analysis of the MD simulations discussed above.
Although the use of a coarse-grained model with a uniform interresidue interaction parameter results in the loss of detailed information about smaller scale, local fluctuations (which are usually side-chain specific), it has been shown (Atilgan et al., 2001; Bahar et al., 1997) that slow, large-scale collective motions are mainly dependent on the tertiary structure of the protein and may be adequately described by a simple network model.
Despite these advantages, it is unclear whether, in general, NMA-based methods are entirely suitable for studying motions related to the function of ligand-gated ion channels. There are likely to be significant energy barriers separating the closed and opened states due to mechanical constraints imposed by the extracellular domain, resulting in anharmonicity of the energy surface. However, in the present case, such energy barriers are missing due to the exclusion of the LBD in the model. The energy surface may therefore be anticipated to exhibit a lower degree of anharmonicity compared to that of the full length receptor, and hence motions in the TM domain which contribute to channel function are more likely to be amenable to study by NMA methods. Obviously, barriers to gating are still likely to exist owing to internal interactions within the TM domain, and an actual channel opening motion may not be revealed by NMA. Nevertheless, previous studies on ion channels such as MscL (Valadie et al., 2003) and KcsA (Shen et al., 2002) have shown that NMA is capable of revealing motions relevant for gating despite the lack of channel opening motions. In this work, therefore, the key assumption is that the collective motions of the regions with the lowest deformation energies, in the absence of the LBD mechanical constraints, are relevant for biological function.
We note that for the MD simulation, the M1 and M3 helices were restrained. However for the ANM model no such restraints were applied. Thus the results obtained from NMA are unlikely to correspond exactly with those from the MD, in that the NMA will reveal motions of the outer helices. However, the motions of the M2 helices (as with all parts of the protein) are governed entirely by their connectivities to other regions, namely the M1 and M3 helices, as well as neighboring subunits. These are constant, and are exactly the same whether the outer helices are motionally restrained or not. Thus the character of the motions of the M2 helices relative to a fixed point (e.g., the pore center) should not differ significantly. For the purposes of comparison with the MD/PCA results, we focus exclusively on the motions of the "free" components (i.e., the M2-M3 loops, the M2 helices, and M2-M1 loops).
The lowest frequency normal modes of both the initial structure and an MD-acquired 10-ns structure were investigated. We note that for the MD structure, those obtained from several other timeframes within the last 5 ns of the trajectory showed similar results (not shown). As usual, the first six modes correspond to eigenvectors with zero eigenvalues (i.e., translations and rotations) and are discarded; the lowest frequency normal mode is mode 7. For the 1OED structure, mode 7 (Fig. 9 A) describes mainly translational fluctuations of the M2-M3 loops, together with the upper part of the M2 helices, with respect to the pore center. This is asymmetric in terms of the pentameric bundle, with the M2's motions directly opposite of those of the - and -subunits. The ßM2 moves along a line roughly tangential with respect to the pore. Qualitatively similar motions are identified for the 10-ns MD structure (not shown), however, here chain (A) moves in an opposite direction to those of ß(B) and (E), whereas (C) and (D) show translation tangential to the pore. The nature of the translational motions results in a bending motion for all of the M2 helices, with the hinge point located approximately at the helix center (i.e., near V255; Fig. 9, B and C). This highlights the possibility that the intrinsic flexibility of the M2s can also arise from the nature of the protein topology; there is lower flexibility at the bottom half of the M2s due to greater connectivities with neighboring subunits, whereas the opposite is true for the upper half. Coupled with the freedom of the M2-M3 loops, the extracellular half may therefore undergo greater motions, resulting in helix bending even without hinge-bending motifs (e.g., involving Gly or Pro) near the kink point. Similar motions are described by mode 8 for 1OED, and modes 8, 9, and 10 for the 10 ns MD snapshot. Other normal modes also exhibit motions which resemble those identified by PCA. Mode 9 for the cryo-EM structure suggests symmetric and unidirectional rotation of the M2 helices, whereas asymmetric rotations (involving only three out of five M2 helices) were observed for the MD structure in modes 11 and 12. However, these rotations do not result in significant pore widening. We note that higher frequency modes for both structures reveal relatively localized fluctuations, such as those restricted to single M2-M3 loops and helix ends, and are not considered here
|
1 µs for the channel opening conformational transition (Chakrapani and Auerbach, 2005). Given recent progress in achieving such times for atomistic simulations of peptide folding (Duan and Kollman, 1998; Simmerling et al., 2002), it is reasonable to assume appropriately coarse-grained channel simulations should be able to achieve such a timescale in the near future. |
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
How are these simulation results related to experimental studies of nAChR and related members of the Cys-loop receptor channel family? The kinking of the helices to constrict the pore in the vicinity of the L9' and V13' hydrophobic rings is consistent with a body of data that places the gate in this vicinity (summarized in e.g., Bertrand et al., 1993; Corringer et al., 2000; Lester, 1992; Lester et al., 2004; Panicker et al., 2002). The majority of the nAChR data have been interpreted in terms of a gate at L9' although it may be difficult to be precise to within one turn of the M2 helix in positioning a gate by mutagenesis and labeling studies. The asymmetry of the M2 helix conformational dynamics would seem to be broadly consistent with the "conformational wave" model of nAChR gating (Cymes et al., 2002; Grosman et al., 2000; Mitra et al., 2004) although as discussed above, caution must be exercised in interpreting the implications of the dynamics exhibited in the current simulations due to the discordant timescales of MD simulations and channel gating in reality. We have therefore correlated the results of our (short timescale) atomistic simulations with the outcome of more coarse-grained simulations using a Gaussian network model.
It is useful to place these results in the context of related theoretical and simulation studies of nAChR. There have been a number of modeling (Kim et al., 2004; Sankararamakrishnan and Sansom, 1995) and simulation studies of pores formed by just the M2 helix bundle (Law et al., 2003; Saiz and Klein, 2002; Saiz et al., 2004), inspired by structural and functional data on channels formed by an M2 helix peptide (Montal et al., 1993; Montal, 1995; Oiki et al., 1988; Opella et al., 1999). However, the EM structure, despite the limitations of resolution, indicates that the packing of the M2 helices is modified by the presence of the outer (M1, M3, M4) helix bundle. The structure of the intact TM domain supports the model of a hydrophobic gate in the center of the M2 helix bundle. The feasibility of such a hydrophobic gating model is supported by MD simulations of water (Beckstein et al., 2001; Beckstein and Sansom, 2003) and of ions (Beckstein and Sansom, 2004) in simplified models of ion channels and by simulations of the putative hydrophobic gate of the MscS mechanosenstive channel (Anishkin and Sukharev, 2004). Recent continuum electrostatics calculations (Corry, 2004) also support such a gating model for the nAChR although the results of such calculations should be treated with some caution in light of comparisons of continuum electrostatics and atomistic PMF calculations of barriers to ion permeation in simple model pores (Beckstein et al., 2004). Thus, the results of the current simulations are consistent with and extend the emergent theoretical model of hydrophobic gating in the nAChR.
The indications of asymmetric conformational dynamics in the M2 helix bundle are significant, especially as comparable asymmetries have been observed in MD simulations of the extracellular ligand-binding domain of the homopentameric 7 nAChR (Henchman et al., 2003). Furthermore, recent kinetic analysis of mutants in the M4 helices of the nAChR (Mitra et al., 2004) have been interpreted in terms of movement of the -subunits before the - and ß-subunits (in mouse nAChR the -subunit replaces the -subunit of the nAChR). Of course, the timescale of the current MD simulations (10 ns) falls several orders of magnitude short of the timescale of channel activation in response to acetylcholine binding (1 ms). However, it may be that the short timescale intrinsic flexibility of the M2 helices reveals at least some aspects of the dynamics of the gate which are modulated within the intact receptor-channel protein by coupling to the wave of conformational change propagated down the protein (Cymes et al., 2002; Grosman et al., 2000) after binding of the agonist to the receptor (i.e., gatekeeper) domain. It is of interest that M2 bending is seen in the coarse-grained calculations, suggesting that in part this may reflect the environment which the remainder of the TM domain presents to these helices.
It is important to consider two limitations of this simulation study. One is the absence of the ligand-binding domain and intracellular vestibule from the model, necessitating the use of restraints to maintain the integrity of the outer (M1, M3, M4) helix bundle. The other is the use of an octane slab to approximate the membrane environment, rather than a lipid bilayer. In particular, the lipid bilayer headgroups will influence water ordering at the membrane-water interface.
To address this, we have recently completed unrestrained MD simulations of the TM domain in a lipid bilayer (A. Hung and M. S. P. Sansom, unpublished), which yield several results that are in qualitative agreement with those discussed in this study. In particular, these latter studies indicate that the minimum pore radius is significantly reduced compared to that of the initial EM structure and that several M2 helices exhibit a significantly kinked conformation, with the hinge point residues coinciding with those described in this study. We are therefore reasonably confident that the behavior of M2 described above is unlikely to be a simulation artifact due to the use of position restraints. However, more detailed analysis and comparison of these and of other recent simulations (Xu et al., 2005) will be needed to be clear as to the possible influence of protein-lipid interactions.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Submitted on September 12, 2004; accepted for publication February 8, 2005.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONSIMULATION METHODOLOGYRESULTS AND DISCUSSIONCONCLUSIONSACKNOWLEDGEMENTSREFERENCES |
|---|
Adcock, C., G. R. Smith, and M. S. P. Sansom. 2000. The nicotinic acetylcholine receptor: from molecular model to single channel conductance. Eur. Biophys. J. 29:29–37.[CrossRef][Medline]
Amadei, A., A. B. M. Linssen, and H. J. C. Berendsen. 1993. Essential dynamics of proteins. Proteins Struct. Funct. Genet. 17:412–425.[CrossRef][Medline]
Amiri, S., K. Tai, O. Beckstein, P. C. Biggin, and M. S. P. Sansom. 2005. The a7 nicotinic acetylcholine receptor: molecular modelling, electrostatics, and energetics. Mol. Membr. Biol. In press.
Anishkin, A., and S. Sukharev. 2004. Water dynamics and dewetting transitions in the small mechanosensitive channel MscS. Biophys. J. 86:2883–2895.
Atilgan, A. R., S. R. Durell, R. L. Jernigan, M. C. Demirel, O. Keskin, and I. Bahar. 2001. Anisotropy of fluctuation dynamics of proteins with an elastic network model. Biophys. J. 80:505–515.
