INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Many classes of ion channels are formed by bundles of parallel -helices surrounding a central water-filled pore. Examples of such channels include the mechano-sensitive channel MscL (Chang et al., 1998), the nicotinic acetylcholine receptor (Unwin, 1995), and channels formed by peptides such as LS2/LS3 (Lear et al., 1988) and alamethicin. Although several high-resolution structures of ion channels have now been determined (Chang et al., 1998; Doyle et al., 1998; Dutzler et al., 2002), peptide channels, because of their simplicity, remain an attractive model system to study ion channel properties. They also provide an opportunity to study self-assembly of transmembrane helix bundles, a process believed to be important in the folding of membrane proteins (White and Wimley, 1999; Popot and Engelman, 2000).

Alamethicin (Alm) is a 20-residue Aib-rich channel-forming peptide, a member of the family of peptaibols. It has been intensively studied by experimental (e.g., Woolley and Wallace, 1992; Sansom, 1993a; Cafiso, 1994; Duclohier, 2001) and computational methods (Kessel et al., 2000; Tieleman et al., 2001a; Tieleman and Sansom, 2001). It occurs in two native forms, the R_f30 form used in this study (Ac-Aib-Pro-Aib-Ala-Aib-Ala-Gln-Aib-Val-Aib-Gly-Leu-Aib-Pro-Val-Aib-Aib-Glu-Gln-Phol) and the R_f50 form, in which the Glu at position 18 is replaced by a Gln, making the peptide electrically neutral. In addition to these native forms, there are several natural as well as designed mutants, including covalently linked dimers of alamethicin (Woolley et al., 1997) and peptides in which the Aib residues have been replaced by Leu (Sansom, 1993a). Alm forms channels with well-defined conductance levels (Fig. 1) that are generally thought to correspond to channels formed by different numbers of peptides.

View larger version (34K): [in this window] [in a new window]   FIGURE 1   Single channel current trace of an alamethicin channel in a diphytanoyl phosphatidylcholine bilayer with 0.5 M KCl and a 125-mV transmembrane potential. Arrows on the right give the presumed number of helices per bundle for each level of conductance. Data from Mellor and Sansom (Sansom, 1991). Also shown are the starting models for the N4-N8 simulations.

Breed et al. (1997) have constructed models of channels formed by four to eight alamethicin helices. Although there are considerable experimental data that in general support the helix bundle model for Alm channels (for review, see Woolley and Wallace, 1992; Sansom, 1993a; Cafiso, 1994; Bechinger, 1997; Duclohier, 2001), it remains unclear whether the molecular models that have been developed are sufficient to explain the conductance properties of the channels. In a previous study we conducted short (2 ns) simulations of a model of the hexameric bundle in a POPC lipid bilayers (Tieleman et al., 1999a). In this paper we use MD simulations in an attempt to examine the conformational stability of N = 4 to 8 helix bundle models and compare some properties of the simulated models with experimental data. The simulations are of duration comparable with the mean passage time of a single ion through such a channel (~5 ns is equivalent to an ionic conductance of 250 pS at 125 mV). We then extend the simulations for the model (N = 6) that best satisfies our computational and experimental validation, to explore the conformational dynamics over a timescale equivalent to the passage of four or five ions. We also explore the sensitivity of the behavior of this model to changes in simulation protocol to help define optimal simulation conditions for ion channels in general. In this fashion we deliver a validated "best guess" model for at least one conductance level of the Alm channel, which we hope will prove suitable as the basis for future more in depth calculations of channel electrostatics and permeation models (Kuyucak et al., 2001; Tieleman et al., 2001c).

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Starting structures

The creation of the alamethicin starting structures, from the models of Breed et al. (1997) and a pre-equilibrated POPC bilayer, has been described in Tieleman et al. (1998). In brief, channel models were assembled from the alamethicin monomer x-ray structure using a distance-restrained in vacuo MD simulation approach. The principal assumptions in constructing these models, derived from experimental data, were that: 1) the helices within a bundle are all parallel to one another and 2) the Gln-7 sidechain (and thus the polar face (see Breed et al., 1997)) of each amphipathic Alm helix was directed toward the lumen of the channel. As a representative example, a snapshot of the Alm N8 bundle is given in Fig. 2.

View larger version (88K): [in this window] [in a new window]   FIGURE 2   Snapshot of the N8 system starting structure after a 25-ps run with position restraints on the peptide. Gln-7 and Phl-20 are drawn as ball and stick, the rest of alamethicin as ribbons, and the phosphatidylcholine headgroup is highlighted in the lipids.

pKa calculations

In a number of channels, both peptide (e.g., Influenza A M2, Forrest et al., 2000) and protein (e.g., the nicotinic receptor, Adcock et al., 1998) it has been suggested that ionizable residues directed toward the pore lumen, and thus toward one another in a radially symmetric assembly, may have local pKa values that differ from those of the isolated amino acids in bulk solution. Furthermore, in a previous study of the N6 channel, the ionization state of the Glu-18 residues had a marked influence on the conformational stability of the bundle, even in relatively short simulations (Tieleman et al., 1999a). We have therefore attempted to estimate the most likely ionization state of the Glu-18 sidechains in the Alm helix bundles. The pKa calculations were based on numerical solution of the linearized Poisson-Boltzmann equation to estimate the relative free energies of the ionized and protonated states of the Glu-18 sidechains. The details of these calculations have been described in Tieleman et al. (1998). They suggest that the average total charge for the Glu-18 residues in the N5-N8 systems is 1 as opposed to 5 to 8. For N4 the most likely charge state was 0. Therefore, we adjusted the protonation state of the Glu-18 residues accordingly.

Simulation details

