INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Potassium (K) channels are integral membrane proteins (IMPs) that permit K⁺ to cross membranes at near diffusion-limited rates (~10⁷ ions s¹ channel¹), but which are highly selective for, e.g., K⁺ over Na⁺. They form a large and diverse family of membrane proteins and are present in the plasma membrane of almost all cells of a wide range of organisms. For example, voltage-gated K channels of neurons play a central role in the repolarizing phase of action potentials, as well as in controlling electrical excitability (Hille, 1992). K channels are also present in the membranes of nonexcitable mammalian cells (e.g., the  cells of the pancreas; Heron et al., 1998) and in bacterial cell membranes (Schrempf et al., 1995). The x-ray structure of a bacterial K channel (KcsA from Streptomyces lividans) (Doyle et al., 1998) has revealed the pore domain to be formed by a bundle of eight transmembrane (TM) helices, i.e., four M1 helices plus four M2 helices, into which is inserted the tetrameric P-loop bundle, which determines K⁺ ion selectivity. Sequence comparisons and toxin binding experiments (MacKinnon et al., 1998) suggest that the same basic structure is found in the central pore-lining domain of most, if not all, K channels, including voltage-gated (Kv) and inward rectifier (Kir) channels.

Unfortunately, such high-resolution three-dimensional structures are known for only a small number of IMPs, reflecting difficulties of protein expression and crystallization. Thus it would be desirable to be able to predict the structures of IMPs with an acceptable level of accuracy. The TM regions of IMPs appear always to adopt a well-defined secondary structure, in most cases -helical. The two-stage model of membrane protein folding (Popot and Engelman, 1990) proposes that TM -helices are intrinsically stable when in a bilayer-spanning environment. This suggests that a two-stage approach may be adopted in the prediction of the structure of membrane proteins: 1) prediction of the structure of TM helices and 2) modeling the packing together of these helices to form an intact TM bundle. Given that the basic structural motif of a symmetrical bundle of TM helices surrounding a central pore seems to be common to a number of ion channels (Oiki et al., 1990), such an approach to structure prediction seems to merit further attention. However, before attempting prediction it may be important to understand the structures and dynamics of the basic folding domain, i.e., the TM helices.

A number of different algorithms have been developed (von Heijne, 1992; Jones et al., 1994; Persson and Argos, 1994, 1997; Rost et al., 1995, 1996; Cserzo et al., 1997) that predict the number and location of TM helices within the sequence of a membrane protein with a reasonable degree of accuracy. However, for subsequent modeling studies it is important that the exact location within an amino acid sequence of a TM helix is known. Results with a simple viral ion channel (Forrest et al., 1999) suggest that molecular dynamics (MD) simulations may aid in the definition of TM helices. Furthermore, studies of the dynamics of isolated TM helices in a bilayer environment may provide clues to their structural and dynamic properties when they are part of an intact membrane protein (Woolf, 1997, 1998). Furthermore, the pore-lining helices of K channels provide an opportunity for a comparative study of the structure and dynamics of homologous helices from distinct yet related channel proteins.

The bacterial K channel (KcsA) is tetrameric, with each subunit containing two TM -helices (M1 and M2; see Fig. 1). The central pore is formed by a bundle of eight -helices, with the four M1 helices on the outside and the four M2 helices on the inside of the bundle. The P-loop, which lies between the M1 and M2 helices in the amino acid sequence, is inserted into the extracellular mouth of the M1/M2 helix bundle, where it constricts the opening of the pore and forms the selectivity filter of the channel. The M1 and M2 helices of KcsA are homologous to the S5 and S6 helices of Kv channels and are believed to be structural equivalents of the M1 and M2 helices of Kir channels. However, except for the P-loop, KcsA versus Kir sequence homologies are weak, and sequence alignment of these regions may present difficulties.

View larger version (26K): [in this window] [in a new window]   FIGURE 1   Sequences of KcsA-M1 (A), KcsA-M2 (B), Kv-S5 (C), Kv-S6 (D), Kir-M1 (E), Kv-M2 (F). For the Kv and Kir sequences the grey and black shaded regions indicate, respectively, the extent of the "long" and "short" models used in the MD simulations. Note that in the MD simulations the helix termini were blocked by N-terminal acetyl and C-terminal amide (i.e., -CONH2) groups. The secondary structure predictions and, in the case of KcsA helix, the extent of the x-ray structure as indicated by DSSP analysis are given below the sequences. Details of the secondary structure prediction methods are given in the main text.

