INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Ion channels are formed in lipid bilayers by integral membrane proteins and enable selected ions to move rapidly (~10⁷ ions s¹ channel¹) down their electrochemical gradients. Ion channels are important in numerous cellular processes, principally electrical signaling (Hille, 1992), and other functions, such as uncoating of viral genomes (Sansom et al., 1998). To understand the physical events underlying the biological properties of channels, one must characterise their structures and their dynamic behavior. However, because ion channels are membrane proteins, we remain ignorant of many of their three-dimensional structures. Indeed, high-resolution crystallographic structures are known for only two ion channels, namely KcsA, a bacterial K⁺ channel (Doyle et al., 1998), and MscL, a mechanosensitive channel (Chang et al., 1998). Furthermore, in both cases it seems that the x-ray structures may correspond to the closed (or largely closed) forms of the channels. This lack of structural data reflects a more general problem for membrane proteins. Although integral membrane proteins are thought to comprise ~20-30% of most genomes (Arkin et al., 1997; Boyd et al., 1998; Wallin and von Heijne, 1998), high-resolution structures have been solved for only a small number of membrane proteins.

Whereas simple systems such as single helices have been simulated in lipid bilayers (Belohorcova et al., 1997; Forrest et al., 1999; Shen et al., 1997; Tieleman et al., 1999b; Woolf, 1997), molecular dynamics (MD) simulations of functional helix bundles in lipid bilayers are still nontrivial (Edholm et al., 1995; Randa et al., 1999; Tieleman et al., 1999a). The M2 protein from influenza A provides a useful and simple system for the development of such techniques. The M2 channel forms proton channels within the viral membrane, which are activated by low pH. It is involved in two stages of viral replication and is essential for virus function. Electrophysiological studies have shown that M2 can be blocked by antiviral drugs such as amantadine (Chizhmakov et al., 1996; Wang et al., 1993). M2 is a small membrane protein, consisting of 97 residues with the TM segment located toward the N-terminus. Spectroscopic (CD and solid-state NMR) studies indicate that the TM segment of M2 is -helical (Duff et al., 1992; Kovacs and Cross, 1997). The functional channel is formed by parallel homotetramers of TM helices (Sakaguchi et al., 1997). A number of studies on synthetic peptides corresponding to the TM helix of M2 suggest that they successfully mimic the intact protein in terms of its 1) -helical conformation (Duff et al., 1992; Kovacs and Cross, 1997; Kukol et al., 1999), 2) ion channel formation (Duff and Ashley, 1992), and 3) tetramerization (Salom et al., 1999). Thus M2 is a relatively simple membrane protein that lends itself to analysis by simulation, complementing a considerable body of experimental data (Sansom, 1998).

Molecular modeling studies of the M2 protein have been carried out by two research groups (Pinto et al., 1997; Sansom et al., 1997). Comparisons of independently derived models suggest convergence to a common structure, which is suitable for investigation in more detail (Forrest et al., 1998). Furthermore, this model structure resembles those derived from experimental data, from solid-state NMR (Kovacs and Cross, 1997), and from infrared spectroscopy (Kukol et al., 1999). The modeling studies suggest that a column of waters can form down the channel, interrupted only by the ring of four His³⁷ residues. Significantly, His³⁷ has been suggested to play a key role in channel activation by low pH (Wang et al., 1995). This evidence suggests a possible mechanism for proton transport via formation of a proton wire (Schmitt and Voth, 1998), with the His³⁷ acting as a selectivity filter (Shuck et al., 1999). However, modeling studies of the channel have so far only modeled the M2 helix bundle in vacuo and included interhelix distance restraints to maintain bundle integrity. Such in vacuo simulations are limited in that they do not include the phospholipid bilayer and water environment experienced by the protein. Simulations by Zhong et al. (Newns et al., 1999; Zhong et al., 1998) have attempted to address the limitations of in vacuo systems, by simulating the M2 helix bundle within a bilayer mimetic (an octane slab and solvated on either side with TIP3P waters). However, the absence of the phospholipid headgroups for such simulations may be problematic. Thus to complement other recent protein/lipid/water simulations (Belohorcova et al., 1997; Randa et al., 1999; Shen et al., 1997; Tieleman and Berendsen, 1998; Tieleman et al., 1999a; Woolf, 1997), made possible by improvements in molecular dynamics algorithms and in computing power, we have carried out simulations of the M2 helix bundle in an explicit phospholipid/water environment.

Recent MD simulations of different lengths of monomeric M2 helices have been carried out in a 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphatidylcholine (POPC) bilayer, to investigate the most probable length of the transmembrane section (Forrest et al., 1999). These results suggest that a 22-residue section of M2 is -helical within the bilayer environment. So tetrameric bundles of 22-residue helices been generated, and simulations have been carried out on these when embedded in an explicit lipid bilayer. The results have been compared with those of an 18-residue four-helix bundle (as investigated in previous in vacuo simulations; Sansom et al., 1997), also within a bilayer environment. We have also investigated a channel model of 18-residue helices generated on the basis of a global searching MD protocol restrained by experimental data from Fourier transform infrared-attenuated total reflection (FTIR-ATR) (Kukol et al., 1999) and a second 22-residue helix bundle model, the conformation of which was biased toward this.