Table 1 lists all simulations. Simulations were carried out with Gromacs (Berendsen et al., 1995). The lipid parameters are based on Berger's (Berger et al., 1997), the protein parameters are a modified version of GROMOS87 with all-atom aromatic residues and modified carbon-water oxygen repulsion as described before (Tieleman et al., 1999a), combined with the SPC water model. Simulations N4-N8 and N6-20ns used a twin-range cutoff for Coulomb of 1.0/1.8 nm and a single cutoff for Lennard-Jones interactions of 1.0 nm. For N6-PME and N6-0.5M we used a cutoff of 0.9 nm for Coulomb and Lennard-Jones interactions and Particle Mesh Ewald (Essmann et al., 1995) to calculate the remaining electrostatics contributions on a grid with 0.12-nm spacing. For N6-RF, we used in addition to the 1.0/1.8-nm twin-range cutoff an analytical reaction field correction that assumes that beyond 1.8 nm each atom sees a continuum dielectric environment with a relative dielectric constant of 80. Test simulations with position restraints on the protein were carried out as described in Table 1. Temperature and pressure were controlled using the weak coupling algorithm (Berendsen et al., 1984) with _T = 0.1 ps, T = 300 K, _p = 1.0 ps with a pressure of 1 bar independently in three dimensions. The time step used was 2 fs with a neighbor list update every 10 steps. Coordinates were saved every picosecond. Bond lengths were constrained with the LINCS algorithm (Hess et al., 1997).

Essential dynamics

The essential dynamics analyses use a method recently developed by Hess to characterize how much of the space sampled by a simulation overlaps with the space sampled by a second simulation (or the second half of the first simulation) (Hess, 2001). This measure is given by the overlap between the covariance matrices of two simulations, A and B:
If the covariance matrices A and B describe a completely different sampled space, they are orthogonal and the overlap is 0. If they describe the same space, the overlap is 1. To estimate how much of the major essential modes from the covariance analysis describes functional motion and how much of it resembles random diffusion, its cosine content can be calculated. Hess has shown that for random diffusion, without potential, the first few principal components are cosines with the number of periods equal to half the principal component index (Hess, 2000). The cosine content c_i of principal component i is given by:
with T the length of the simulation, and p_i(t) the amplitude of the motion along eigenvector i at time t. The cosine content varies between 0 (no similarity to a cosine) to 1 (a perfect cosine).

General analysis

Most analyses were done using Gromacs programs; HOLE (Smart et al., 1996) was used to calculate pore radius profiles and to estimate channel conductances; molecular graphics were made with Molscript (Kraulis, 1991).

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Structure and dynamics of the helix bundles

To gain an impression of the conformational stability (on a multinanosecond timescale) of the different models in this study the structural drift, measured as root mean square deviation (RMSD) from the starting model was measured (Fig. 3). Let us first consider the tetrameric bundle (N4). Within the first 500 ps the RMSD for N4 jumps to a value of 0.3 nm, whereas those for the other models rise more slowly to ~0.2 nm over this period. Visual examination of the model shows that the N4 bundle model becomes quite irregular with one peptide losing most of its helicity. This is in marked contrast to the behavior of, e.g., the hexameric model N6 (Fig. 3 A) Furthermore, N4 shows substantially greater conformational drift over such a time period than we have seen for other tetrameric channel models (such as the LS2 proton channel (0.17 nm after 4 ns) (Randa et al., 1999) and the Influenza A M2 proton channel (Forrest et al., 2000)) that were simulated under very similar conditions. This leads us to suspect that a tetrameric Alm helix bundle channel does not form a stable channel or at least that our starting model does not represent even a metastable packing arrangement of 4 Alm helices. Of course, it remains conceivable that tetrameric bundles of Alm, although not forming channels, act as a seed or nucleation for the formation of large (channel) assemblies. Indeed, preliminary results from direct simulations of self-assembly of monomeric Alm helices have revealed some such nucleation events (Tieleman, unpublished data). In contrast to N4, the RMSD curves for the other simulations (N5 to N8) gradually increase to plateaus after 4 ns and appear to have leveled off. For N5, N7, and N8, the plateau RMSD values are between 0.25 and 0.3 nm. This value is comparable with that seen in simulations of membrane proteins based on medium resolution (~0.3 nm) x-ray structures (Tieleman and Berendsen, 1998; Shrivastava and Sansom, 2000) or plausible homology models (Capener et al., 2000). Remarkably, the simulation of N6 shows a lower RMSD than the other systems with a value of less than 0.2 nm after 4 ns. This suggests a degree of conformational stability comparable with that seen, e.g., in simulations of OmpA, a high-resolution -barrel membrane protein (Bond et al., 2002), lending an additional degree of plausibility to this model of the hexameric channel.

View larger version (28K): [in this window] [in a new window]   FIGURE 3   (A) RMSD from the starting model for N4 (solid) and N6 (dashed) with snapshots at 0.5, 2, and 3.5 ns. (B) RMSD from the starting model for N5, N7, and N8. Data are running averages over 100 ps and are calculated from the backbone atom positions after fitting on all C- carbons.

Secondary structure analysis

A more local measure of conformational stability is provided by examination of the time evolution of the secondary structures of the peptides in each bundle (Fig. 4). In general, the C-terminal segments of the peptides are somewhat less helical and distortions of the initial helices tend to develop in the center of the helices. This is as anticipated on the basis of previous simulation studies (Gibbs et al., 1997; Tieleman et al., 2001b) and correlates with some experimental data (Gibbs et al., 1997). In terms of secondary structure N4 is the least stable bundle, with substantial unfolding in the second helix. In both N4 and N5 one of the helices shows substantial 3₁₀ helix content. In N5 significant structural fluctuations in the middle of the helix (close to the Gly-x-x-Pro motif (Sansom and Weinstein, 2000)) are observed in three of five helices, again in agreement with previous simulation and experimental studies (Gibbs et al., 1997; Tieleman et al., 2001b). The same feature is present in helices in N6, N7, and N8. Once again, N6 is rather more stable compared with the other bundles. Thus, both at a local and a more global level the N4 helix bundle seems to be the least stable and the N6 bundle the most stable.