The primary aim of this study is to analyze the structural dynamics in a full lipid bilayer and water environment of isolated M1 and M2 helices from KcsA and of the equivalent helices from a Kv and a Kir channel (Shaker and Kir6.2, respectively). In particular, we have compared the behavior of the homologous pore-lining helices from the three K channel species to see whether any clues to possible functional roles of their dynamic properties can be revealed.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Secondary structure prediction

Transmembrane helices of the Shaker were predicted using the following programs: TopPred2 (http://www.biokemi.su.se/~server/toppred2; von Heijne, 1992), TMAP (http://www.embl-heidelberg.de/tmap/tmap  info.html; Persson and Argos, 1994, 1997), DAS (http://www.biokemi.su.se/~server/DAS; Cserzo et al., 1997), and PHDhtm (http://www.embl-heidelberg.de/predictprotein; Rost et al., 1995, 1996). The same programs were used to predict the KcsA M1 and M2 helices. Additional programs used for Kir6.2 were memsat (http://globin.bio.warwick.ac.uk/~jones/memsat.html; Jones et al., 1994), TMHMM (http://www.cbs.dtu.dk/services/TMHMM-1.0/; Sonnhammer et al., 1998), and TMPred (http://www.isrec.isb-sib.ch/software/TMPRED  form.html; Hofmann and Stoffel, 1993).

Modeling of TM -helices

Initial models of helices were generated via restrained MD simulations in vacuo, using the simulated annealing protocol described in previous papers (Kerr et al., 1994, 1996). In brief, the model generation consists of two stages. In stage I, the C atoms of the helices are fixed and annealing is started at 1000K, with the gradual introduction of covalent terms and subsequently of the van der Waals potential. In stage II, intrahelical distance restraints are introduced (to maintain an -helical backbone), and electrostatic potentials are gradually incorporated into the potential energy function. This two-stage procedure results in, e.g., 25 different structures for a single helix, differing in side-chain conformations. One of these structures was selected at random and used as the starting point for MD simulations in a bilayer.

Molecular dynamic simulations

The simulation method employed has been described in some detail in previous papers (Tieleman et al., 1999c; Forrest et al., 1999). A fully equilibrated palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer consisting of 128 lipid molecules was used as the starting point for all simulations. A short in vacuo MD simulation was performed in which a radial force was applied to exclude lipid atoms from a cylinder of radius 0.7 nm. A single POPC molecule, the fatty acyl tails of which overlapped with the peptide, was removed. This gave a bilayer with 127 POPC molecules and a central hole large enough to accommodate a TM helix. A single model helix (generated by the restrained in vacuo MD method described above) was introduced into the cylindrical cavity within the bilayer. This system was energy minimized and then solvated with SPC water (using 30-40 water molecules per lipid molecule). For those systems with a net charge on the peptide (assuming all side chains to be in their ionized state), counterions (Na⁺ or Cl) were added by a procedure that systematically positioned the ion instead of each water molecule in turn, eventually selecting the ion position with the lowest potential energy. After addition of the ion(s), the system was once again energy minimized, followed by a short equilibration simulation of 100 ps, during which peptide atoms were restrained to their initial coordinates.

Molecular dynamics simulations were run using GROMACS (Berendsen et al., 1995). A twin-range cutoff was used for longer-range interactions: 1.0 nm for van der Waals interactions and 1.7 nm for electrostatic interactions. A production simulation of 1 ns was run for each helix. The time step was 2 fs, with the LINCS algorithm used to constrain bond lengths. We used NPT conditions (i.e., constant number of particles, pressure, and temperature) in the simulation. A constant pressure of 1 bar independently in all three directions was used, with a coupling constant of _p = 1.0 ps (Berendsen et al., 1984). This allows the bilayer/peptide area to adjust to its optimum value for the force field employed. Water, lipid, and protein were coupled separately to a temperature bath at 300 K, using a coupling constant of _T = 0.1 ps.

The lipid parameters were as in previous MD studies of lipid bilayers (Berger et al., 1997; Marrink et al., 1998). These lipid parameters give good reproduction of the experimental properties of a DPPC bilayer. The water model used was SPC (Hermans et al., 1984; van Gunsteren et al., 1996), which has been shown to behave well in lipid bilayer/water simulations (Tieleman and Berendsen, 1996).

General

Model building was performed using Xplor 3.1 (Brünger, 1992) with the Charmm (Brooks et al., 1983) param19 parameter set. Molecular dynamics simulations were run using GROMACS (http://rugmd0.chem.rug.nl/~gmx/gmx.html) (Berendsen et al., 1995). Essential dynamics and domain motion analysis (using DYNDOM) were performed as described by Hayward and Berendsen (1998). Graphics visualization was done using Quanta v4.1 (Molecular Simulations, Waltham, MA) and Rasmol v2.6 (Roger Sayle, Glaxo-Wellcome). MD simulations were performed on a 10-node SGI Origin 2000, taking ~10 days of cpu time on a single R10000 processor.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Secondary structure predictions

Fig. 1 shows the regions of the Kcsa, Shaker, and Kir6.2 sequences corresponding to the pore-lining helices (M1 and M2 in Kcsa and Kir, and S5 and S6 in Shaker). Alongside the sequences are shown the TM helix predictions from the different programs. It can be seen that although the same core regions are identified by most of the algorithms, there are significant differences in the lengths of the TM helices. In general, the TM helices predicted by TMAP were longer than those predicted by the other three methods. Models of KcsA-M1 (28 residues) and KcsA-M2 (29 residues) were based on the TMAP predictions, as these gave the best agreement with the reported structure of KcsA (note that at the time this study was initiated the coordinates of KcsA were not available and so the exact lengths of M1 and M2 were not known). Two sets of models of Kv-S5 and Kv-S6 were generated. The first ("short") set was based on the minimum length predictions of the core TM helices (DAS with 19 residues for S5S; TopPred2 with 22 residues for S6S). In addition, two longer models of these helices were also modeled by introducing four residues at both the N- and C-terminals. This resulted in models S5L (27 residues) and S6L (30 residues). Note that both of these extended models corresponded quite closely with the TMAP predictions. For Kir6.2 some additional TM prediction methods were employed. In a manner similar to that for the Kv helices, both "short" and "long" helix models were constructed. Kir-M1S and Kir-M2S corresponded to the consensus core of all seven predictions (of 19 and 22 residues, respectively), whereas Kir-M1L and Kir-M2L corresponded to the maximum extent of the TM helix across all of the predictions. The latter gave helix lengths of 26 residues and 29 residues for M1L and M2L, respectively. Note that the longer helix predictions are more consistent with the Kcsa x-ray structure and with the results of a statistical analysis of 77 TM helices in nine membrane proteins of known x-ray structure, which suggest a mean length of 25 (± 5) residues (Ulmschneider and Sansom, unpublished results).

Bilayer models

The TM helices were inserted into the POPC bilayer so as to optimize peptide/bilayer interactions. Current theories of TM helices suggest that charged and amphipathic aromatic residues should be present in the regions corresponding to the water/bilayer interface, i.e., in the region of the phospholipid headgroups. This greatly aided positioning of the helices in the bilayer. So KcsA-M1 has two amphipathic aromatic residues (W4 and Y22) and two charged residues (R4 and E28), while KcsA-M2 has two aromatic residues (W2 and F29) plus a single charged residue (R4). Because the "short" S5S and S6S helices contain neither amphipathic aromatic residues nor charged residues, they were positioned so that they symmetrically spanned the hydrophobic core of the bilayer. In contrast, the "long" S5L helix has three charged residues (R3, E4, and E27), and S6L has an amphipathic aromatic residue (W2) and a single charged residue (K4) at its N-terminus. As with KcsA-M1 and KcsA-M2, these helices were positioned so as to place the charged and amphipathic aromatic residues in the interfacial region. As can be seen from Fig. 2, the short TM helix of S6S, for example, is barely able to span the hydrophobic core of the lipid bilayer, unlike the long S6L helix. In the case of the Kir helices a similar situation arose, although it was somewhat less clear-cut in that the Kir-M1S helix, although too short to fully span a bilayer, did contain aromatic residues close to either end.

View larger version (32K): [in this window] [in a new window]   FIGURE 2   Snapshots of Kv-S6S (A) and Kv-S6L (B) taken at the end of the respective simulations. The helices are shown as ribbons and the P atoms of the lipid headgroups as van der Waals spheres. In B the side chains of the helix are included to illustrate how residues W2 and K4 interact with the lipid headgroups (see main text and Fig. 9). The approximate extents of the bulk water (w), polar lipid headgroup (p), and hydrophobic bilayer core (h) regions are indicated. In C and D the structure of the non-TM segment (see Table 1 and Discussion) is shown at the start and end of a 1-ns simulation. In all four diagrams the N-terminus is at the top of the diagram.

Having positioned the TM helices in the bilayer and run a short (100 or 300 ps) equilibration simulation in which the protein backbone atoms were restrained (see Methods for details), we ran a 1-ns production MD simulation for each helix (Table 1).

Drift during the simulations

The final values of the overall C root mean square deviations (RMSDs) (see Table 1) are all on the order of 0.1-0.2 nm from the starting structures, indicating relatively small deviations over the course of the 1-ns simulations. The largest deviation, 0.22 nm, from the t = 0 structure is for the S6L helix. Interestingly, a somewhat higher average final deviation is seen for the inner (i.e., M2/S6) helices (0.16 ± 0.03 nm) than for the outer (i.e., M1/S5) helices (0.12 ± 0.02 nm). Note that overall deviations of this magnitude are seen in simulations starting from experimental x-ray structures of membrane proteins (Tieleman and Berendsen, 1998; Shrivastava and Sansom, 2000). However, this does not imply that subtle differences in structural dynamics are not observed between the different simulations. For example, visual examination of the structures of S6S and S6L at the end of their simulations (Fig. 2) suggests that the shorter TM segment of S6S undergoes distortions at either end in attempting to span the hydrophobic bilayer core and satisfy its hydrogen-bonding propensities, while the longer S6L model, despite the presence of a central kink, forms an otherwise undistorted -helix. Thus more detailed analysis of the simulation results seems to be warranted.

Secondary structure

Secondary structure as a function of time was analyzed using DSSP (Kabsch and Sander, 1983) (Fig. 3). As anticipated, the terminal residues are more susceptible to loss of -helicity than those residues embedded in the hydrophobic core of the bilayer. This is particularly the case for the "short" TM helix models (Kv-S5S and Kv-S6S, Kir-M1S and Kir-M2S), where the requirement to form H-bonds to the ends of the TM segments leads to their being stretched, with concomitant loss of -helicity. This has been observed in other simulations of TM helices for which both "short" and "long" models were explored (Fischer et al., manuscript submitted for publication). We therefore excluded the "short" simulations from more detailed analysis.

View larger version (74K): [in this window] [in a new window]   FIGURE 3   Secondary structure analysis, using DSSP (Kabsch and Sander, 1983), in all 10 simulations. (A) KcsA-M1. (B) KcsA-M2. (C) Kv-S5S. (D) Kv-S6S. (E) Kv-S5L. (F) Kv-S6L. (G) Kir-M1S. (H) Kir-M2S. (I) Kir-M1L. (J) Kir-M2L. Black squares represent -helical residues, dark gray represents 310-helix, light gray represents turn, and white squares represent random coil.

DSSP analysis reveals that both KcsA-M1 and KcsA-M2 remain in an -helical conformation throughout the simulation. The results are more complex for the (long) Kv helices. For S5L there is transient and limited loss of -helicity around residues S20 and S21. For S6L, the DSSP analysis supports the earlier suggestion (Kerr et al., 1996; Kerr and Sansom, 1996) of two helical regions separated by a region of dynamic flexibility around residue 18. In particular, if one examines the S6L secondary structure trajectory (Fig. 3 F), it can be seen that a switch from -helix to turn conformation occurs at ~300 ps and is sustained throughout the remainder of the simulation. This clearly suggests that S6L behaves like two -helices linked by a nonhelical region. This is also supported by analysis of the time-averaged backbone torsion angles ( and ) for S6L (data not shown). The DSSP results for the Kir helices reveal a similar pattern. In particular, Kir-M1S shows marked loss of helicity at either end, whereas the longer Kir-M1L helix remains stable throughout the simulation. In contrast, Kir-M2S does not show any marked loss of helicity at its ends, probably because it is (just) long enough to span the hydrophobic core of the bilayer. Further analysis will focus on Kir-M1L and Kir-M2L.

The DSSP results may be extended by looking at the time-averaged percentage helicity as a function of residue number for the six "long" helix simulations (Fig. 4). For KcsA-M1 and KcsA-M2 these profiles are relatively flat, although there is some indication of lower helicity in the center of M2 near residue G14. Kv-S5L is -helical for most of its length, apart from a dip in helicity near S20 and S21. In Kv-S6L a very dramatic drop in helicity occurs around residues 17 and 18. Note that this is just one helix turn before the P-V-P sequence motif. This supports the DSSP analysis in suggesting a clearly nonhelical region in the middle of S6L when it is in a lipid bilayer. Interestingly, the PVP motif is highly conserved in Kv channel sequences (it is PIP, a similar motif, in the Kv2 (Shab) family). The profile for Kir-M1L resembles that of KcsA-M1 in that it is more or less flat. In contrast, the Kir-M2L profile suggests some drop in percentage helicity at either end of the TM segment.

View larger version (27K): [in this window] [in a new window]   FIGURE 4   Percentage helicity versus residue number for simulations. (A) KcsA-M1. (B) KcsA-M2. (C) Kv-S5L. (D) Kv-S6L. (E) Kir-M1L. (F) Kir-M2L.

Overall, the picture emerging from the secondary structure analysis is 1) KcsA-M1 is a relatively rigid -helix; 2) KcsA-M2 remains -helical throughout the simulation, but with an indication of flexibility in the vicinity of G14; 3) Kv-S5L is largely -helical, with a flexible C-terminus; 4) Kv-S6L behaves like two -helices joined by a nonhelical linker just before the Pro-Val-Pro sequence motif; and 5) both Kir-M1L and Kir-M2L are largely -helical.