Replacing in vacuo interhelix distance restraints with an explicit lipid membrane environment results in systems that are strongly dependent on the initial model structure. The dynamic behavior of the helix bundles may be interpreted as revealing "breathing" motions of the closed channel. With this in mind, we investigate the stability of the different systems and the structure and dynamics of the pore region of the M2 protein. Furthermore, the behavior of the water within the channels is studied to further understanding of the proton transport mechanism.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Generation of M2 bundle models

Four models of the tetrameric M2 helix bundle have been investigated (see Table 1). Models 18merA and 22merA were generated using restrained in vacuo MD and a simulated annealing (SA-MD) protocol as previously described by Sansom et al. (1997). Briefly, this involves generation of C templates for idealized tetrameric bundles of -helices with the helices orientated such that residues Ser³¹, Gly³⁴, and His³⁷ line the central pore (Bauer et al., 1999; Holsinger et al., 1994; Shuck et al., 1999). These C templates were used in a two-stage SA-MD method, using n  n + 4 distance restraints to maintain -helical backbone conformations and interhelix distance restraints to maintain the integrity of the four-helix bundle. In the case of the 22merA model, noncrystallographic symmetry restraints (Brünger, 1992) were introduced to increase the fourfold rotational symmetry of the bundles. The target distances of the interhelix restraints were slightly increased in 22merA relative to 18merA. Each run of the SA-MD procedure yielded an ensemble of 25 structures, from which the structure with the highest fourfold symmetry was selected as the starting point for extended MD simulations in a bilayer/water environment. The sequences used (Table 1) were taken from the Weybridge strain of influenza A virus.

Model 18merB used the sequence from the Udorn strain of the virus (see Table 1), which differs from the Weybridge strain sequence at two positions (I28V and F38L) (Wang et al., 1993). The model was generated to agree with site-directed infrared dichroism spectroscopy studies (Kukol et al., 1999) on the transmembrane domain of M2 reconstituted in dimyristoylphosphocholine (DMPC) vesicles. These data were utilized in an in vacuo molecular dynamics protocol as described in detail by Kukol et al. (1999). Initial idealized four-helix bundle templates were generated, and a symmetrical MD search was carried out whereby the helices were rotated in a stepwise fashion. This involved a simulated annealing procedure incorporating the spectroscopically determined restraints and additional restraints maintaining the helix geometry and interhelix distances. The resultant ensemble of structures formed clusters at certain helix rotational angles, from which average structures were calculated for each cluster, and a further simulated annealing procedure was then applied. During all MD and energy minimization stages, the helix tilt and rotational pitch angles from the experimental data were incorporated into the overall energy function in the form of a parabolic penalty term.

Finally, model 22merB (Table 1) has the same sequence as 22merA. However, its structure was generated by weakly restraining its C backbone to the experimentally derived structure of 18merB during the final dynamics stage of the SA-MD, resulting in an experimentally restrained 22-residue model. Thus 18merA and 22merA are models based solely on published mutation data, whereas models 18merB and 22merB are based on spectroscopic information.

pK_A calculations

The ionization states of the members of the ring of Asp²⁴ residues at the N-terminal mouth of both 22-residue models were investigated using a protocol for calculating pK_A values of rings of ionizable side chains in ion channel models (Adcock et al., 1998). First, G_BORN, the contribution to pK_A shift due to the protein and bilayer environment, and G_BACK, the contribution due to interaction of the residue with nontitratable charges, were used to calculate pK_A,INTRINSIC:
where pK_{A, MODEL} is the pK_A of an isolated amino acid in free solution.

Second, the pK_{A, INTRINSIC} value was used to calculate the probability of a residue existing in its ionized state, p(x):

where  = RT¹ and x is an N-element state vector, the elements of which are either 0 or 1, depending on whether the residue is un-ionized or ionized, respectively.  = 1 for a basic residue and  = +1 for an acidic residue. G_{i, k} is the screened Coulombic interaction energy between pairs of ionizable residues i and k. The values of p(x) were used to generate titration curves, from which absolute pK_A values were obtained.

Setup of bundle/bilayer/water simulation systems