View larger version (21K): [in this window] [in a new window]   FIGURE 4   Secondary structure, calculated by DSSP (Kabsch and Sander, 1983), as a function of time for N4, N5, N6, N7, and N8. Blue corresponds to -helix. Only -helix, 310-helix, bend, and turn are observed. White horizontal bands separate the helices.

Pore properties

Having examined the conformational stability of the various bundles, let us now explore the consequences for stable pore formation. Snapshots of the various bundle structures at the end of the 4-ns simulations (Fig. 5) highlight the backbone structures and the water columns (for N5 to N8) inside the channels.

View larger version (52K): [in this window] [in a new window]   FIGURE 5   Snapshots after 4 ns of simulation of the water columns in the N4, N5, N6, N7, and N8 channel.

Let us first consider N4. It is evident that the structure of N4 at 4 ns is rather irregular. Furthermore, N4 does not contain a continuous water column. This would suggest that N4 does not form a stable channel, and that the lowest conductance level shown in Fig. 1 may therefore correspond to N5. This supports the earlier suggestion based on model building and in vacuo simulations (Breed et al., 1997).

For the other bundles (N5 to N8) well-defined columns of water are present. These columns contain ~57 (N5), 63 (N6), 84 (N7), and 95 (N8) water molecules. Thus, each of the bundles N5 to N8 can be expected to provide adequate solvation to remove the Born energy barrier associated with movement across a lipid bilayer (Parsegian, 1969).

Several sets of measurements of alamethicin conductance levels are available (Hanke and Boheim, 1980; Sansom and Mellor in Sansom, 1991; Vodyanoy et al., 1993; You et al., 1996). One can attempt to relate these measured conductances to the pore models. There has been considerable effort in recent years to predict single channel conductances on the basis of channel structures (for review, see Tieleman et al., 2001c). A relatively simple approach is based on integration along the pore axis of the calculated electrical resistances of a series of electrolyte-filled cylinders that fit with the pore (Smart et al., 1997). This approach has worked well with a number of peptide and other channel models (Smart et al., 1998; Law et al., 2002). To this end, the pore radius profiles were calculated using the program HOLE (Smart et al., 1993, 1996), which uses a Monte Carlo algorithm to find the maximal radius spherical probe that will fit in the channel at a given depth. It also calculates the solvent accessible area using a Connolly surface and calculates the radius of a circle with the same area. Profiles for both methods are shown in Fig. 6 for the N5 to N8 bundles. (The N4 bundle has a "pore" with a radius of ~0.01 nm, i.e., much too small to accommodate an ion and so was not considered further).

View larger version (23K): [in this window] [in a new window]   FIGURE 6   Pore radius profiles, averaged over snapshots at 10-ps intervals for N4, N5, N6, N7, and N8. Dashed lines show the radius been determined using a spherical probe, solid lines show the radius using the "Connolly" option of Hole. In this case, the radius given is an equivalent radius calculated from the area determined using a probe of radius 1.4 A.

For the N5 simulation, the narrowest region of the pore is found around the ring of Glu-18 residues. We note in passing the suggestion that a ring of five glutamate sidechains is present in the pore of the cation selective nicotinic acetylcholine receptor (Adcock et al., 1998). For the other models (i.e., N6 to N8) the narrowest part of the pore corresponds to the ring of Gln-7 sidechains. In the original models, minimal radii are N5 0.08 nm, N6 0.17 nm, N7 0.37 nm, and N8 0.55 nm, compared with, e.g., a value of 0.13 for the Pauling radius of a K⁺ ion. The minimal pore radius of N5, initially smaller than that of K⁺, increases significantly and remains near 0.2 nm over the last 2 ns of the simulation. Interestingly, this is close to the minimal radius of N6 for most of the simulation of the latter. Using HOLE, the conductance in pS, assuming 1 M KCl, for the different levels is as follows: N4-114 pS, N5-921 (16), N6-1247 (10), N6-20 ns-1214, N7-1555 pS, N8-1852 pS. Statistical errors are only a few percent but there probably are larger systematic errors. The conductance was calculated from snapshots spaced 10 ps apart and averaged. For this type of simplified calculation a scaling factor of five has been suggested to account approximately for the reduced mobility of ions in the narrow pores as well as other differences between the bulk solution and the solution in the narrow pores (Smart et al., 1997). This would give a conductance for N6 of ~250 pS, close to experimental levels (see Discussion).

We have also examined the properties of water in the Alm channel simulations. As a well-defined peptide-channel, alamethicin provides a relatively simple model for more complex ion channels. In previous simulations we have found that parallel helix bundles strongly orient water molecules inside the pore and that the diffusion coefficient of water is decreased relative to bulk values (for review, see Tieleman et al., 2001c). Such modified properties of water within ion channels are of importance in, e.g., parametrization Brownian dynamics methods for simulating ion channel current-voltage curves (Kuyucak et al., 2001). An interesting question that we can address using our N4-N8 models is whether water structure depends on the size of the channel. The molecular basis for the water ordering inside parallel helix bundles appears to be the strong electric field due to the combined helix dipoles of the parallel helices (Sansom et al., 1996). Thus, the effect is electrostatic in nature and might be expected to depend on relatively long-ranged interactions (see below). In Fig. 7 the water ordering is given for water in N5-N8. A similar ordering is found in all four channels. At the entrance and exit of the channel water molecules orient with their dipole moments antiparallel to the dipole of the helix bundle, whereas within the channel a strong ordering of water is observed, with the water dipoles antiparallel to the helical dipoles.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   Orientational order of water molecules along the pore axis for N5, N6, N7, and N8. The N4 channel does not form a continuous channel, but the orientation of the water molecules is similar. The peaks at ~1.5 and 5.5 nm are due to ordering of water molecules against the lipid interface (Tieleman et al., 1997) with the hydrogen atoms pointing toward the bilayer center and the oxygen atoms pointing toward the water layer.