Structural fluctuations

The root mean square fluctuations (RMSFs) of the C atoms from their time-averaged positions (Fig. 5) are relatively small (ranging from 0.04 nm to 0.16 nm) but show some variation along the length of the helix. For KcsA-M1 and KcsA-M2 the central cores of the helices have C RMSFs of ~0.05 nm, which rise slightly for the terminal residues. For the Shaker helices the pattern is a bit more complex. S5L shows a steady rise in C RMSF from residue 20 onward, suggesting that the C-terminal half of the helix is a little more flexible than the N-terminal half. However, this effect is small. In contrast, S6L shows relatively high (>0.1 nm) fluctuations for residues within the TM helix in addition to those at the termini. In particular, S6L shows a region of dynamic flexibility around the middle of the TM segment. The RMSF profiles for Kir-M1L and Kir-M2L are very similar to those of the corresponding KcsA helices. It is interesting to note that a common pattern seems to emerge if one averages across all six helices: the termini and (to a lesser extent) the center of the helix have a higher RMSF than the two regions flanking the center. We note that Shen et al. (1997) suggested a degree of flexibility in the center of an Ala₃₂ -helix spanning a DMPC bilayer.

View larger version (22K): [in this window] [in a new window]   FIGURE 5   C RMS fluctuations from mean structure versus residue number for (A) KcsA-M1 (bold solid line), Kv-S5L (narrow solid line), and Kir-M1L (dotted line) and (B) KcsA-M2 (bold solid line), Kv-S6L (narrow solid line), and Kir-M2L (dotted line).