The set-up of the bundle/bilayer/water system for prolonged MD simulations was essentially as described by Tieleman et al. (1999b). M2 helix models were embedded in a preequilibrated lipid bilayer consisting of 110 molecules of POPC. A hole was generated in the bilayer by a short MD simulation during which a cylindrical radial force was applied to repel lipid molecules. The four-helix bundle was then placed within the hole. Each system was fully solvated with SPC waters and then energy minimized. System details are given in Table 1. In the 22merA and 22merB systems, there was an overall charge of 1 due to the unprotonated Asp²⁴ residue. Thus a Na⁺ counterion was added in each case by replacing a single water molecule at positions corresponding to the lowest Coulombic energy of the ion. Each system was once more energy minimized, followed by an MD equilibration stage of 100 ps, during which the backbone atoms of the protein were restrained to their initial positions. Production runs consisted of a further 2 ns of unrestrained MD. An example snapshot of one system is shown in Fig. 1.

View larger version (117K): [in this window] [in a new window]   FIGURE 1   Snapshot of the protein/bilayer system for 18merA at t = 2 ns. Only the C backbone of the protein is displayed, the carbonyl oxygens of the POPC molecules are shown as space-filling atoms, and water has been omitted for clarity. w represents the water region, p represents the phospolipid headgroup region, and h represents the extent of the hydrophobic region of the system.

MD simulations

MD simulations were carried out as described previously (Forrest et al., 1999; Tieleman et al., 1999a,b), using periodic boundary and constant pressure conditions. A constant pressure of 1 bar was applied independently in all three directions, using a coupling constant of t_P = 1.0 ps (Berendsen et al., 1984), allowing the bilayer/protein area to adjust to an optimum value. Water, lipid, and peptide were coupled separately to a temperature bath (Berendsen et al., 1984) at 300 K, using a coupling constant t_T = 0.1 ps. Long-range interactions were dealt with by using a twin-range cutoff: 1.0 nm for van der Waals interactions and 1.7 nm for electrostatic interactions. The time step was 2 fs, with LINCS (Hess et al., 1997) used to constrain bond lengths, and the force field was based on GROMOS 87 (Hermans et al., 1984).

The SPC water model (Berendsen et al., 1981) was used, as this has been suggested to be preferable to SPC/E when interfaces are simulated as discussed in detail by van Buuren et al. (1993). The lipid parameters are based on those of Berger et al. (1997). The latter authors simulated a bilayer of 64 dipalmitoylphosphatidylcholine molecules and showed good agreement of simulated and experimental lipid densities (and hence area per lipid) and hydrocarbon chain order parameters. The same POPC blayer and parameters have been used in previous MD studies of peptide/bilayer systems (Forrest et al., 1999; Randa et al., 1999; Tieleman et al., 1999a,b).

Computational details