Sensitivity analysis for N6

Based on the preceding analysis, two features suggest that the N6 model merits more detailed analysis. First, the bundle model exhibits unusual conformational stability on a 4-ns timescale, comparable with that of experimentally determined channel structures. Second, the calculated pore conductance (~250 pS after scaling by a factor 5) appears in good agreement with experimental data, although some caution is required due to the approximate nature of the conductance estimate. Therefore, we have explored the extent to which the behavior of the N6 model is influenced by the simulation protocol and also how it behaves in a simulation longer (by ~4×) than the mean passage time of an ion through the channel. A comparable analysis of sensitivity to simulation conditions has recently been performed for a K channel model (Capener and Sansom, 2002).

To this end we have performed additional simulations of the N6 model, using different algorithms to treat long-range electrostatic interactions (namely twin-range cutoffs, particle mesh Ewald (PME) and a reaction field (RF)and a simulation that includes 0.5 M salt solution in the bulk solution). It is known that the treatment of such interactions in simulations may have effects on the behavior of water and ions (e.g., Allen and Tildesley, 1987; Feller et al., 1996; Hess, 2002), peptides (e.g., Weber et al., 2000), and interfaces (e.g., Feller et al., 1996; Tobias et al., 1997), so given the complexity of the N6 system, we wished to see what the overall effect of changes in simulation algorithm were and whether they in any way modified the biological conclusions to be drawn from such simulations.

Comparison of the structural drift (i.e., RMSD versus time; Fig. 8) over 4 ns for the various N6 simulations reveals that this is lowest (~0. 18 nm after 4 ns) for the cutoff simulation but is only a little higher for the other simulation conditions (e.g.,~0.2 nm for the RF simulation). Furthermore, secondary structure analysis using DSSP (data not shown) did not reveal any significant differences between the various N6 simulations. Thus, the conformational stability of the N6 bundle seems to be relatively robust to the simulation protocol used, suggesting that this stability may be inherent in the initial model for the conductance level of the channel.

View larger version (21K): [in this window] [in a new window]   FIGURE 8   RMSD from the starting model for N6 with different simulation protocols. Data are running averages over 100 ps and are calculated from the backbone atom positions after fitting on all C- carbons.

The algorithm used to treat long-range electrostatics does appear to have a significant effect on the behavior of water within the Alm N6 pore. Fig. 9 compares the ordering of water under different conditions. Fig. 9 A shows the density profile of water, protein, lipid, and octane to aid orientation. In Fig. 9 B, the ordering of water in 4-ns simulations with cutoff, PME, and reaction field and a 2-ns simulation with PME that includes 0.5 M NaCl solution in the bulk phase are compared. The difference between the cutoff simulation and the other three is striking. The maximal degree of ordering of water within the pore is significantly lower in N6-PME, N6-RF, and N6-0.5M than in N6. The glutamate residues are located between 4 and 5 nm. In that region the ordering along the z axis is much less, due to the local influence of the glutamate ring. The differences between N6-PME and N6-0.5M are relatively small but consistent: water ordering in the presence of salt is somewhat less everywhere in the system, which is expected due to the shielding of long-range electrostatic interactions by salt.

View larger version (42K): [in this window] [in a new window]   FIGURE 9   Comparison of orientational order of water molecules along the pore axis using different electrostatics algorithms and environments. The density profile in A is given for comparison and to aid interpretation of the water orientation and is calculated from simulations PME-pr and PME-oct-pr (Table 1). (B) Water orientation in simulations N6, N6-PME, N6-RF, and N6-0.5M. (C) Water orientation in simulations with PME in the POPC bilayer and in an octane slab, with harmonic position restraints on all nonhydrogen atoms of alamethicin. PME2 in both cases stands for a system with an additional 4 nm of water added. In those cases, the x-scale extends to 11 nm but has been truncated after 7 nm in the graph. PME3-oct-pr is a large octane slab of ~18 × 18 nm rather than 6 × 6 nm. (D) Simulations comparing cutoffs and a reaction field correction in lipid and octane.

We have attempted to probe in more detail the origin of the difference in water ordering within the pore between different simulation protocols (Fig. 9, C and D). To dissect out effects on water ordering from possible effects on helix bundle structure, we carried out a set of simulations in which the channel atoms were harmonically restrained to their initial positions (with a force constant of 1000 kJ mol¹ nm² on nonhydrogen atoms), but the water was free to reorient (see Table 1).