Hinge bending

There has been some discussion of the role of hinge-bending motions in the context of possible gating models of channels (Kerr et al., 1996). Furthermore, Perozo et al. (1998) have provided experimental evidence for motions of the pore-lining M2 helices during KcsA gating. It is therefore of interest to look for possible hinge regions in K channel helices. One way of visualizing this, which has been employed in, for example, NMR studies of peptide helix conformations (Dempsey et al., 1991) and in MD studies of the simple channel-forming peptide alamethicin (Tieleman et al., 1999c), is to superimpose snapshots of the peptide structure taken from the MD trajectory on, for example, the N-terminal half of each helix. It should be emphasized that this is simply a way of displaying the consequences of hinge-bending motions and does not, for example, imply greater stability of the N-terminal half of the helix. From examination of the structures thus superimposed (Fig. 6) it would seem that KcsA-M1 and Kv-S5L (i.e., the outer helices) both show less bending motion than KcsA-M2 and Kv-S6L (i.e., the inner helices). For Kir6.2 the difference is less pronounced, with both helices showing some bending motion. Thus for both M2 helices, and for the corresponding S6L helix from Shaker, there is strong evidence for hinge bending.

View larger version (19K): [in this window] [in a new window]   FIGURE 6   C traces of KcsA-M1 (A), Kv-S5L (B), Kir-M1L (C), KcsA-M2 (D), Kv-S6L (E), and Kir-M2L (F). For each simulation 11 C traces, corresponding to structures saved at intervals of 100 ps, are shown. In each case the C atoms were superimposed for the N-terminal half of the helix.

