Side-chain conformation at the selectivity filter shapes the permeation free-energy landscape of an ion channel
Edited by Richard W. Aldrich, The University of Texas at Austin, Austin, TX, and approved June 27, 2014 (received for review May 14, 2014)
Significance
Our mutagenesis work on the ring of glutamates in the charge-selectivity filter region of the muscle nicotinic acetylcholine receptor led us to propose that the rate at which ions permeate through the open channel depends not only on the number of these glutamates, but also on the conformation of their side chains. Because our inferences were made on the basis of electrophysiological observations, however, we decided to test the plausibility of these ideas using computer simulations. Remarkably, the simulations gave ample credence to all aspects of our proposal and allowed us to gain insight into the effect of specific glutamate rotamers on single-channel conductance.
Abstract
On the basis of single-channel currents recorded from the muscle nicotinic acetylcholine receptor (AChR), we have recently hypothesized that the conformation adopted by the glutamate side chains at the first turn of the pore-lining α-helices is a key determinant of the rate of ion permeation. In this paper, we set out to test these ideas within a framework of atomic detail and stereochemical rigor by conducting all-atom molecular dynamics and Brownian dynamics simulations on an extensively validated model of the open-channel muscle AChR. Our simulations provided ample support to the notion that the different rotamers of these glutamates partition into two classes that differ markedly in their ability to catalyze ion conduction, and that the conformations of the four wild-type glutamates are such that two of them “fall” in each rotamer class. Moreover, the simulations allowed us to identify the mm (χ1 ≅ –60°; χ2 ≅ –60°) and tp (χ1 ≅ 180°; χ2 ≅ +60°) rotamers as the likely conduction-catalyzing conformations of the AChR’s selectivity-filter glutamates. More generally, our work shows an example of how experimental benchmarks can guide molecular simulations into providing a type of structural and mechanistic insight that seems otherwise unattainable.
Sign up for PNAS alerts.
Get alerts for new articles, or get an alert when an article is cited.
The role an ion channel can play in the physiology of a cell is dictated by the rate at which ions permeate, the type of ions that permeate, the stimulus that gates the channel open, the rates of interconversion among all conductive and nonconductive conformations, and the channel’s level of expression in the membrane. In this paper, we are concerned with the chemical determinants of the rate at which cations permeate through the muscle nicotinic acetylcholine receptor (AChR), an archetypal neurotransmitter-gated ion channel.
Several “rings” of negatively charged residues decorate the walls of the ion-permeation pathway of the muscle nicotinic AChR. Of these, it is the ring of four glutamates and one glutamine in the first (N-terminal) turn of the pore-lining M2 transmembrane α-helices (the “intermediate ring of charge” at position –1′; Fig. 1A) that lowers the energetic barrier to cation permeation the most (1). Recently, on the basis of single-channel currents recorded from mutant AChRs, we proposed that only two of the four glutamates in the ring contribute to set the size of the unitary currents, and that these glutamates are deprotonated even at pH 6.0 (2). This is in stark contrast with the situation in voltage-dependent Ca2+ channels (CaV channels) and cyclic-nucleotide-gated channels (CNG channels), for example, where all four selectivity-filter glutamates have been suggested to contribute (directly or indirectly) to the formation of one (in CaV channels) or two (in CNG channels) proton-binding sites that are largely protonated at pH 6.0 (3, 4). Moreover, our results led us to propose that the difference between the muscle-AChR glutamates that catalyze cation permeation and those that do not is the conformation adopted by their side chains in the open channel (2). Furthermore, current recordings from mutant constructs bearing only two glutamates in the ring suggested that the interconversion between these alternate conformers depends on the membrane potential with membrane hyperpolarization (that is, making the inside of the cell more negative) favoring the conformation(s) that catalyzes cation conduction (2). From here, we reasoned that the cation-stabilizing conformation(s) must place the (negatively charged) carboxylate group farther (in terms of “electrical distance”) from the cytoplasmic end of the channel than does the nonstabilizing conformation; this led us to refer to the former conformation(s) as “up” and to the latter as “down,” merely for simplicity. Structural details aside, these results supported the notion that the ability of the glutamates at position –1′ to affect the channel’s conductance to different extents has to do with their conformational flexibility.
Fig. 1.
Essentially, nothing is known about the conformation of these glutamates in muscle or nonmuscle AChRs from the application of direct structural approaches, and the insight provided by the structural models of their bacterial (cation-selective) counterparts seems to be of little value. Indeed, all currently available structural models of ELIC, a channel from Erwinia chrysanthemi, correspond to a nonconductive state (for example, refs. 5–7), and in GLIC—a bacterial ortholog from Gloeobacter violaceus for which an open-channel model has been proposed—the glutamates are shifted to an adjacent position of the primary sequence (position –2′ rather than position –1′; refs. 8, 9), and thus are expected to occupy a different location in the M2 α-helix. On the other hand, the structural model of the AChR from Torpedo (a muscle-type AChR) in the open-channel conformation [Protein Data Bank (PDB) ID code 4AQ9; ref. 10] was generated from data analyzed at too low a resolution for the side chains to be located, let alone for their conformations to be determined. Moreover, the single-channel conductance of the homomeric serotonin type-3A receptor (5-HT3AR, another cation-selective member of the pentameric ligand-gated ion-channel superfamily), for example, is so small (∼0.5 pS; ref. 11) that the applicability of structural models of this channel to interpret the cation-conduction properties of the muscle AChR would not be a foregone conclusion. Finally, the glutamate-gated Cl– channel (GluCl), an invertebrate ortholog for which an open-channel model has also been proposed (PDB ID code 3RHW; ref. 12), does not have glutamates at the intracellular end of M2, a feature that is conserved among all known anion-selective members of the superfamily. As a result, it seems prudent to state that the only insight into the structural organization—and certainly, the protonation state—of the muscle-AChR’s ring of glutamates at position –1′, thus far, comes from our single-channel electrophysiology studies (2).
To test our proposal about the pronounced effect of glutamate side-chain conformation on the AChR’s single-channel conductance within a framework of stereochemical rigor, we set out to perform molecular simulations; no other approach seemed to offer the possibility of studying this phenomenon in a more systematic manner and with the required level of atomic detail. To this end, we applied all-atom molecular dynamics (MD) and Brownian dynamics (BD) simulations to a homology model of the mouse-muscle open-channel AChR’s transmembrane pore that we validated against extensive experimental observations. In remarkable agreement with experiments, the simulations pointed to an asymmetric arrangement of glutamate side-chain conformations as the basis for the markedly different contribution of the various subunits to the single-channel conductance. Furthermore, the simulations allowed us to make predictions as to the specific rotamers adopted by the four glutamates in the wild-type channel.
Results
A Structural Model of the Open-Channel AChR.
The main objective of this study was to test the notion that the conformation adopted by the glutamates of the intermediate ring of charge is a major determinant of the muscle-AChR’s single-channel conductance. To address this question using molecular simulations, we needed a structural model of the receptor in the open-channel state. Although such a model has recently become available from tubular crystals grown from isolated postsynaptic membranes of the electric organ of Torpedo fish (10), the low resolution of the electron crystallography images (6.2 Å) precludes the assignment of the location and conformation of the side chains. Thus, in an attempt to correct for these uncertainties while still taking advantage of the valuable insight provided by these data (indeed, the tubular crystals contained plasma-membrane embedded AChRs), we generated a homology model (Fig. 1) combining the latter with NMR-based models of the α4 and β2 AChR subunits (15). To optimize the model, we used side chain-conformation prediction algorithms (16) and a variety of all-atom MD simulations (17, 18), as indicated in detail in Materials and Methods. To obtain a model of convenient size for all-atom molecular simulations while still retaining the most critical parts of the protein for ion conduction, we modeled only the M1–M3 stretch of each adult-type mouse-muscle AChR subunit; this heteropentamer consists of two α1, one β1, one δ, and one ε subunits. Because the ε subunit contains an extra amino acid at position –3′ of the M1–M2 linker that is not conserved in the β2 subunit after which ε was modeled, we omitted this residue (Gly-250) from the homology model. All known anion-selective members of the superfamily from vertebrates also contain an additional residue at position –3′ of each subunit, and the existing structural data (12, 19) indicate that this extra amino acid is accommodated within the preceding M1–M2 linker without affecting the location of position –1′, in the first turn of the M2 α-helix, relative to the pore (Fig. 1C). Thus, we considered that omitting this residue from our model’s ε subunit is unlikely to have a drastic effect on our simulations. Note also that the ε subunit contains a glutamine rather than a (negatively charged) glutamate at position –1′.
Fig. 1C shows the distances from Cα atoms to the central axis of the channel for the pore-lining M2 α-helices of our mouse model, the model of Torpedo’s muscle-type AChR (PDB ID code 4AQ9; ref. 10) and the model of the GluCl α-subunit homopentamer from Caenorhabditis elegans (PDB ID code 3RHW; ref. 12). For the mouse and Torpedo AChR models, the plotted distances are nearly the same, which is not surprising because the latter was used to model the former. Two differences stand out, however, when these distances are compared with those computed from the GluCl structure. On the one hand, the distance corresponding to position –2′ is shorter in GluCl, which (along with the presence of a proline at this position instead of the glycine present in the AChR) contributes to a narrower pore constriction (Fig. 1D). Gratifyingly, this difference is fully consistent with earlier estimations of minimum pore size of the anion-selective (20, 21) and cation-selective members of the superfamily (22, 23) made on the basis of electrophysiological observations. On the other hand, the Cα atom at position 9′ of the AChR is closer to the pore axis than is the Cα atom at position 10′ by ∼1.5 Å (on average), whereas these two Cα atoms are almost equidistant from the central axis in GluCl. That position 9′ in the open-channel conformation of the muscle AChR is more fully exposed to the aqueous pore lumen than is position 10′ is strongly suggested by the experimental finding that a lysine substituted at position 9′ blocks the single-channel current to a larger extent and that a histidine at this position displays a more bulk-like pKa value (24). Thus, we conclude that the noted structural differences between the models of GluCl and the muscle AChR in the open state are amply justified.
The final step of the series of procedures that led to a structural model consisted of ∼170 ns of equilibrated MD simulation. We used the CHARMM36 force field (25) with an applied electrical potential of –100 mV (negative on the intracellular side of the channel) and with harmonic restraints applied on the protein-backbone atoms so as to maintain a stable pore geometry. The goal of this simulation was to let the side chains relax toward their minimum free-energy conformation(s). As the final time points (“frames”) were approached, the side chains of the two α1-subunit (chains A and D) glutamates adopted χ1 and χ2 torsional angles of approximately –65° and 180°, respectively (Fig. 2), which correspond to one of the staggered-angle combinations of dihedral angles for the two sp3–sp3 bonds of the side chain. Hereafter, we use the nomenclature of the “penultimate rotamer library” (ref. 26; see also SI Appendix, Fig. S1), according to which this rotational isomer (rotamer) is designated as “mt” for “minus” (χ1 = –60° ± 45°) and “trans” (χ2 = 180° ± 45°). Similarly, toward the end of the simulation, the torsional angles of the β1-subunit glutamate side chain were χ1 ≅ –70° and χ2 ≅ –55° (that is, an mm conformation), whereas those of the δ-subunit side chain were χ1 ≅ –75° and χ2 ≅ +70° [an mp conformation, where “p” (χ = +60° ± 45°) is for “plus”]. To learn about the plausibility of these side-chain conformations, we decided to delve deeper into the stereochemistry of side chains in the first turn of an α-helix.
Fig. 2.
The glutamates of the AChR’s intermediate ring occur in the first turn of an α-helix (13, 15), a region where four residues in a row—typically referred to as positions Ncap, N1, N2, and N3 (27, 28)—can form only one of the two i→i–4 (donor→acceptor) backbone–backbone hydrogen bonds that are possible at positions in the interior of α-helices. At these four positions, the backbone amide groups cannot form hydrogen bonds with backbone carbonyl oxygens; instead, nearby side chains often act as hydrogen-bond acceptors. Indeed, the fact that AChRs contain an aspartate at position –5′ and a glutamate at position –1′ (that is, the Ncap and N3 positions, respectively) suggests that these side chains can accept each other’s backbone amide group to form reciprocal hydrogen bonds. This is a common motif (known as the “capping box”) that occurs whenever the Ncap and N3 positions are occupied by side chains that are good hydrogen-bond acceptors (27, 29–35). It follows, then, that the torsional free-energy landscape of the N3 glutamate side chains may well be biased toward conformations that allow the formation of these backbone–side-chain hydrogen bonds with the Ncap aspartates.
To learn about the conformation of N3 glutamates in capping boxes of known structure, we surveyed the PDB for structures solved to 2.0 Å or better resolution containing α-helices [identified using STRIDE (36), as implemented in Visual Molecular Dynamics (VMD; ref. 37)] with glutamate at position N3 and threonine, serine, asparagine, or aspartate at position Ncap (these are the four most common amino acids at position Ncap of the capping-box motif); none of these structures corresponded to a membrane protein, however, let alone to an ion channel. Analysis of 200 unique such α-helices (Fig. 3 and SI Appendix, Table S1) revealed that, in most cases (66.5%), glutamate side chains at position N3 adopt an mt conformation (Fig. 4), a rotamer where the side chain is positioned in such a way that it can accept a hydrogen bond from the backbone-amide nitrogen of the Ncap residue while minimizing steric clashes. This is the rotamer that the α1-subunit side chains tend to adopt toward the end of the ∼170-ns MD simulation (Fig. 2). In fewer α-helices (6.0%), N3 glutamates adopt an mm conformation (Fig. 4). This rotamer precludes the formation of the hydrogen bond between the N3 carboxylate oxygens and the Ncap amide group, but instead, it allows the interaction of the glutamate side chain with the surrounding solvent. This is the type of conformation the β1-subunit side chain adopts toward the end of the simulation (Fig. 2). In this rotamer, and in the particular context of the AChR, the negatively charged carboxylate seems to be optimally positioned to interact with both the water inside the pore and the passing cations, whereas the backbone-amide nitrogen can still donate a hydrogen bond to the Ncap side chain. As for the δ subunit, the glutamate side chain adopts an mp conformation (Fig. 4), a rotamer that was observed in only 1 of the 200 glutamates at position N3 that we analyzed. In this rotamer, the backbone amide and the side-chain carboxylate are positioned in such a way that they can form an intraresidue hydrogen bond, but the near-absence of this type of conformation in the set of surveyed structures suggests that the patterns of hydrogen bonds allowed by the other rotamers are far more stabilizing. It has already been noted that the CHARMM36 all-atom protein force field may overestimate the energy associated with the intraresidue electrostatic interaction between the side chain and the backbone-amide hydrogen of glutamate (25). For comparison, inspection of the dihedral angles adopted during the ∼170-ns MD simulation by the side chains of the naturally occurring glutamates at the other (C-terminal) end of the M2 α-helix—a region where no backbone-amide group is left without a backbone carbonyl-oxygen hydrogen-bond partner—revealed that the mt rotamer dominates the torsional landscape followed by mm; the occupancy of the mp rotamer was negligible (SI Appendix, Fig. S2). In light of these ideas, and as an initial way of dealing with this discrepancy, we decided to dismiss the mp rotamer of the glutamate at position N3 of the δ subunit and considered instead the next most frequent type of conformation adopted by this side chain in our MD simulation: the mm rotamer (Fig. 2C).
Fig. 3.
Fig. 4.
BD-Computed Current–Voltage Properties.
We set out to compute the single-channel current–voltage (i–V) relationships of structural models using BD simulations. To this end, several different frames occurring toward the end of the MD simulation (Fig. 2) and containing both α1-subunit glutamates in an mt conformation and both the β1- and δ-subunit glutamates in an mm conformation were selected at random. For one of these time points (Fig. 5A), the BD-estimated conductance at 27 °C with 150-mM KCl bathing both ends of the channel was ∼147 pS (Fig. 5B), which is in good agreement with the experimentally estimated value of ∼138 pS for the wild type, and ∼160 pS for the mutant without the εGly at position –3′, at ∼22 °C, in the cell-attached configuration ([K+]pipette ≅ 147 mM; ref. 2). Importantly, out of a total of 2,443 ion crossings (in both directions) computed to generate this curve, 2,440 corresponded to K+ (and the other three crossings corresponded to Cl–), and an i–V curve computed with a 1:10 ratio of KCl concentrations across the membrane ([KCl]out = 15 mM; [KCl]in = 150 mM) displayed a reversal potential of ∼–62.5 mV (at 27 °C; Fig. 5C), quite close to the value of –59.5 mV (–55.0 mV, if activity coefficients were taken into account) that would be expected from a perfectly cation-selective channel. Furthermore, we found that replacing the glutamates of the α1 subunits and the glutamine of the ε subunit with alanines leaves the BD-computed single-channel conductance largely unaffected (∼138 pS; Fig. 5B), in complete agreement with our experimental observations for this triple mutant (∼142 pS; ref. 2). Similarly, we found that replacing the four glutamates and one glutamine of the intermediate ring with five alanines lowers the BD-computed single-channel conductance to ∼29 pS and introduces some inward rectification, with the rectilinear portion of the i–V plot at negative membrane potentials projecting onto the voltage axis at ∼–52 mV (Fig. 5B). In close quantitative agreement, cell-attached single-channel patch-clamp measurements indicate that the conductance of this all-alanine mutant AChR is ∼31 pS and that the rectilinear portion of the i–V plot at negative membrane potentials projects onto the voltage axis at ∼–40 mV. Hereafter, we refer to this particular time point of the equilibrated MD as “the working model” (Fig. 5A).
Fig. 5.
The Effect of Side-Chain Conformation.
To better understand the effect of different side-chain conformations on the rate of cation conduction, we generated additional structural models, each having the four glutamates at position –1′ in a given rotamer in the background of the working model described above, and computed the corresponding values of single-channel conductance using BD simulations. We limited these computations to the six rotamers of glutamate that are present at least once in our culled library of 200 α-helices (Fig. 3C and SI Appendix, Table S1). As shown in Fig. 6 A and B, we found that the BD-computed single-channel conductance values are (mean ± SE, in pS) 231 ± 3, 193 ± 7, 119 ± 7; 92 ± 8, 78 ± 7, and 53 ± 11 for the models with all four glutamates in the tp, mm, tt, mt, mp, and pt rotamers, respectively. Thus, it seems as though the contributions of the different rotamers to the single-channel conductance partition, roughly, into two groups: tp and mm contribute more (we refer to these conformations as “conduction-catalyzing”), and the other four rotamers—notably including the mp rotamer, which the MD simulations favor at the δ subunit—contribute less. Extending these results to our working model of the wild-type muscle AChR (Fig. 5A), this would mean that the α1-subunit glutamates (which adopt the mt rotamer) contribute to the amplitude of the single-channel currents to a lower extent than do the glutamates of the β1 and δ subunits (which adopt the mm rotamer). To further test this idea, we replaced the α1-subunit glutamates in our model with alanines while leaving the β1- and δ-subunit glutamates unmodified, and computed the corresponding i–V curve; the single-channel conductance turned out to be ∼136 pS, which is very close to the value of ∼147 pS computed by BD for the (four-glutamate) working model. Importantly, the negligible effect of the α1-subunit glutamates on single-channel conductance found through simulations is in striking agreement with our previous experimental findings (2). For the sake of completion, we also computed the i–V curve corresponding to a model where the β1- and δ-subunit glutamates were replaced with alanines while leaving the α1-subunit glutamates unmodified; the single-channel conductance turned out to be ∼61 pS, comparatively close to the value of ∼29 pS computed for the model containing no glutamates at all at this position.
Fig. 6.
The Effect of the Number of Glutamate Side Chains.
To gain insight into the effect of the number of glutamates at position –1′ on single-channel conductance, we again used BD simulations. Starting with a model having five alanines at position –1′, we added glutamates (one at a time) in the mm rotamer. This is the rotamer adopted by the glutamate side chains of the β1 and δ subunits of our working model, and one of the glutamate rotamers that contributes the most to the single-channel conductance according to our BD simulations (Fig. 6 A and B). The computed i–V plots (Fig. 6C) indicate that each added mm glutamate increases the single-channel conductance: in going from zero to one, from one to two, and from two to three, each additional glutamate increases the conductance by 50, 57, and 57 pS, respectively. However, in going from three to four and from four to five mm glutamates, the conductance increases by only 14 and 23 pS, respectively. Experimentally, we found that the single-channel conductance increases by ∼50 pS in going from zero to one (in any subunit), and by ∼60 pS in going from one to two glutamates in the β1 and δ subunits (2). Addition of two more glutamates in the α1 subunits (of course, without control over or knowledge of the conformation adopted by the introduced side chains) did not increase the conductance any further (2). Indeed, this is the finding that previously led us to hypothesize that the four wild-type glutamates make drastically different contributions to the single-channel current amplitude (2). Moreover, whereas the BD-computed conductance of a model with five glutamates in the mm rotamer is 230 pS, that of a model with three mm glutamates (in β1, δ, and ε) and two mt glutamates (in both α1 subunits) is 185 pS. Remarkably, the experimentally estimated conductance of the mutant having five glutamates in the intermediate ring (and without the Gly at position –3′) is also 185 pS (2).
Inspection of the BD-computed ion-permeation potential-of-mean-force (PMF) profiles (Fig. 6D) revealed that a model having an intermediate ring with five alanines poses a broad free-energy barrier to cation conduction that peaks at around position 2′, three residues away from position –1′ toward the extracellular side. We found that each additional glutamate (in the mm rotamer) modeled at position –1′ lowers the PMF in this region in such a way that the peak of the profile eventually shifts from around position 2′ to position 9′; in the wild type (black line), the PMF profile peaks at position 9′.
Overall, the striking similarity between the i–V properties computed by BD simulations and those estimated experimentally from the wild-type channel and several mutants provides ample support to our structural model of the muscle AChR’s pore in the open-channel conformation.
The mp Rotamer of Glutamate.
The assignment of the mm conformation to the δ-subunit glutamate gave rise to a structural model whose BD-computed i–V properties are very similar to those estimated experimentally from the wild-type AChR. Furthermore, with this choice of rotamer, the resulting arrangement of side-chain conformations is such that two of the glutamates are predicted to contribute to the single-channel conductance much more than the other two, in complete agreement with what we have inferred from experiments (2). However, during ∼170 ns of MD simulation, the probability of the glutamate side chain of the δ subunit adopting the mm rotamer was only ∼0.03, whereas the probability of occupying the mp rotamer was ∼0.97; almost no other rotamer was visited (Fig. 2; for comparison, in the same simulation, the β1-subunit glutamate adopted the mm rotamer with a probability of 0.95).
To address the possibility that this discrepancy represents an artifact of the limited duration of the MD simulation, we performed another simulation, this time starting with all four glutamate side chains in the mm rotamer (SI Appendix, Fig. S3). In the first tens of picoseconds, this initially symmetric arrangement of four glutamates in a conduction-catalyzing conformation evolved into a configuration with only one glutamate in the mm rotamer (none adopted the tp), much like during the first simulation (Fig. 2). In the present case, however, it was one of the α1-subunit glutamates that adopted the mm rotamer; the other three glutamates were, largely, mp. We repeated this simulation three more times (reaching a total of 430 ns; SI Appendix, Figs. S4–S6), all starting with four mm glutamates, and in all cases the simulations pointed to a final configuration of essentially only one glutamate adopting a conduction-catalyzing conformation (SI Appendix, Fig. S7). We also ran a 100-ns MD simulation using the AMBER ff99SB* force field (39) for the protein, instead of the CHARMM36 force field that we used in all other cases, and found that this change has essentially no effect on the results. Indeed, only one of the four glutamates adopted a conduction-catalyzing conformation (SI Appendix, Fig. S8).
Inspection of all of the MD simulations presented thus far indicates that the number of glutamate side chains adopting the conduction-catalyzing, mm or tp rotamers is essentially one, which contrasts with the number of two such rotamers predicted by BD and electrophysiological recordings. Furthermore, the MD simulations indicate that the δ-subunit glutamate adopts the rare mp rotamer most of the time. The mp rotamer of glutamate is stabilized by the intraresidue hydrogen bond formed between the carboxylate side chain and the backbone amide, and it has been noted that the importance of this interaction—relative to all other possible interactions that these moieties can establish in the first turn of an α-helix—may be overestimated in MD simulations (25). Thus, we wondered what would happen to the δ-subunit glutamate if the relative stability of the mp rotamer were lowered. To answer this question, we ran four more MD simulations under conditions expected to either weaken this intraresidue bond or strengthen other, competing interactions. First, we removed the negative charge from this glutamate by “mutating” it to glutamine. To preserve the four-glutamate–one glutamine composition of the wild-type ring, we also “mutated” the ε-subunit glutamine to glutamate. We ran a 100-ns MD simulation starting with the same frame as that used for our first simulation (Fig. 2), the only difference being that the chemical identities (but not the conformations) of the δ- and ε-subunit side chains were swapped. Notably, the side chain of glutamine in the δ subunit changed conformation quickly and occupied the mm rotamer with a probability as high as 0.765 after the first 25 ns of simulation (Fig. 7), an increase by a factor of ∼21 relative to that of a glutamate in δ (Fig. 2C). On the other hand, the side chain of glutamate in the ε subunit adopted the mp rotamer nearly all of the time (SI Appendix, Fig. S9). Importantly, in our simulations of the wild-type channel, the side chain of the ε-subunit glutamine rarely dwelled in the mp rotamer.
Fig. 7.
The other three simulations run under conditions aimed at destabilizing the mp rotamer of glutamate led essentially to the same result. Certainly, removing the positive charge from the backbone-amide hydrogen atom (SI Appendix, Fig. S10), increasing the transmembrane potential from –100 mV to –1,000 mV (SI Appendix, Fig. S11), and fixing a K+ inside the pore (between positions –1′ and 2′; SI Appendix, Fig. S12) decreased the occupancy of the mp rotamer and increased the occupancy of the (conduction-catalyzing) mm rotamer, with the occupancies of all other rotamers remaining very low. Thus, the unexpectedly high occupancy of the mp rotamer in some of our simulations may indeed be the result of the challenge posed by the particular pattern of hydrogen bonds in the first turn of an α-helix.
Discussion
Molecular simulations provided us with the chemical and structural rigor that our electrophysiological observations, naturally, lacked. For example, we can now classify the different rotamers of glutamate on the basis of their predicted effect on the amplitude of the single-channel current: according to the BD simulations, the negatively charged side chain of the muscle-AChR’s glutamates of the intermediate ring is best positioned to catalyze the permeation of cations when adopting an mm or a tp conformation. It is tempting to speculate, then, that the two glutamate side chains that contribute to set the size of the muscle-AChR’s unitary current interconvert rapidly (more rapidly than we can detect electrophysiologically) between the mm and tp rotamers. Similarly, the two glutamate side chains that do not contribute appreciably to the currents would interconvert rapidly between rotamers other than mm or tp. Clearly, we cannot rule out the possibility that interconversions between these two classes of rotamer also occur, but if they did, they would have to take place in such a way that at any point in time two of the glutamates adopt one class of rotamer and the other two glutamates adopt the other class of rotamer. Thus, if one of the glutamate side chains crossed rotamer classes, then another one would be expected to do so in the opposite direction quickly thereafter.
Regarding the BD simulations, we find the agreement with experiments to be remarkable, especially considering that this method treats the electrostatics as a continuum and uses static structures for the computation of the i–V curves. We surmise that the fact that cations permeate the channel with their hydration shell nearly intact (as observed during our MD simulations) explains, at least in part, the success of BD in the case of the AChR. However, we would like to point out that the marked numerical agreement between the BD-computed and the experimentally determined single-channel conductance values has its limitations. Indeed, as has been the case for the structural models of other ion channels, our BD simulations failed to reproduce the relationship between single-channel conductance and cation concentration that is observed experimentally for the muscle AChR. Although the electrophysiologically estimated conductance is already at a maximum at 150-mM concentration of monovalent cations (40), the BD-computed values keep increasing even at 1 M KCl (SI Appendix, Fig. S13A). Thus, the numerical agreement may be regarded as somewhat fortuitous; the numbers would not have been so similar if we had chosen to run the experiments and the simulations at a KCl concentration higher than 150 mM. Also, although our BD simulations reproduced the experimentally determined small effect of mutating the ring of aspartates at position –5ʹ to alanines on the single-channel conductance, the simulations greatly overestimated the effect of mutating the acidic side chains at position 20ʹ (SI Appendix, Fig. S13B). Finally, our model lacks the entire N-terminal extracellular domain, the long intracellular M3–M4 linker, the M4 transmembrane segment, and the short C-terminal tail, the first two of which are expected to affect the single-channel conductance of the muscle AChR, at least to some extent (41, 42). Despite these caveats, however, and even if the numerical agreement had only been approximate, the notion that the BD simulations estimated the effect of mutations to the ring of –1ʹ glutamates correctly lends credence to the predictions made by these simulations for situations that cannot be controlled experimentally, such as the effect of defined glutamate rotamers at this particular position of the M2 α-helix.
It may be argued that we could have estimated the number of ion crossings at different voltages using a more detailed computational approach, such as all-atom MD, instead. However, it seems unlikely that a study of the sort we are presenting here—with several replicates for each measurement, a large number of analyzed mutant constructs, and a variety of explored side-chain conformations—could have been completed in a reasonable timeframe using MD even with special-purpose hardware of the kind used, for example, in ref. 43. Note that ion permeation at experimentally accessible voltages (say, between –200 mV and +200 mV) is a relatively slow process, and thus, simulations much longer than those we have performed here to analyze the time course of side-chain dihedral angles would be required to make meaningful comparisons between simulated and experimentally obtained values.
We examined the six most likely rotamers of glutamate in the context of the AChR’s open-channel structural model (Fig. 4) in an attempt to identify what causes the different conformations of glutamate to contribute to the single-channel conductance to different degrees. We analyzed different geometric variables and found that the BD-computed conductance decreases monotonically with the distance between the carboxylate oxygens and the center of mass of the five-leucine ring at position 9ʹ (Fig. 8), that is, the position at which the ion-permeation PMF profile of the wild-type channel peaks (Fig. 6 B and D). On the other hand, no simple relationship exists between the computed conductance and the distance between the carboxylate oxygens and the pore’s central axis (Fig. 8), a distance that disregards the position along the z axis. Thus, it seems as though the terms “up” and “down,” which we tentatively used in previous work to denote the glutamate conformations that make a large (up) and small (down) contribution to the single-channel conductance (2), are actually quite reasonable and more appropriate than, for example, “in” and “out.”
Fig. 8.
We are convinced that as our understanding of ion-channel structure and function becomes increasingly more sophisticated, our appreciation of the functional impact of side-chain conformational flexibility will also increase, and hence, so will our reliance on atomic-resolution theoretical approaches of the sort we used here. Certainly, we speculate that side-chain flexibility may also be functionally relevant in the case of the glutamates in the selectivity filter of CaV channels, CNG channels, and the Orai Ca2+ release-activated Ca2+ channel. In fact, recent MD simulations of the selectivity filter of bacterial voltage-dependent Na+ channels, which contains four glutamates, have pointed to a role of side-chain flexibility in ion permeation (44, 45) although, clearly, it would be reassuring to obtain experimental support for this notion, much as we did for the muscle AChR (2).
Overall, the present study provides a compelling example of how protein function can be shaped profoundly by details that a mere inspection of static structural models may miss. Above everything else, this work highlights the value of single-channel electrophysiology as a powerful tool in structural biology, the fundamental importance of using experimental benchmarks when performing computer simulations, and the pressing need for increasingly more accurate (yet practical) theoretical treatments of protein electrostatics.
Materials and Methods
Generation of a Muscle-AChR Structural Model.
We generated a homology model of the transmembrane pore of the adult mouse-muscle AChR that includes the M1, M2, and M3 transmembrane α-helices, and the two intervening linkers of each subunit (α1-subunit residues: 208–300; β1 subunit: 219–311; δ subunit: 222–314; and ε subunit: 217–309, with residue 250 being omitted; see below). To model the α1 subunit, we used an NMR-based structural model of the α4 AChR subunit (PDB ID code 2LLY; ref. 15), whereas to model the β1, δ, and ε subunits, we used an NMR-based structural model of the β2 AChR subunit (PDB ID code 2LM2; ref. 15). Residue substitutions were introduced to the models of the α4 and β2 subunits so as to match the primary sequences of the mouse α1 (accession number: P04756), β1 (P09690), δ (P02716), and ε (P20782) subunits. Subsequently, these models were structurally aligned to the open-channel model of the muscle-type AChR from Torpedo (PDB ID code 4AQ9; ref. 10) using VMD (37). Because the ε subunit contains an extra amino acid at position –3ʹ (Gly-250) of the M1–M2 linker that is only shared by γ subunits (the fetal counterpart of ε), and because the structure of the ε subunit was modeled on the basis of the structure of a β subunit, we omitted this ε-subunit residue from our homology model. To assign the conformations of the mutated side chains in the initial model, we used SCWRL4 (16), a side chain-conformation optimization software. Using the CHARMM36 force field (25) for lipids, proteins, and ions, this initial model was then placed in a 105 × 105-Å2 membrane of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine using the membrane-builder plugin in VMD (37). The structure was then solvated with TIP3P water and 150 mM KCl for a total of 80,226 atoms (also using VMD); the dimensions of the simulation box were 105 × 105 × 86 Å3. Using NAMD version 2.9 (17) for all all-atom MD simulations, this structure was minimized for 30,000 steps during which only the acyl chains of the phospholipids were allowed to move. This step was followed by 10,000 steps of minimization with no fixed atoms. The structure was then heated to 300 K over 1.5 ns using equally spaced 50-K increments; during the first nanosecond of this simulation, harmonic constraints of 5 kcal·mol–1⋅Å–2 were placed on backbone atoms to allow for side-chain reorientation. To increase the accuracy of the model, the molecular dynamics flexible-fitting algorithm (MDFF; ref. 18) within NAMD (17) was applied using the portion of the electron-density map of the open-channel structural model of the muscle-type AChR from Torpedo (electron microscopy code EMD-2072; ref. 10) that corresponds to the region of the protein being modeled. After 9.5 ns, the structure converged to a stable state whose Cα atoms overlap almost exactly with those of the Torpedo’s open-channel model (Fig. 1C). After MDFF, equilibrium simulations were performed for ∼170 ns using a constant temperature (300 K) and a constant pressure (1 atm). All simulations used a 1-fs timestep with nonbonded interactions calculated every 2 fs and full electrostatics calculated every 4 fs. A 10-Å switching distance and a 12-Å cutoff for nonbonded interactions were used. Long-range electrostatic interactions were calculated using the particle mesh Ewald method with a 1-Å grid spacing and with periodic boundary conditions being applied in all dimensions. A constant temperature was maintained using Langevin dynamics with a damping coefficient of 1 ps–1. All simulations used harmonic restraints of 0.5 kcal·mol–1⋅Å–2 on the protein-backbone atoms so as to maintain a stable pore geometry. The length of the simulation box along the z axis was 86 Å, and hence an electric field of 1.195 × 10−3⋅V⋅Å–1 along the z axis was applied to mimic a transmembrane potential of 100 mV. Unless otherwise stated, a potential of –100 mV (negative on the intracellular side of the channel) was applied during all MD simulations (46).
Umbrella-Sampling MD Simulations.
To compute the ion-permeation PMF along the central axis of the wild-type pore, we performed umbrella-sampling simulations starting with the last time point of the equilibrated MD simulation illustrated in Fig. 2. One hundred and four “umbrella windows” were used with a spacing of 0.5 Å for a total distance of 52 Å. Using the collective-variables module in NAMD, a harmonic restraint of 5 kcal·mol–1⋅Å–2 was placed on an individual ion at every window so as to sample the entire transmembrane region. Each window was simulated for 1.5 ns, and the first 300 ps were removed from the analysis to account for equilibration. The umbrella windows were then combined using the Weighted Histogram Analysis Method (WHAM, as implemented in http://membrane.urmc.rochester.edu/content/wham; ref. 47) to generate the corresponding single-ion PMF profile; the transmembrane potential during these simulations was zero. We also calculated the ion-permeation PMF for the wild-type model using BD simulations (see below) and found that the two profiles are quite similar (SI Appendix, Fig. S14). Thus, because BD-generated PMFs can be obtained at a much lower computational cost, we applied BD for all other PMF profiles computed in this study.
BD Simulations.
The equilibrated protein structure was used as an input for BD simulations using the Grand Canonical Monte Carlo/Brownian dynamics (GCMC/BD) simulation software (48). Input files for the GCMC/BD simulations were generated using the CHARMM Graphical User Interface (49). Simulations were carried out with a grid spacing of 0.5 Å and a simulation box with dimensions of 85 × 82 × 102 Å3 with a 5-Å buffer region. In all simulations, a 150-mM KCl solution bathed both sides of the membrane unless otherwise indicated. The diffusion coefficients of K+ and Cl− in bulk solution were taken as 0.196 Å2/ps and 0.203 Å2/ps, respectively (50). The diffusion coefficients corresponding to these ions inside the channel’s pore were taken as one-half of their values in bulk using a 10-Å switching function on both ends of the channel to mimic the effect of confinement (49). The dielectric constants of both the membrane and protein regions were set to 2, and that of all aqueous regions was set to 80. Both the pore and the membrane were taken to be 30 Å in length. All simulations were run for 1 μs, and each simulation was replicated five times with different velocity seeds at all membrane potentials. BD simulations were also used to calculate multiion PMF profiles using the expression –RTln(Ci/Co), where Co is the charge density of the bulk solution (for 150 mM KCl, Co was taken as 9.03 × 10−5 e·Å–3) and Ci is the average charge density over the entire simulation for each point along the pore’s central axis.
Acknowledgments
We thank K. I. Lee at The University of Kansas for discussions on Brownian dynamics and J. Vermaas at the University of Illinois at Urbana–Champaign for helpful discussions concerning umbrella-sampling simulations. Simulations were performed using supercomputing facilities provided through Extreme Science and Engineering Discovery Environment Grant MCB130181 and the Taub cluster of the Computational Science and Engineering Program at the University of Illinois at Urbana–Champaign. This work was supported by grants from the US National Institutes of Health (Grant R01-NS042169 to C.G. and Grant T32GM008276 to T.J.H.).
Supporting Information
Appendix (PDF)
Supporting Information
- Download
- 2.55 MB
References
1
K Imoto, et al., Rings of negatively charged amino acids determine the acetylcholine receptor channel conductance. Nature 335, 645–648 (1988).
2
GD Cymes, C Grosman, The unanticipated complexity of the selectivity-filter glutamates of nicotinic receptors. Nat Chem Biol 8, 975–981 (2012).
3
MJ Root, R MacKinnon, Two identical noninteracting sites in an ion channel revealed by proton transfer. Science 265, 1852–1856 (1994).
4
X-H Chen, I Bezprozvanny, RW Tsien, Molecular basis of proton block of L-type Ca2+ channels. J Gen Physiol 108, 363–374 (1996).
5
RJ Hilf, R Dutzler, X-ray structure of a prokaryotic pentameric ligand-gated ion channel. Nature 452, 375–379 (2008).
6
G Gonzalez-Gutierrez, et al., Mutations that stabilize the open state of the Erwinia chrisanthemi ligand-gated ion channel fail to change the conformation of the pore domain in crystals. Proc Natl Acad Sci USA 109, 6331–6336 (2012).
7
R Spurny, et al., Pentameric ligand-gated ion channel ELIC is activated by GABA and modulated by benzodiazepines. Proc Natl Acad Sci USA 109, E3028–E3034 (2012).
8
N Bocquet, et al., X-ray structure of a pentameric ligand-gated ion channel in an apparently open conformation. Nature 457, 111–114 (2009).
9
RJ Hilf, R Dutzler, Structure of a potentially open state of a proton-activated pentameric ligand-gated ion channel. Nature 457, 115–118 (2009).
10
N Unwin, Y Fujiyoshi, Gating movement of acetylcholine receptor caught by plunge-freezing. J Mol Biol 422, 617–634 (2012).
11
N Hussy, W Lukas, KA Jones, Functional properties of a cloned 5-hydroxytryptamine ionotropic receptor subunit: Comparison with native mouse receptors. J Physiol 481, 311–323 (1994).
12
RE Hibbs, E Gouaux, Principles of activation and permeation in an anion-selective Cys-loop receptor. Nature 474, 54–60 (2011).
13
N Mnatsakanyan, M Jansen, Experimental determination of the vertical alignment between the second and third transmembrane segments of muscle nicotinic acetylcholine receptors. J Neurochem 125, 843–854 (2013).
14
OS Smart, JG Neduvelil, X Wang, BA Wallace, MS Sansom, HOLE: A program for the analysis of the pore dimensions of ion channel structural models. J Mol Graph 14, 354–360, 376 (1996).
15
V Bondarenko, et al., NMR structures of the transmembrane domains of the α4β2 nAChR. Biochim Biophys Acta 1818, 1261–1268 (2012).
16
GG Krivov, MV Shapovalov, RL Dunbrack, Improved prediction of protein side-chain conformations with SCWRL4. Proteins 77, 778–795 (2009).
17
JC Phillips, et al., Scalable molecular dynamics with NAMD. J Comput Chem 26, 1781–1802 (2005).
18
LG Trabuco, E Villa, K Mitra, J Frank, K Schulten, Flexible fitting of atomic structures into electron microscopy maps using molecular dynamics. Structure 16, 673–683 (2008).
19
PS Miller, AR Aricescu, Crystal structure of a human GABAA receptor. Nature, 2014).
20
J Bormann, OP Hamill, B Sakmann, Mechanism of anion permeation through channels gated by glycine and γ-aminobutyric acid in mouse cultured spinal neurones. J Physiol 385, 243–286 (1987).
21
K Fatima-Shad, PH Barry, Anion permeation in GABA- and glycine-gated channels of mammalian cultured hippocampal neurons. Proc Biol Sci 253, 69–75 (1993).
22
TM Dwyer, DJ Adams, B Hille, The permeability of the endplate channel to organic cations in frog muscle. J Gen Physiol 75, 469–492 (1980).
23
J Yang, Ion permeation through 5-hydroxytryptamine-gated channels in neuroblastoma N18 cells. J Gen Physiol 96, 1177–1198 (1990).
24
GD Cymes, Y Ni, C Grosman, Probing ion-channel pores one proton at a time. Nature 438, 975–980 (2005).
25
RB Best, et al., Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles. J Chem Theory Comput 8, 3257–3273 (2012).
26
SC Lovell, JM Word, JS Richardson, DC Richardson, The penultimate rotamer library. Proteins 40, 389–408 (2000).
27
LG Presta, GD Rose, Helix signals in proteins. Science 240, 1632–1641 (1988).
28
JS Richardson, DC Richardson, Amino acid preferences for specific locations at the ends of α helices. Science 240, 1648–1652 (1988).
29
JA Bell, WJ Becktel, U Sauer, WA Baase, BW Matthews, Dissection of helix capping in T4 lysozyme by structural and thermodynamic analysis of six amino acid substitutions at Thr 59. Biochemistry 31, 3590–3596 (1992).
30
S Dasgupta, JA Bell, Design of helix ends. Amino acid preferences, hydrogen bonding and electrostatic interactions. Int J Pept Protein Res 41, 499–511 (1993).
31
ET Harper, GD Rose, Helix stop signals in proteins and peptides: The capping box. Biochemistry 32, 7605–7609 (1993).
32
JW Seale, R Srinivasan, GD Rose, Sequence determinants of the capping box, a stabilizing motif at the N-termini of α-helices. Protein Sci 3, 1741–1745 (1994).
33
EA Zhukovsky, MG Mulkerrin, LG Presta, Contribution to global protein stabilization of the N-capping box in human growth hormone. Biochemistry 33, 9856–9864 (1994).
34
R Aurora, GD Rose, Helix capping. Protein Sci 7, 21–38 (1998).
35
S Penel, E Hughes, AJ Doig, Side-chain structures in the first turn of the α-helix. J Mol Biol 287, 127–143 (1999).
36
D Frishman, P Argos, Knowledge-based protein secondary structure assignment. Proteins 23, 566–579 (1995).
37
W Humphrey, A Dalke, K Schulten, VMD: Visual molecular dynamics. J Mol Graph 14, 33–38, 27–28 (1996).
38
P Emsley, B Lohkamp, WG Scott, K Cowtan, Features and development of Coot. Acta Crystallogr D Biol Crystallogr 66, 486–501 (2010).
39
RB Best, G Hummer, Optimized molecular dynamics force fields applied to the helix-coil transition of polypeptides. J Phys Chem B 113, 9004–9015 (2009).
40
S Elenes, M Decker, GD Cymes, C Grosman, Decremental response to high-frequency trains of acetylcholine pulses but unaltered fractional Ca2+ currents in a panel of “slow-channel syndrome” nicotinic receptor mutants. J Gen Physiol 133, 151–169 (2009).
41
JA Peters, et al., Novel structural determinants of single-channel conductance in nicotinic acetylcholine and 5-hydroxytryptamine type-3 receptors. Biochem Soc Trans 34, 882–886 (2006).
42
SB Hansen, H-L Wang, P Taylor, SM Sine, An ion selectivity filter in the extracellular domain of Cys-loop receptors reveals determinants for ion conductance. J Biol Chem 283, 36066–36070 (2008).
43
MØ Jensen, V Jogini, MP Eastwood, DE Shaw, Atomic-level simulation of current-voltage relationships in single-file ion channels. J Gen Physiol 141, 619–632 (2013).
44
N Chakrabarti, et al., Catalysis of Na+ permeation in the bacterial sodium channel Na(V)Ab. Proc Natl Acad Sci USA 110, 11331–11336 (2013).
45
C Boiteux, I Vorobyov, TW Allen, Ion conduction and conformational flexibility of a bacterial voltage-gated sodium channel. Proc Natl Acad Sci USA 111, 3454–3459 (2014).
46
B Roux, The membrane potential and its representation by a constant electric field in computer simulations. Biophys J 95, 4205–4216 (2008).
47
S Kumar, JM Rosenberg, D Bouzida, RH Swendsen, PA Kollman, Multidimensional free-energy calculations using the weighted histogram analysis method. J Comput Chem 16, 1339–1350 (1995).
48
W Im, S Seefeld, B Roux, A Grand Canonical Monte Carlo-Brownian dynamics algorithm for simulating ion channels. Biophys J 79, 788–801 (2000).
49
KI Lee, et al., Web interface for Brownian dynamics simulation of ion transport and its applications to beta-barrel pores. J Comput Chem 33, 331–339 (2012).
50
B Hille Elementary Properties of Ions in Solution. Ion Channels of Excitable Membranes (Sinauer, 3rd Ed, Sunderland, MA), pp. 309–345 (2001).
Information & Authors
Information
Published in
Cover image: Pictured is a histological cross-section of an Arabidopsis thaliana root tip, illustrating the overall radial symmetry of tissue layers. The vascular cylinder displays bilateral symmetry with two phloem poles (empty, unstained cells) mirrored across the xylem axis (turquoise cells). Antia Rodriguez-Villalon et al. defined the signaling pathways that determine whether protophloem cells develop fully into plant vascular tissue or fail to differentiate, as shown in the dark blue cell at the 4 o'clock position in the vascular cylinder. See the article by Rodriguez-Villalon et al. on pages 11551–11556. Image courtesy of Bojan Gujas.
Classifications
Submission history
Published online: July 21, 2014
Published in issue: August 5, 2014
Keywords
Acknowledgments
We thank K. I. Lee at The University of Kansas for discussions on Brownian dynamics and J. Vermaas at the University of Illinois at Urbana–Champaign for helpful discussions concerning umbrella-sampling simulations. Simulations were performed using supercomputing facilities provided through Extreme Science and Engineering Discovery Environment Grant MCB130181 and the Taub cluster of the Computational Science and Engineering Program at the University of Illinois at Urbana–Champaign. This work was supported by grants from the US National Institutes of Health (Grant R01-NS042169 to C.G. and Grant T32GM008276 to T.J.H.).
Notes
This article is a PNAS Direct Submission.
Authors
Competing Interests
The authors declare no conflict of interest.
Metrics & Citations
Metrics
Citation statements
Altmetrics
Citations
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
Cited by
Loading...
View Options
View options
PDF format
Download this article as a PDF file
DOWNLOAD PDFLogin options
Check if you have access through your login credentials or your institution to get full access on this article.
Personal login Institutional LoginRecommend to a librarian
Recommend PNAS to a LibrarianPurchase options
Purchase this article to access the full text.