First, we compare the ordering of water using PME (Fig. 9 C) between the same channel model in a lipid environment and in membrane mimetic octane slab in an attempt to isolate the contributions of the lipid headgroups from the contributions of the protein and bulk solvent. Also shown is the water orientation in the same two systems but with an additional 4-nm-think layer of water added to address the issue of possible artificial periodicity effects in PME (Hunenberger and McCammon, 1999; Weber et al., 2000). Finally, the water orientation in a system with a water layer of the same thickness but with nine times as much octane is shown. It appears lipid headgroups do not influence the orientation of water inside the channel much as this is dominated by the charges on the protein. Nonetheless, there are some differences (e.g., at 2.5 nm and between 4 and 4.5 nm) that might be caused by the long-range effects of the lipid headgroups. There also appears to be a somewhat larger effect of adding extra water in the case of octane. The zwitter-ionic headgroups are likely to contribute to shielding the interactions between images in PME, whereas octane has no charges. However, the effect of adding an additional 4-nm layer of water does not appear to be very strong, suggesting periodicity artifacts are not a major concern in these simulations. Increasing the size of the octane slab does not significantly alter the water orientation inside the channel, but does decrease the average orientation in bulk because only part of the bulk water is affected by the protein charges. When comparing simulations that use PME to simulations using cutoffs or reaction field (Fig. 9 D), water has a small net orientation everywhere in the system. This is consistent with what would be expected from a macrodipole oriented along the z axis in a dielectric medium. Without PME, this effect is not reproduced. In N6-0.5M this effect is not present either, because there the ions in the salt solution redistribute to compensate for the helical dipoles. When the protein is restrained, both cutoff and reaction field give similar results for water ordering, both rather different from what is obtained with PME (Fig. 9 D).

Overall, these results imply that PME or a comparable method is preferred for treating electrostatic interactions. They also imply that water orientation inside a narrow channel is mostly determined by the protein and only to a small extent by interactions with lipids and solvent.

Long N6 simulations

Having established that the N6 model appeared to be that in best agreement with the experimental conductance data, and that the stability of this model was not sensitive to the simulation protocol, we wished to explore the dynamic behavior of this on a longer timescale than that for mean passage of an ion. A simulation of 20 ns (~4 × the mean passage time) was therefore performed. In particular, we wanted to examine fluctuations in the bundle structure and pore radius profile that might occur on this larger timescale. We are interested in fluctuations in pore structure on timescales comparable with those of the permeation of ions, as these may play a role in the process of permeation, via phenomena such as stochastic resonance (Doering and Gadoua, 1992). A fuller understanding of such phenomena is likely to be of importance if more coarse-grained simulations of ion flux through channels are to accurately predict physiological data.

The structural drift of the N6 simulation (Fig. 10 A) (which was run independently for the 20-ns simulation from the earlier 4-ns simulation) showed a similarly low RMSD to that in the shorter simulation. Thus, the RMSD shows several plateaus but is less than 0.25 nm during most of the 20 ns. This gives us confidence that more detailed analysis of the longer timescale dynamics of the N6 bundle is worthwhile.

View larger version (22K): [in this window] [in a new window]   FIGURE 10   (A) N6-20ns RMSDs plus snapshots of the bundle at 2.5, 10, 17.5 ns. (B) Secondary structure of N6-20ns as calculated by DSSP.

Because the N6-20ns simulation is based on a model of the hexameric Alm bundle and not a "true" (e.g., crystallographic) structure, we also calculated the RMSD of every snapshot structure along the trajectory with respect to every other snapshot structure along the same trajectory (Fig. 12). A similar calculation can show the differences between each of the individual helices of the hexameric bundle by fitting each helix at each time to every other helix at every other time point (Figs. 11 and 12). Each helix has the lowest deviation when compared with itself but a slightly higher deviation when compared with other helices. In N6-20ns, helix two differs significantly from the other five helices. This is due to a change in secondary structure in the middle of helix two, where some residues do not have -helical geometry during most of the 20 ns. Although the / angles remain in the -helical region of the Ramachandran plot, the  dihedral angle distribution for Gly-12 is shifted toward 0° compared with the other  dihedrals (details not shown). Also in N6-20ns, the C terminus of one helix partially unfolds for a few nanoseconds and refolds into a stable -helix.

View larger version (32K): [in this window] [in a new window]   FIGURE 11   Summary of helix motions in N6-20ns (A-C). (A) Lateral tilt of the helices, as it were within the wall of the cylinder formed by the bundle. (B) Radial tilt, as it were perpendicular to the wall of the cylinder. (C) Distance between each helix axis and the central axis of the bundle, calculated as the average of all helix axes.

View larger version (148K): [in this window] [in a new window]   FIGURE 12   RMSD matrix for the N6-20ns simulation. The upper left and lower right half of each matrix contains different information. The scale on all four sides of a matrix is in nanoseconds. Each point in the matrix is a color-coded RMSD value for the backbone atoms fitted to the backbone atoms of the structure indicated on the axis. The upper left half shows the RMSD of the backbone of each helix with respect to all the other helices. The triangles are a comparison of a helix with itself at different times during the simulation; the squares are comparisons between different helices. The lower right half of each matrix is the RMSD for backbone atoms fitted on backbone atoms between two different times for the full bundle. The left color legend corresponds to the upper left half, the right to the lower right half.

The helix bundle geometry can also be followed by a simple geometrical description. For each helix, we take residues two to five at the C terminus and residues 17 to 20 at the N terminus and use these to define an axis. The average axis of all helices is the bundle axis. Then we define the distance between a helix axis and the central axis as the distance between their midpoints, the lateral tilt of a helix axis with respect to the "cylinder wall" formed by all helices, and the radial tilt of a helix axis away from the bundle. In Fig. 11, these properties are summarized for N6-20ns. The radial and lateral tilt varies between approximately 5 and +5°. The distance from the helix axes to a central axes changes significantly over time and differs per helix. Combined with the essential dynamics results below these results suggest that in this simulation motions of individual helices and within individual helices dominate the total dynamics. Significant rearrangement of helices is possible on a 20-ns time scale, but it is also clear that the tilt orientations and the distances are not sampled well within 20 ns. Thus, there are fluctuations on a timescale comparable with that of ion permeation that might be expected to influence the mean permeation energy profile experienced by that ion. However, we are unable to sample the longer timescale fluctuations that might be expected to contribute to, e.g., excess open channel noise (Sigworth, 1985).