Such hinge-bending motions may be analyzed more quantitatively by combining essential dynamics (Amadei et al., 1993) with the DYNDOM analysis of Hayward and Berendsen (1998) to reveal the major components of the intrahelix motions. The results of this analysis (Table 2) show that for both M2 helices the hinge is close to a glycine residue (G19 in KcsA, G13 in Kir6.2), whereas for S6L it is one turn of the helix before the highly conserved PVP sequence motif. In all three cases the hinge-bending motion results in a significant excess of interdomain over intradomain motion. Not surprisingly, this ratio is greatest for S6L, where the PVP motif is also associated with pronounced loss of helicity in the middle of the TM segment. The time dependence of bending about the central hinge of S6L can be observed by following the end-to-end distance of the helix (Fig. 7). At the outset of the simulation the helix is relatively straight. At 300-350 ps a pronounced kink develops, which is then maintained for the remainder of the simulation. In contrast, the end-to-end distance of S5L does not show any such dramatic changes. Of course, to fully characterize the motion of S6L as "hinge bending," one would like to see several oscillations. This was indeed the case in earlier in vacuo simulations (Kerr et al., 1996). It is perhaps not surprising that the lipid environment has slowed the frequency of such oscillations. Much longer simulations would therefore be required to observe them more fully.

View larger version (25K): [in this window] [in a new window]   FIGURE 7   End-to-end distance of the S5L (dotted line) and S6L (solid line) helices versus time. Snapshots of the S6L C trace are shown for t = 100 and 500 ps (as indicated by the gray arrows).