MD simulations were carried out on a 10-processor, 195-MHz R10000 Origin 2000 and on a 72-processor, 195-MHz R10000 Origin 2000 and took ~8 days per processor per 1-ns simulation. Simulation and analysis were carried out using the GROMACS (Berendsen et al., 1995) suite (http://rugmd0.chem.rug.nl/~gmx/gmx.html). Pore radius profiles were calculated using HOLE (Smart et al., 1993). Other analysis used in-house programs. Peptide bond order parameters were calculated for N-H bonds using the following formula:
where  is the angle of the backbone N-H bond to the bilayer normal. Initial models were generated using Xplor (Brünger, 1992) or CNS (Brünger et al., 1999). Structures were examined using Quanta (Biosym/MSI) and Rasmol, and diagrams were drawn using MolScript (Kraulis, 1991).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Helix bundle models

The initial structures of the four protein models are shown in Fig. 2. All of the systems can be seen to contain a ring of Ser³¹ residues and a ring of His³⁷ residues, all of the side chains of which are oriented toward the pore. In the 22mer models, 22merA and 22merB, a ring of Asp²⁴ residues is present at the N-terminus. It is possible to see a greater degree of supercoiling in model 18merB, the structure of which was based on experimental evidence and which was generated using a protocol different from that of the other models.

View larger version (69K): [in this window] [in a new window]   FIGURE 2   Snapshots of the protein models at t = 0 ns, displayed using ribbons. Residues Asp24, Ser31, and His37 have been drawn explicitly. In each case the N-termini (extracellular end) of the helices are at the top of the diagram.

pK_A's of Asp²⁴ residues

The N-terminal mouth of the M2 four-helix bundle contains a ring of Asp²⁴ residues (Fig. 2). Similar rings of acidic side chains are found in the C-terminal mouth of bundles of alamethicin helices (Tieleman et al., 1999a) and either mouth of the pore-lining M2 bundle from the 7 nicotinic receptor (nAChR) channel (Adcock et al., 1998). Studies of the latter channels have indicated that the close proximity of the residues and their position at the C-terminus of the aligned -helix dipole result in an increase in their pK_A values. This would mean that they would not all be fully ionized at pH 7. For influenza A M2, it is more difficult to predict the ionization state(s) of the Asp²⁴ side chains. Their location at the N-terminus of an -helix dipole would promote ionization, which interactions across the Asp²⁴ "ring" are likely to suppress. Thus it is important to estimate ionization states of these residues within their particular protein/bilayer environment.

The titration curves of the four Asp²⁴ residues for an M2 helix bundle are shown in Fig. 3. The pK_A values of each Asp residue in four models (generated as preliminary structures for the 22merA model) are given in Table 2. Note that the corresponding residues of different helices do not have identical pK_A values. This reflects the lack of exact rotational symmetry of the side chains in each bundle. All of the values calculated are much higher than the experimentally determined value of pK_A = 4.4 for an isolated aspartic acid molecule in solution (Stryer, 1988).

View larger version (12K): [in this window] [in a new window]   FIGURE 3   Titration curves for the Asp24 residues of one of the four structures on which pKA calculations were carried out (model 1 in Table 2). The Asp24 residues come from helices H1 (solid line), H2 (dotted line), H3 (dashed line), and H4 (long-dashed line). The four models were taken from the ensemble output from SA-MD in the preliminary stages of generation of the 22merA model. The pKA values are measured by definition at p(x) = 0.5.

Although the accuracy of such pK_A calculations is limited (Warwicker, 1999), a general pattern does emerge. The values in Table 2 suggest that all of the residues will be uncharged at pH 7. However, averaging across the titration curves of all four structures on which the calculation was carried out, the net charge of the Asp²⁴ ring is ~ 1, i.e., only one of the four Asp²⁴ residues will be fully ionized, while the other three will remain protonated. This assignment of ionization states was used in the subsequent bilayer MD simulations. Previous simulations of M2 helix bundles in a bilayer mimetic environment (Newns et al., 1999) did not incorporate the influence of such effects and were carried out with all ionizable residues (with the exception of His³⁷) in their fully charged states. It is likely that the resultant electrostatic force could cause some repulsion of the helices. Such an effect was observed in preliminary simulations of POPC-embedded M2 helix bundles with all four of their His³⁷ residues fully protonated. (Forrest and Sansom, unpublished results).

Bundle stability

The relative stabilities of the various helix bundle models within the bilayer were compared using the root mean square deviations (RMSDs) of their C backbones from the initial structure for each model over time (Fig. 4). All four simulations exhibit a initial increase in RMSD over the first 100 ps or so, after which the RMSDs continue to gradually increase. The final RMSDs are all in the range of 0.2-0.25 nm, values comparable to those found for bilayer simulations carried out on x-ray starting structures of proteins, such as the OmpF porin trimer (Tieleman and Berendsen, 1998), the bacterial K⁺ channel, KcsA (Shrivastava and Sansom, manuscript submitted for publication), and FhuA (Sansom, unpublished results). The one exception is model 22merB, the final RMSD of which is ~0.3 nm. This would appear to be due to a substantial conformation change in the first 50 ps, in which its helices tend toward a dimer-of-dimers structures (see below) and after which the RMSD levels out in a way similar to those of the other three models.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   C RMSD versus time for the four models.

Another indicator of structural stability is the root mean square fluctuation (RMSF) of the C atoms from their average values over the duration of the simulation (Fig. 5). All four models exhibit greater fluctuations at their helix termini, possibly because of their increased number of possible interactions with water and lipid headgroups. Interestingly, these fluctuations seem to be most marked for the 18merB and 22merB simulations, i.e., those based on the experimentally derived models. However, the values of RMSF for the core residues are predominantly less than 0.1 nm, indicating a general stability of all four models.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   C RMSFs for the four models. The extents of the four helices in each bundle are indicated with dotted lines.

Overlaying the C backbones of the structures of 200-ps snapshots through the simulation (Fig. 6) demonstrates the marked fluctuations in structure of helices 2 and 4 from 18merB observed in the RMSF plots (Fig. 5). It is also possible to see the fluctuations in the central residues of helix 3 from 22merB, reflecting a tendency for this helix to bend in toward the central cavity of the helix bundle and confirming the large RMSF values in Fig. 5. This is exemplified by the snapshot at t = 2 ns given in Fig. 7. Such behavior is not observed for the shorter 18merB model. The C traces also demonstrate the lack of a marked left-handed supercoil conformation in 22merA, in contrast to the experimentally derived model, 18merB, which shows a distinctive left-handed supercoil. However, weakly restraining the C backbone of 22merB to that of 18merB during model-building does not appear to have induced a similarly strong effect. The results for both 22mer models suggest that the longer helices do not favor such a tightly coiled conformation, in contrast to the experimental tilt data. This situation is maintained during the bilayer simulation, possibly enhanced by the more complex interactions of the N-terminal protein side chains with their peptide/lipid headgroup/water environment.

View larger version (59K): [in this window] [in a new window]   FIGURE 6   Superimposed C traces of snapshots every 200 ps from the (A) 18merA, (B) 22merA, (C) 18merB, and (D) 22merB simulations. Uppermost are views taken in the x-y plane, with the N-termini at the top. Underneath these are views down z.

View larger version (57K): [in this window] [in a new window]   FIGURE 7   Snapshots of the systems in ribbon format at t = 2 ns. Again the N-termini are shown at the top.

Bundle structure

The observations taken from snapshots of the system can be examined in more detail by measurement of bundle features such as helix tilt, interhelix distances, and helix crossing angles. Experimental studies (solid-state NMR (Kovacs and Cross, 1997) and FTIR-ATR (Kukol et al., 1999)) have also yielded data on the tilt of the M2 helices relative to the bilayer normal. These can be compared to the tilt angles of each helix during the simulation (Fig. 8). Note that the initial tilts vary from ~10° to ~30° for the different systems, i.e., they vary between different starting models. For 18merA and 22merB, the tilt angles average 20° through the simulation and exhibit very little fluctuation around that value. However, by the end of the simulation, 22merB is more tilted than 18merA. 22merA is notable because the values are spread between ~7° and ~30°, indicating a rigid body tilting motion of the helix bundle as a whole, as seen by visual examination of the structures over time. Thus some of the helices would become less tilted than others with respect to the bilayer normal. The helices in 18merB are tilted to a much greater extent, fluctuating around 30°. In each case, the tilts appear to reflect the degree of supercoiling observed in the snapshots and C traces (Figs. 2, 5, and 6).

View larger version (18K): [in this window] [in a new window]   FIGURE 8   Helix tilt values with respect to the bilayer normal versus time for each helix in the four models. Solid line, helix 1; dotted line, helix 2; dashed line, helix 3; long-dashed line, helix 4.

It is interesting to note that the experimental values for tilt angles are 33 ± 3° for solid-state NMR and 32 ± 6° for FTIR, which are somewhat higher than the values from the simulations. However, helix tilt angles are difficult to calculate, being sensitive to the method used, if nonideal helices are involved. This may be explored further by plotting the order parameters of the N-H bonds from the backbone of each residue (S_NH) averaged over the duration of the simulation (Fig. 9). A value of S_NH = 1 indicates that the NH bond is parallel to the bilayer normal, and therefore S_NH decreases as the tilt increases. The S_NH values for the current simulations average around 0.7 but show large fluctuations around that value, following the periodicity of the helix. A similar pattern is seen for a plot of order parameters for the model at the beginning of simulation 18merB, taken from Kukol et al. (1999) (data not shown). Furthermore, there appears to be a gradual increase in S_NH with residue number, indicating that the helices become less tilted toward the C-terminus and that they contain a degree of flexibility.

View larger version (21K): [in this window] [in a new window]   FIGURE 9   Order parameters (SNH) for the backbone N-H bonds versus residue number for the four helix bundles averaged over the length of the simulation and over all four helices. Solid line, 18merA; dotted line, 18merB; dashed line, 22merA; long-dashed line, 22merB.

Stable packing of helices is often mediated by "ridges-in-grooves" interdigitation of side chains. As discussed by Chothia (Chothia et al., 1981), such packing may be identified via analysis of helix crossing angles, i.e., the angle between the long axes of adjacent two helices. Crossing angles, , were calculated for each pair of adjacent helices as a function of time (Fig. 10). In general these plots mirror those for the helix tilts. The 18merB model appears to have the largest helix crossing angles, with an average over all helix pairs of 40°. However, the four lines diverge considerably during the last nanosecond, suggesting a rearrangement of helix packing within this structure. In comparison, 18merA displays smaller crossing angles (average of 28°) but would appear to be more stable. Each of the 22mer models exhibits behavior similar to that of the other. The average over all pairs of helices in 22merA is 25° and for 22merB it is 30°. In both cases, the lines bifurcate. In the case of 22merA, helix 3 moves away from the rest of the bundle (Fig. 7), which affects the crossing angles. However, the bifurcation is most marked in the 22merB simulation, which is likely to be due to the bending motion of helix 3 in this system (Figs. 6 D and 7).

View larger version (17K): [in this window] [in a new window]   FIGURE 10   Helix crossing angles () against time for each pair of adjacent helices. Solid line, helix 1 to helix 2; dotted line, helix 2 to helix 3; dashed line, helix 3 to helix 4; long-dashed line, helix 1 to helix 4.

It is also useful to look at changes in the cross-sectional shape of the helix bundles. Fig. 11 shows distances between helices lying opposite one another across the bundles as functions of time. If the two distances (i.e., between helices 1 and 3 and between 2 and 4) are similar, the helices form the vertices of a square, whereas if the values differ greatly, then the bundle has a trapezoidal cross section, which might be interpreted as the formation of a "dimer of dimers" (Zhong et al., 1998). 22merB shows the most dramatic example of the latter behavior, even though it is initially square in cross section. The formation of this dimer of dimers would appear to explain the sharp increase in RMSD observed in the first 200 ps of the simulation (Fig. 4). Then 18merA forms a dimer of dimers after ~1 ns. The experimentally derived model, 18merB shows some evidence of trapezoidal cross section between ~200 and ~1500 ps, but then reverts to a square cross section during the latter part of the simulation. The 22merA model has much greater interhelix distances than all other models, fluctuating around 1.7 nm for the final nanosecond. The two distances stay relatively similar (i.e., the bundle has a square cross section) during 2 ns of simulation. However, the size of the interhelix distances suggests that the helices are not within van der Waals contact along the lengths of their helices. Thus formation of a tetramer or of a dimer of dimers does not appear to be correlated with the length of the helices. As seen in the in vacuo simulations of Sansom et al. (1997), M2 tends to form a pyramid shape, where its N-termini are closer together than its C-termini (Fig. 6, A, B, D). Because current calculations use the average C positions along the entire length of each helix, this asymmetry explains the differences in values.

View larger version (16K): [in this window] [in a new window]   FIGURE 11   Interhelix distance values versus time for the opposing pairs of helices. Solid line, helix 1 to helix 3; dotted line, helix 2 to helix 4. In each case, the helices were defined using the averages of the C positions over the length of the helix.

The bundles are not affected by independent motions of the helices in the z dimension (in the bilayer normal). Calculation of the average z coordinates of the four His³⁷ residues over time (data not shown) indicates that helix length has no systematic effect on the relative positions of the helices in the bilayer.

Pore and pore water behavior

It has been suggested that the M2 protein may act as a proton channel by providing a column of waters through which H₃O⁺ transfer may occur, via the proton wire (or Grotthüs) mechanism (Gutman and Nachliel, 1997). Thus it is of interest to observe the behavior of the pore and the pore waters in these simulations. Pore radius profiles (Fig. 12) have been calculated using HOLE (Smart et al., 1997) at the beginning and the end of the simulations. In general, the M2 bundles have rather narrow pores. Often the radius of the pore drops below that of a water molecule (~0.16 nm) and the channel is thus effectively occluded by the side chains at these positions. In the case of 18merA there are two large occlusions at both entrances to the channel at the start of the simulation. The occlusions flank a central cavity of radius ~0.3 nm and length ~0.5-1.0 nm. This pocket is seen to house four water molecules throughout the simulation. The movements of these waters appear to be confined to within the pocket, and only in the last 200 ps of the simulation do they exchange with bulk waters, forming a single continuous column of water. Given the experimental evidence on the role of the His³⁷ in activation (Wang et al., 1995), it might be expected that the C-terminus would provide a barrier to proton transport at pH 7, when the channel is likely to be closed. However, 18merA opens up at the C-terminus and remains closed at the N-terminus, as can be seen from the plot for t = 2 ns. It is likely that this situation is transient, as only a single column of waters is formed (and thus would be easily broken) and as the other helix bundle models show predominantly the opposite behavior. Lengthier simulations are necessary to clarify this point.

View larger version (23K): [in this window] [in a new window]   FIGURE 12   Pore radius profiles, calculated at the start, t = 0 ns (solid line) and the end, t = 2 ns (broken line) of the simulation. Gray bars indicate the approximate extents of the bundles, running N-terminus (left) to C-terminus (right). Solid black lines superimposed on the gray bars indicate the approximate His37 positions. The radius of water (~0.16 nm) is indicated by a dotted line.

The pore of model 22merB also appears to be occluded during most of the simulation. It has a pocket of five waters until ~1400 ps, after which a similar single-water column forms, which connects the pocket waters with bulk, this time at the N-terminal mouth. Fig. 13 shows the pore radius profile of 22merB at t = 1400 ps (dotted line) and compares it to those obtained for the start and finish of the simulation. The open N-terminus is maintained for ~500 ps, before it closes up once more, resulting in a pocket containing 10 water molecules. Thus the 18merA and 22merB models appear to exhibit similar behavior.

View larger version (29K): [in this window] [in a new window]   FIGURE 13   Pore radius profiles, showing important time points in the simulations. (A) 22merA with its N-terminus open at t = 0 ps (solid line), at t = 400 ps (dotted line) where its C-terminus is opening, and at t = 2 ns (thicker line). (B) 18merB at t = 0 ps (solid line), t = 50 ps (dotted line), and t = 1 ns (thicker line). (C ) 22merB at t = 0 ps (solid line), t = 1400 ps (dotted line), and t = 2 ns (thicker line).

The 18merB and 22merA models also show some similarities in their pore radius profiles. Both models have large pore radii, particularly at their N-termini. In the case of 22merA, this end is open for the entirety of the simulation, and its C-terminus opens up after ~300 ps (Fig. 13). This reflects its large interhelix distances (Fig. 11) and the motion of one of the helices away from the rest of the bundle (Fig. 7). The N-terminus of 18merB opens up after ~50 ps, as shown by the dotted line in Fig. 13. Its C-terminus becomes considerably wider at ~1 ns, as shown by the thick line in Fig. 13. The columns of water observed contain two to three parallel chains of "water wires" at any one time. However, occlusions remain around the position of the ring of His³⁷ residues, preventing formation of a continuous column of water along the length of the channel, suggesting that the models are of the "closed" forms of the channel, as proton transport could not occur.

The models whose pores are open during the simulations, i.e., 22merA and 18merB, are also the systems whose interhelix distances are similar for opposing pairs of helices (see above). Furthermore, in the models whose pores are occluded (18merA and 22merB), the bundles have a trapezoidal cross section, i.e., two opposing helices move toward each other, into the pore, and the other two helices move apart, elongating the bundle in the other direction. This behavior is illustrated for 22merA in Fig. 7. Thus the "breathing" motions of the M2 protein as it fluctuates between tetrameric and dimer-of-dimers conformations appear to be coupled to formation/breaking of a continuous water column. The helix bundles demonstrate considerable motions, even in the absence of a low pH and consequently in what one might expect to be closed forms of the channel.

As mentioned above, the proposed mechanism for proton transport is via a water wire Grotthüs-type mechanism (Pomes and Roux, 1998). For such a mechanism to occur efficiently, the waters within the pore would need to be effectively "frozen," i.e., showing much reduced diffusion compared to the bulk water. Indeed, within the channel, water diffuses sufficiently slowly that its diffusion coefficient cannot be estimated on a ~10-ps time scale (data not shown). A similar effect is seen in other proton channels, including gramicidin (Pomes and Roux, 1998) and a channel formed by four helices of the synthetic leucine-serine peptide known as LS2 (Randa et al., 1999). The behavior of individual pore waters can also be investigated by following their z coordinates over time (Fig. 14). It can be seen that waters found in the pore undergo very small stepping motions along z and fluctuate around particular positions for time scales of approximately hundreds of picoseconds. Similar behavior is seen for gramicidin, LS2, and synthetic nanotubules (Engels et al., 1995). In comparison, a water molecule that is released into bulk solution (Water 5140) can undergo much greater fluctuations in z and more frequent stepping motions. Again, a very similar behavior was observed for the waters in the LS2 proton channel (Randa et al., 1999).

View larger version (21K): [in this window] [in a new window]   FIGURE 14   Trajectories of the z coordinates of water molecules from 22merA. Waters 5406, 5141, and 5524 undergo small fluctuations around specific z coordinates with small jumps between those positions. Water 5140 is seen to leave the pore at ~750 ps, so that its motions begin to undergo much greater fluctuations with more frequent transitions. The extent of the bilayer is marked by dotted lines.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Comparison with experimental results

It is instructive to compare the simulation results with those of experiments. Solid-state NMR (Kovacs and Cross, 1997) and isotopic-labeling IR spectroscopy data (Kukol et al., 1999) have provided independent evidence of the tilts of the helices of synthetic peptides of M2 in bilayer environments. Solid-state NMR gives values of 33 ± 3°, whereas IR gives 31.6 ± 6.2° with respect to the bilayer normal. In the current work, models 18merA, 22merA, and 22merB all contain slightly less tilted helices (~20° with fluctuations of ~5°) than the experimental data. These values appear to be dependent on the starting conformation; the helices of model 18merB, the initial tilt angles of which were set at 30°, remain at about this value and are in better agreement with the experimental evidence, as might be expected, because its structure was derived from this evidence. Weakly restraining the 22merB to this conformation does not result in such a strongly tilted (or supercoiled) structure, suggesting that the longer helices do not favor such a conformation. There are several important factors, which might cause the experimental values to be somewhat higher than those in the simulation studies. First, tilt values calculated from experimental data were based on an assumption that the helix is rigid. However, calculation of the backbone N-H bond order per residue suggests that the helices are fairly flexible, especially toward the termini of the helices. Thus in the 25-residue-long synthetic peptides (Ser²²-Leu⁴⁶) used for the experiments, this disorder is likely to be greater still. It is not clear how the inclusion of such intrahelix flexibility would influence the interpretation of the spectroscopic data in terms of an overall helix tilt angle. Second, in both experiments, the lipid bilayers were composed of DMPC lipids, which are shorter by ~0.3 nm than the POPC lipids used here. The effect of this may be to cause the tilt (or supercoil) to be increased to compensate for hydrophobic mismatch (de Planque et al., 1998). Further NMR studies on synthetic M2 peptides in DOPC bilayers may help provide an answer (Kovacs et al., 1999). Third, the protein environments during the experiments are at a relatively low level of hydration, which might encourage the helices to tilt more than would otherwise be the case, to optimize their H-bonding interactions. Finally, the time scales of such experiments are around a microsecond, suggesting that some averaging may occur compared to the nanosecond time scale used here. Given these differences, it is perhaps more significant that both studies and experiment agree that the helices are tilted, rather than that they differ over the exact size of the tilt angle, particularly as during the simulations the tilt angles cover a large range of values from ~10° to ~40°.

Different models

MD simulations have shown no obvious effect caused by changing the length of the helices from 18 to 22 residues. It might be expected that the introduction of three polar groups (Ser²², Ser²³, and Asp²⁴) at the N-terminus would strongly influence the interactions of the protein with the lipids and possibly act to stabilize the helix bundle. However, there is no evidence for this in the current work. On the contrary, large RMSDs and considerable distortion of the symmetrical coiled-coil conformation are observed for the 22merB model. Instead, the evidence suggests that 22merA and 18merB are paired in their behaviors, and 18merA and 22merB are also similar. Thus it is difficult to conclude that extending the helices has a consistent effect on the outcome of the simulations.

A further distinction between the models is that the sequence of model 18merB is taken from a different strain of influenza A virus, effectively introducing mutations at positions 28 and 38 of Ile to Val and Phe to Leu, respectively. The latter mutation might well be expected to cause an increase in the pore radius at that position (C-terminal to the His³⁷ residue). Indeed, the pore radius profiles generally show two minima at around this position (z = 0 to +1 nm), and the more C-terminal minimum is larger in the profile of 18merB. However, this does not cause any unexpected behavior in terms of the opening frequency of the C-terminus. The Ile-to-Val mutation might be expected to have less of an effect, because of the similarity of their side chains. Furthermore, these residues are observed to point toward the lipid rather than into the pore. As expected, no significant effect is observed.

The use of different models in multiple short simulations has given trajectories that vary in their stabilities, conformations, and dynamics. One might suggest that if the simulations were extended, by perhaps an order of magnitude, the trajectories might begin to demonstrate overlapping behaviors. More realistically, the different conformations may be allowing better sampling of the range of motions experienced by the system, which would require very long time scales (perhaps on the order of hundreds of nanoseconds) to be observed for a single start structure.

Implications for M2 function

These simulations appear to reveal dynamic fluctuations of the tetrameric bundle of M2 helices, even when all four His³⁷ residues are deprotonated (i.e., the presumed closed state). The ability of waters to form a continuous pore appears to be dependent on the breathing motions of the helix bundle. Under these conditions we see two types of behavior: 1) the channel is occluded at both termini with a pocket of 3-10 trapped waters or 2) the channel contains a column of waters that is sometimes broken at one or the other mouth. Experimental evidence suggests the protonation of one or more His³⁷ side chains is responsible for opening of the channel. Indeed, such channel gating is observed after sequential protonation of two His³⁷ residues for simulations of M2 helix bundles in octane (Newns et al., 1999). Such protonation may be expected to favor the second conformation of the channel over the first.