We have attempted to assess more rigorously the situation with respect to sampling of conformational fluctuations in the Alm N6 helix bundle. A covariance analysis on the backbone coordinates (not taking into account two residues at either terminus) of N6-20ns simulation over 20 ns only provides limited information on possible collective motions of the helix bundle. Most of the motion seems to be intrahelical. The first two eigenvectors, which would show the most important collective motions, look like cosine functions. Using the method of Hess (2000), the cosine content estimate is 85% and 75%, respectively, for the first two eigenvectors. Hess (2000) has shown that this corresponds to mostly random diffusive motion on a flat potential surface. The first eigenvector consists mainly of motions of the helix termini and changes in the distance of the helices to the center. In comparison, the cosine content of the first two eigenvectors in simulations of the small globular proteins HPr and lysozyme in water is less than 50% in a 15-ns simulation (Hess, 2001). To estimate the degree of sampling, we calculated the overlap as defined in the Materials and Methods between 4-ns subsets of the 20-ns simulation (Table 2). This overlap is generally low, but comparable with the same calculation on 8-ns segments of 40-ns simulations of HPr and lysozyme, where the average overlap between any combination of 8-ns trajectory parts for HPr was 0.45 and 0.43 (for two different 40-ns simulations) and for lysozyme 0.44 and 0.48 for two different 40-ns simulations.

A covariance analysis on the individual helices can be used to estimate how much of the same phase space is sampled by the helices. Using the overlap measure described in Materials and Methods, the 15 possible pairs of helices in the N6-20ns six-helix bundle show an overlap of between 0.6 and 0.75. This indicates that the dynamics of the different helices are comparable, but on this timescale not yet identical. Thus, if there are any collective motions of an alamethicin helix bundle (e.g., one might imagine collective motions that distorted the cross-section of the pore from circular to elliptical, thus reducing its conductance) they are likely to occur on a timescale much longer than 20 ns.

To investigate possible conductance fluctuations we have also examined the minimal pore radius as a function of time for this simulation (Fig. 13). Although the average minimal radius of N6 is 0.22 nm, the fluctuations on a longer time scale are large, and in fact at least twice in 20 ns the minimal radius drops below 0.14 nm. Based on the rather large fluctuations in minimal pore radius and pore radius profiles, geometrical methods to link the structures from the simulations to channel conductance would need to incorporate the dynamic changes in the channel. For a narrow pore, perhaps such as N4, taking an average pore radius profile or an average structure might hide fluctuations that do allow passage of ions, even if the average radius is too small. Taking different structures and averaging over some estimate of the conduction based on these different structures would be a significant improvement.

View larger version (24K): [in this window] [in a new window]   FIGURE 13   (A) RMSD as a function of time for N6-20ns, calculated from the backbone atoms after fitting on C-s. (B) Minimal radius of the channels as function of time, calculated using the Connolly option in Hole, for (A) N6-20ns and (B) N4-N8.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

General context

It would be very useful to be able to model membrane proteins based on their sequence and some low-resolution structural data. Alamethicin is an attractive model system to test this process on because of its known monomeric structure, the presumably simple structure of bundles formed by alamethicin and the abundance of experimental data. Single-channel measurements in particular in theory provide excellent data to compare models with because they are extremely sensitive to changes in structure and probe a single molecule. A long-term goal of our work is to be able to calculate the conductance levels of alamethicin channels, predict their selectivity, and explain phenomena like rectification and pH-dependence of the selectivity from molecular models. This will require both accurate model building and higher-level simulations such as Brownian dynamics to make the connection to the observed electrophysiological properties. Understanding the flexibility of helices and the motions they undergo within helical aggregates is important to understand conformational changes, for example linked to proline residues, that might be involved in ion channel gating (Tieleman et al., 2001b) or transport through various transporter proteins. The motions of helices within a helical bundle on the timescale of our simulations are also relevant to the development of coarser models of protein mechanics. Finally, simulations of simple systems like alamethicin are useful to understand the effect of different simulation methods.

To what degree can we address these issues using the set of simulations in this paper? The structure of the individual alamethicin helices remains mostly stable in all simulations, although small fluctuations occur. It has been noted before that on a 1- to 100-ns time scale there is no difference between starting with an ideal helix model or a crystal structure model, but the environment (single helix in water, methanol, membrane, bundle) has some influence (Tieleman et al., 1999b; Tieleman et al., 2001a). This seems consistent and provides some confidence that simulations of other helical peptides give reasonable structures. The bundle structures themselves show significant changes over 4 ns and in one case over 20 ns. Interpreting these changes is more difficult. It seems a bundle of four alamethicin peptides is not stable, but it cannot be ruled out that this is due to at least to some extent to an inaccurate starting model. The N6 bundle appears to be the most stable. It is possible that this is simply the best model out of the five alamethicin channels, although there is no obvious reason for this based on the modeling procedure, which incorporated the same assumptions for all of N4-N8 (Breed et al., 1997), or on the simulation procedure that is also the same for all of N4-N8. The only conclusive tests of this would be to get an experimental structure of the channels, which is unlikely given the transient nature of the channels, or to reliably calculate a critical piece of experimental data, such as the conductance under specific conditions

Limitations of the models

There are several limitations to the models presented here. The models are built using a vacuum modeling procedure, simulated using the approximations present in the potential function and force field, and limited by the timescale we can reach with the simulations. At present they agree with available data, but part of the reason is that there are no experimental data that we can compare to with enough accuracy to improve the models. An accurate method of linking models to measured conductance and selectivity would provide a significant step toward developing better models. General improvements in simulation methods, computer speed, and force fields could also lead to improved models even if we do not have the experimental data to be certain of that.