Interactions of side chains with their environment

In the context of the two-stage folding model for membrane proteins, it is of interest to examine the interactions between the K channel TM helices and their bilayer environment. In Fig. 8 we examine the total numbers of H-bonds from side-chain atoms to polar atoms of the lipids, i.e., to oxygens of the phosphate, glycerol, and acyl moieties. It can be seen that each helix forms between about one and six H-bonds to the lipid headgroups, although this fluctuates dynamically on a time scale of several hundred picoseconds. The overall numbers of such H-bonds are smaller for the Kir6.2 helices, reflecting the smaller numbers of polar side chains at the termini of Kir-M1 and Kir-M2 than for the other helices. The major H-bonding interactions are listed in Table 3. From this it is evident that persistent H-bonds are formed by polar side chains to both lipid headgroups and water molecules in the interfacial region. In particular, arginine side chains in both KcsA-M2 and Kv-S5L form H-bonds to both lipid and water. Amphipathic aromatic side chains (tryptophan and tyrosine) clearly play an important role in hydrogen bonding to the lipid, as do basic side chains. An example of such interactions is shown in Fig. 9, which shows H-bonding interactions of Kv-S6L with two lipid molecules. The lysine (K4) side chain can be seen to be "snorkeling" up to the headgroup region, where it H-bonds to both a phosphate and an ester oxygen. Such a role for lysines at the ends of TM helices has been suggested by a number of authors (Segrest et al., 1990; Mishra et al., 1994; Monne et al., 1998) and has been observed in MD simulations of simple, synthetic TM peptides (Belohorcova et al., 1997), where H-bonds to headgroup-region water molecules were observed. The tryptophan side chain (W2) forms a H-bond to the ester oxygen of an adjacent lipid molecule.

View larger version (58K): [in this window] [in a new window]   FIGURE 8   Protein/lipid H-bonds versus time for KcsA-M1 (A), KcsA-M2 (B), Kv-S5L (C), Kv-S6L (D), Kir-M1L (E), and Kir-M2L (F).

View larger version (39K): [in this window] [in a new window]   FIGURE 9   Snapshot of peptide/lipid H-bonds for Kv-S6L at the end of the simulation. The K4 side chain forms H-bonds to phosphate (p) and ester (e) oxygens of the lipid headgroup, whereas the W2 side chain forms a single H-bond to an ester oxygen.

Such interactions of amphipathic aromatic side chains with the lipid headgroups have been suggested (Schiffer et al., 1992) to "lock" a transmembrane helix or protein into the lipid bilayer. These interactions, especially for tryptophan side chains (and for indole as a simple analogue of these), have been the subject of a number of experimental (Yau et al., 1998) and computational (Woolf and Roux, 1994, 1996; Woolf, 1998; Grossfield and Woolf, 1998; Arkin and Brunger, 1998) investigations. We examined the motions of three aromatic rings (F1, W2, and F29) of S6L. For this we used the method described by Tieleman et al. (1999b), in which the orientation of an aromatic side chain is defined in terms of two order parameters, S_NORMAL and S_LONG (as defined in the caption to Fig. 10). From the results presented in Fig. 10 it can be seen that the aromatic rings at the water/lipid interfaces undergo rapid (on a time scale from 1 to 10 ps) but small fluctuations about orientations that remain fixed for several hundred picoseconds. Major switches in orientation thus occur infrequently but can clearly be seen for F1 and W2. It would seem that to obtain a statistically valid description of such changes in side-chain/bilayer interactions, simulations at least one order of magnitude longer than those presented here will be required. However, inspection of side-chain rotamers at the start versus the end of the simulations reveals that interconversion of side-chain rotamers does occur on a subnanosecond time scale, albeit rarely. Such interconversions are seen, for example, for arginine side chains (in M2 and S5) and for aromatic side chains (in S6).