Proton translocation mechanism

The behavior of the pore waters in the M2 helix bundle models is seen to be similar to that observed in previous simulations of the synthetic leucine-serine peptide, LS2 (Randa et al., 1999). This peptide was also simulated in an explicit POPC bilayer as a four-helix bundle, in which form it is believed to form a proton pore. In both cases, the waters are effectively "frozen," as confirmed by the trajectories of the z coordinates of pore waters over time and by the extremely low diffusion coefficients for water molecules within the pore region. This suggests that both systems may utilize a Grotthüs-type mechanism for proton conduction.

Simulation methodology

It is useful to consider the limitations of the simulation methods used in this paper. A major concern is the length of the simulations. Although 2 ns is a relatively long simulation time by current standards, it is short when compared to biological time scales, e.g., of proton conduction (~100 ns). However, as discussed above, multiple short runs with starting structures may help to overcome this problem. Another concern is that the treatment of long-range electrostatic interactions with a cutoff may be too approximate. Pure bilayer simulations (Tieleman and Berendsen, 1996) using this method have produced reasonable agreement with experiment, but improved methods (e.g., Ewald summation) have been discussed (Tieleman et al., 1997; Tobias et al., 1997) and are becoming more feasible for large and complex systems. A final concern is the treatment of the ionizable residues. Calculations of pK_A values have provided some insight into the probable states of the aspartate side chains but do not take into account possible interaction with charged lipid headgroups. Furthermore, these simulations do not yet incorporate dynamic changes in protonation state during the course of a simulation, although this has been attempted for gramicidin (Pomes and Roux, 1996; Sagnella et al., 1996).

In conclusion, molecular dynamics simulations of protein embedded in a fully atomistic lipid bilayer have been carried out on tetrameric helix bundles corresponding to different lengths of the transmembrane domain of the M2 protein from influenza A virus. We observe motions on a subnanosecond time scale that may represent breathing motions of the helix bundle. These include conformations where the cross section of the helix bundle forms a trapezium and correspond to a "dimer of dimers." Such motions are well sampled by the use of multiple simulations with different starting conformations, as the simulations are extremely sensitive to the initial structures. The occlusions observed in the pores of the helix bundles suggest that despite these motions, the systems represent closed channels, as would be expected at neutral pH. Analysis of the pore waters suggests that their motion is highly restricted and supports the proposal that proton transport occurs via a water-wire mechanism. Future simulations will focus on the effect of low-pH protonation of the His³⁷ residues on the helix bundle structure and dynamics.