An important simulation step concerns the charge state of the Glu-18 residues and merits some discussion. The parallel helix bundles have two electrostatic features that explain a possibly large shift in pKa values compared with the pKa of glutamate in water. First, the Glu residues are located near the C terminus of the helical peptides, experiencing the strong field caused by the helix dipoles that shifts the pKa of glutamate to higher values. This was also observed in models of the nicotinic acetylcholine receptor (Law et al., 2000). Second, once one Glu is charged, the electrostatic repulsion between the first charged Glu and a second one in close vicinity of the first one causes an additional shift in pKa. In fact, this combination of effects might be simpler to treat compared with the case of an Alm K18 mutant, in which the glutamates have been replaced by lysines. In that case, the two effects counteract each other to some extent, the helix field stabilizing the charged state of lysine, the lysine-lysine interactions destabilizing the charged state (Borisenko et al., 2000). However, even in the present case there remains some uncertainty. Although the continuum calculations suggest a significant shift of pKa values, these calculations have their own limitations, including the use of Poisson-Boltzmann approximations in a narrow channel, the choice of dielectric constants for the pore water and the protein, and the fact that they are based on a single model with a particular arrangement of charged groups. Simulations of an N6 bundle with a charge of 6 were unstable compared with simulations with a charge of 0 (Tieleman et al., 1999a), but this was without salt present. In simulations of Alm K18 octameric bundles, a model with a +8 charge without salt was unstable, but when 0.5 M or 1 M salt was included in the simulation the same model was structurally stable over 10 ns (D. P. Tieleman, M. S. P. Sansom, and G. A. Woolley. 2002, submitted manuscript). Thus, it is conceivable that more than one Glu-18 could be charged. Interestingly, with more than one Glu-18 charged in at least the wider channels, N7 and N8, their radii might have increased somewhat and better agreement of the calculated with the measured conductance might have been obtained. Structures from in particular the N6 simulations with different charge states and salt concentrations could form an interesting basis for a critical evaluation of the sensitivity of the continuum calculations.

Channel conductances

Channels in the lowest conductance state are impermeable to Ca²⁺ and Cl, whereas the larger channels do conduct these ions. Methods based on electrodiffusion theory in various forms have been qualitatively reasonably successful in reproducing current-voltage curves for simple channels, in particular LS3 (Dieckmann et al., 1999). However, even in that case only classes of structures could be considered, and a specific model could not be singled out based on its calculated conductance. Interestingly, the averaged structure of a molecular dynamics simulation gave the best agreement with the experimentally observed conductance levels. Monte Carlo simulations of ions passing through atomic-detailed models are a promising approach to link atomic structures to current-voltage curves and measured selectivities, but they too suffer from a lack of accuracy at this point, even in well-defined pores like OmpF porin (Im et al., 2000; Kuyucak et al., 2001). Equivalent cylinder models (Sansom, 1993b) and other geometry-based methods will be limited in accuracy but have the virtue of simplicity. The initial motivation for modeling alamethicin bundles with different numbers of helices was to attempt to link these structural models to the observed conductance levels in alamethicin channels (Fig. 1). Several sets of measurements are available from the literature. Hanke and Boheim (1980) measured levels in 1 M KCl of 19, 280, 1300, and 2700 pS for the first four conductance levels. Sansom and Mellor measured 88, 550, 1200, 2000, and 2700 pS in 0.5 KCl (Sansom, 1991), corresponding to ~190, 1100, 2400, 4000, and 5400 pS in 1 M KCl. The lowest level of 19 pS was not detected, but it seems reasonable the next three levels in the experiments of Hanke and Boheim correspond to the first three of Sansom. Vodyanoy et al. (1993) measured the effect of polymers on alamethicin channels and presented relative conductances (in 1 M NaCl) for three levels of alamethicin: 2.2, 4.5, and 7.2. These levels appear to be the same as the second, third, and fourth level of Sansom and Mellor. If the 1100 pS level is given a relative value of 2.2, 2400 pS would correspond to 4.8, and 4000 pS to 8. Similarly, if the 1300 pS level of Hanke and Boheim would correspond to 2.2, the relative conductance of the 2700 pS level would be 4.6. In conductance measurements of covalently linked alamethicin dimers, three different levels were found, in which the third was 1 nS in 1 M KCl. This level was suggested to correspond to the hexamer. The first had a very low conductance and the second, presumably corresponding to a pentamer formed by 2.5 dimer had a relatively low lifetime compared with the third level (You et al., 1996). The combined data could be interpreted as the existence of a low conductance state under certain conditions, formed by a tetramer; the next level was always seen corresponding to a pentamer with an approximate conductance of 550 pS (the first level of Sansom and Mellor, and of Vodyanoy et al., 1993); the next level corresponding to a hexamer with a conductance of ~1000 to 1300 pS, and the next level with a conductance of ~2400 to 2700 pS. Higher levels have also been described (Fig. 1).

The calculated "raw" conductance of the different channel models averaged over 4 ns (20 ns for N6-20ns) were: N4-114 pS, N5-921 (16), N6-1247 (10), N6-20ns-1214, N7-1555 pS, N8-1852 pS. These values are based on the assumption of bulk properties for 1 M KCl in the channels, including bulk diffusion coefficient and bulk resistance. Calibration calculations of HOLE suggest that these values should be scaled by approximately a factor of five, mainly based on lower diffusion coefficients of water and ions inside a narrow channel. There is now considerable evidence that the diffusion coefficients for ions in narrow channels is a factor 5 to 10 lower than in bulk (for review, see Tieleman et al., 2001c). Thus, scaling the calculated values by a factor of 5 would yield ~20, 180, 250, 310, 370 pS for N4, N5, N6, N7, N8, respectively. Regardless of a scaling factor, the increase in conductance predicted by this simple method is nearly linear from N5 to N8, which does not agree with the experimental measurements. However, it is currently not possible to determine conclusively whether this is a problem with the models or with the simple method to estimate the conductance. It seems likely that for the larger channels diffusion based methods become more accurate, which suggest that the values for N7 and N8 are on the low side relative to each other and to N6. For N4 and N5, the average radius of the pore (~0.2 nm) approaches the limit of applicability of diffusion-based approaches (Tieleman et al., 2001c). In N4, the actual minimal radius is practically always smaller than an ion radius, but this is not taken into account by the conductance estimate as this simply takes the volume of a slice through the pore.