View larger version (45K): [in this window] [in a new window]   FIGURE 10   Aromatic side-chain orientations for Kv-S6L. The order parameters (solid line, SNORMAL; broken line, SLONG) for the aromatic rings of residues F1 (A), W2 (B), and F1 (C) are shown. SNORMAL is defined in terms of NORMAL, the angle between a normal vector to the plane of the aromatic ring and the perpendicular to the bilayer plane, by SNORMAL = 1/2(3cos2NORMAL  1). SLONG is derived in a similar manner from the angle between long axis of the aromatic side chain (i.e., the vector between C and C) and the bilayer normal. A value of 1 for SNORMAL implies that NORMAL = 0°, i.e., the plane of the aromatic ring is parallel to plane of the membrane, while a value of 0.5 (NORMAL = 90°) indicates the aromatic plane is perpendicular to membrane plane. Similarly, if SLONG is 1, then the side-chain long axis is parallel to the bilayer normal, while a value of 0.5 implies that the long axis of the ring is perpendicular to the bilayer normal.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

What is the continued relevance of MD simulations of isolated TM helices, now that simulations of complex integral membrane proteins (e.g., OmpF) (Tieleman and Berendsen, 1998) and of helix bundles (Tieleman et al., 1999a) are possible? In the context of the two-state folding model for membrane proteins (Popot and Engelman, 1990), one might expect that MD simulations of isolated helices could aid in optimal definition of TM helix lengths. Here we have shown, for pore-lining helices of Kv and Kir channels, that helices longer than the consensus hydrophobic core predictions are needed for stable transbilayer helices to be formed. Indeed, if the shorter helices are embedded in a lipid bilayer, the significant degree of helix/bilayer mismatch (Mouritsen and Bloom, 1984, 1993; Killian, 1998) appears to result in unwinding of the helix termini to allow the peptide to span the bilayer. A similar distortion of "short" TM helices has been observed in MD simulations of simple viral ion channel proteins (Fischer et al., manuscript submitted for publication).

To test whether a non-TM region would or would not form a stable helix when simulated in a transbilayer environment, we took a 21-residue sequence (AVYFAEAGSENSFFKSIPDAF = "non-TM" in Table 1) that corresponds to the end of S5L, the intervening loop, and the beginning of the P-helix from the Shaker Kv sequence. This was built as an -helix, embedded in a POPC bilayer, and simulated for 1 ns in exactly the same manner as for the TM helices discussed in this paper. As can be seen from Fig. 2, C and D, this sequence does not form a stable helix in a TM environment. Instead, this non-TM helix uncoiled fairly early on (100-200 ps) during the simulation and remained thus until the end of the 1-ns period. The loss of helicity (as seen from DSSP analysis; not shown) was greater than that observed for any of the other simulations, as was the final C RMSD (Table 1). It seemed that the loss of helicity was due to the attempt to bury a cluster of hydrophilic side chains in the bilayer environment. This control calculation emphasizes the importance of hydrophobic matching between proteins and lipids for the integrity of membrane protein structure (White, 1994; Killian, 1998).

Toward the end of this study, the x-ray structure of the KcsA channel protein was published (Doyle et al., 1998). This enabled us to compare the M1 and M2 helices at the end of the MD simulations with the corresponding regions in the x-ray structure of the intact protein. For M1 the RMSD was 0.11 nm for C atoms and 0.23 for all atoms. For M2 the corresponding figures are 0.11 and 0.31 nm. The corresponding C RMSD value for comparing the initial (in vacuo generated) helices with the x-ray structure was ~0.15 nm for both helices. This suggests that MD simulation in an explicit lipid/solvent environment may relax an initial -helical model to bring it closer to the crystal structure of the intact protein. However, a more exhaustive set of simulations for all TM helices (now ~100) of known three-dimensional structure would be required to establish this with statistical certainty, which is beyond current computational resources. We also compared the secondary structure and backbone torsion angles of the modeled helices and the crystal structure. These were generally quite similar, with slight deviations from crystal structure observed only at the termini of the M1 helix. For example, the  angle for A27 in M1 is 70° and 49° for the model helix at the start and end of the simulation, whereas for the crystal structure it is 114°. Such deviations at the termini may be a consequence of interactions with the lipid polar headgroups and the water molecules. We also compared helix tilt. The M1 and M2 helices in the crystal structure are considerably tilted with respect to the bilayer normal (by almost 30° in the case of M2). In contrast, the single helices in the present simulation remain more or less parallel to the bilayer normal. This suggests that, as one might expect, helix tilt angles are probably dominated by protein/protein interactions rather than helix/bilayer interactions.

Simulations of isolated helices may also reveal aspects of TM helix dynamics and/or structural stability which, although present in the corresponding intact proteins, are modulated by the remainder of the intact protein. For example, in the current study (and in an earlier in vacuo study; (Kerr et al., 1996) we have shown a pronounced hinge-bending motion of the S6 helix from a Kv channel. It is likely that in an intact membrane protein such conformational transitions may be coupled to changes in the rest of the protein and so might be anticipated to occur on a somewhat longer time scale.

Hinge-bending motions, which seem to be most pronounced in KcsA-M2, Kv-S6, and Kir-M2, may be of importance in the mechanism(s) of channel gating. This is supported by various experimental data. Spin-labeling studies (Perozo et al., 1998, 1999) suggest that activation of KcsA at low pH is linked to movement of opposing M2 helices away from one another. For Kv channels, the results of Armstrong (1971) and Liu et al. (1997), for example, suggest that the S6 helices may move away from one another on channel activation so as to provide relatively unhindered access to the interior of the channel. As discussed above, the PVP motif, which is responsible for the flexibility of Kv-S6, is highly conserved in voltage-gated K⁺ channel sequences, suggesting that it may have some functional significance. Furthermore, simulations of the intact KcsA channel protein (Shrivastava and Sansom, 2000) suggest that small changes in the conformation of M2 may be linked to channel gating. Together, these observations suggest that the dynamic properties of single M2/S6 helices revealed in the current studies may be of some biological relevance. In particular, on the basis of the current studies we predict that the PVP motif in S6 provides a "molecular hinge" that plays a pivotal role in channel gating, analogous to that which has been suggested (Unwin, 1995) for the pore-lining M2 helix of the nicotinic acetylcholine receptor. It would be of interest to see, for example, whether introduction of a similar motif into M2 of KcsA would lead to changes in channel gating that might be possible to characterize via experimental structural studies.

Our studies also provide further simulation data on the nature of lysine/phospholipid and tryptophan/phospholipid interactions. As seen by Belohorcova et al. (1997), for example, lysine residues can "snorkel" (Segrest et al., 1990; Mishra et al., 1994) to the surface to form favorable interactions with phosphates of the lipid headgroups. As now seen in a number of MD simulations (Woolf and Roux, 1994, 1996; Woolf, 1998; Grossfield and Woolf, 1998; Forrest et al., 1999), the NH of tryptophan side chains may form a H-bond to a lipid polar moiety, while the rest of the ring stays in the hydrophobic core. This supports the proposed role of amphipathic aromatic side chains in "locking" integral proteins into a bilayer (Schiffer et al., 1992; Yau et al., 1998).

It is important to consider, albeit briefly, the limitations of our simulations. One limitation is the use of just a single lipid species, POPC. It will be of some interest to see whether our conclusions on the enhanced stability of S6L versus S6S, for example, are robust to changes in lipid headgroup and acyl chain length. Now that MD simulations in a bilayer are not too costly in terms of cpu time, it will be of immense interest to study peptide/bilayer interactions in a range of lipids. As mentioned above, a limitation in terms of exploring, for example, dynamic changes in aromatic ring orientation is the length of the simulations. It should be feasible to extend at least a few simulations by an order of magnitude. Third, long-range electrostatic interactions have been treated approximately (i.e., by use of cutoffs). It will be of interest to see what changes in conclusion, if any, arise from the use of more sophisticated approaches to long-range electrostatics, such as Ewald summation (Tobias et al., 1997; Tieleman et al., 1997).

Finally, our results on isolated pore-lining helices suggest that changes in helix conformation, through hinge-bending motions, may play a role in channel gating. It is now feasible to explore this hypothesis in more detail by simulations, in a lipid bilayer, of the intact KcsA channel protein (Shrivastava and Sansom, 2000) or with homology models of Kir channels (Capener and Sansom, unpublished results).