Electrostatics

A technical point that we have tried to address for alamethicin concerns three different methods to treat electrostatic interactions. This has been an important topic of research in the past with an extensive literature considering liquids (Allen and Tildesley, 1987), interfaces (Feller et al., 1996; Tobias et al., 1997), DNA, and proteins in solution (Levy and Gallicchio, 1998; Wang et al., 2001). A cutoff, even at 1.8 nm, causes artificial structuring of water and extra noise in the simulations (Hess, 2002). The orientation of water in the narrow channel with a high electric field due to the charge distribution on the helices is the most sensitive property calculated in this paper. Both N6-PME and N6-RF give approximately the same water orientation profile, which differs significantly from the profile calculated with the cutoff method. This difference agrees with the results in pure bulk water. The remaining difference between a reaction field and PME are likely due to the incorrect approximation the reaction field approach makes in the inhomogeneous membrane systems and the neglect of long-range electrostatic effects outside the cutoff range due to the lipid headgroups and outlying parts of the protein. However, PME calculates interactions between all periodic images. In aqueous solution, this can lead to artificial stability of helices under certain fairly extreme conditions (Hunenberger and McCammon, 1999; Weber et al., 2000). In the low-dielectric environment present in our simulations, less extreme conditions might already have a significant influence on the dynamics of the proteins. In test simulations on parallel helix bundles (models of Influenza A M2 and others), the RMSD of the bundle was much lower with PME than with cutoffs (Forrest and Sansom, unpublished results). In other test calculations that compare the properties of phosphatidylcholine liquid crystalline bilayers when simulated with PME or cutoffs, significant changes in areas per lipid, and other structural parameters occur in simulations of the order of 25 ns (unpublished data). It seems that PME is preferable over the RF treatment and the twin-range cutoffs. A more systematic characterization and exploration of other alternatives such as lattice summation in two dimensions (Yeh and Berkowitz, 1999) and nonperiodic methods (Warshel and Papazyan, 1998) are desirable for simulations of inhomogeneous lipid membrane systems. An interesting observation is the relatively small effect of salt on water orientation in the N6 channel. Additional simulations will more fully explore the distribution of salt and its effect on water and channel structure in synthetic alamethicin mutants (D. P. Tieleman, M. S. P. Sansom, and G. A. Woolley. 2002, submitted manuscript).

Time scale and correlated motions

Ion permeation in alamethicin channels occurs on a 5-ns time scale. Because at 20 ns our longest simulation is significantly more than that, we hoped to observe concerted motions of the bundles as a whole, and hypothesized that the high-frequency noise observed in the open state of the channels might be linked to fluctuations in the structure on a 50- to 500-ns time scale. It appears that on at least a 20-ns time scale, motions of individual helices dominate the total dynamics of a bundle. Compared with water-soluble proteins, there is no convergence even in up to 20 ns in the first principal components of the motion of N6-20ns. A similar behavior is observed in a 17-ns simulation of a model of a pentameric pore formed by the M2 peptide from the nicotinic acetylcholine receptor (Law et al., 2002). This might be due to longer intrinsic timescales of helix motion in membrane proteins, but it is also possible that the potential energy landscape for a peptide channel is quite flat and dominated by motions of single helices. If we take into account its function, which is just a passive pore that conducts ions on a 5- to 10-ns time scale, this might be a logical explanation. If collective motions of the bundle are not responsible for noise, smaller changes in structure on a 100-ns time scale that have a significant effect on the electrostatic potential could be an alternative source.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Although alamethicin is one of the simplest ion channel models, it remains challenging to determine its structure and link the observed conductance levels reliably to molecular models. Nonetheless, it is becoming increasingly feasible to thoroughly explore molecular models; the time scales of simulations are approaching the time scale of ion permeation in simple channels, and different simulation parameters and conditions can be tested to assess the sensitivity of results to assumptions made in the simulations. In particular, water ordering in the narrow alamethicin channels depends strongly on including long-range electrostatic interactions. The key findings of this study are: we provide simulation-based data that suggest which experimentally observed conductance levels correspond to which peptide aggregation state; N4 appears too small to form a stable channel; collective motions on a 20-ns time scale do not appear to play an important role in ion conduction through alamethicin channels; and our model for N6 is as stable in a 20-ns simulation as simulations based on crystal structures. The model not only shows a detailed structure of a peptide channel and lends confidence to general properties calculated from the simulations with relevance for other membrane proteins, such as water orientation, but it also opens up the possibility of combining mode building with recent developments aimed at describing ion permeation. In the near future an accurate link between modeling and electrophysiology data, including selectivity, rectification and conductance levels should become possible. This will be invaluable for understanding many different ion channels as well as for designing new ones.

There remain several significant differences between our simulation models and a real channel. For instance, we do not know for certain whether the structure of our model is correct, a real channel is bathed in salt solution and interacts with permeating ions, it is possible that the presence of a transmembrane potential difference in the electrophysiology measurements is important for the motions of the bundle, and it has been shown that the likelihood of observing a specific conductance state depends on the type of lipid. Future work will address some of these issues for alamethicin and other membrane proteins.