INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Integral membrane proteins can form transbilayer pores, as ion channels (Unwin, 1989; Hille, 1992) and other transport proteins (Cowan et al., 1992; Engel et al., 1994; Walz et al., 1995). Many classes of ion channels are formed from a bundle of parallel -helices that surround a central water-filled pore and form a pathway through which the ions travel. Examples of this type of channel include the mechanosensitive channel MscL (Chang et al., 1998), the M2 protein from influenza A virus (Forrest et al., 1998), and the nicotinic acetylcholine receptor (Unwin, 1993, 1995; Sankararamakrishnan et al., 1996). The bacterial K⁺ channel KcsA (Doyle et al., 1998) is more complex, the pore being formed by bundles of antiparallel helices into which is inserted a selectivity filter. Depending on the size of the pore and the charge distribution of the protein, these channels are charge selective. Wider pores (~0.9 nm) are less selective, while narrow pores (~0.4 nm) are highly charge selective and can also discriminate between isovalent ions. The mechanisms of these selectivities are poorly understood, except at a qualitative level. Few channel propertiesparticularly those of intrachannel waterare understood in any depth. Because these channels are important in molecular transport, and molecular transport is important in many fields such as drug design and delivery, they are worthy of deeper investigation.

Synthetic ion channels were first investigated in depth in the 1980s (Oiki et al., 1987, 1988; Lear et al., 1988) because of their minimalist yet functionally representative structure. Such channels provide pertinent information in both experimental and theoretical settings, while eliminating some of the more complicated features of natural channels. Properties such as charge selectivity, dipole orientation, and translational/rotational mobility of water molecules are among those studied either theoretically or experimentally. Among the first synthetic ion channels to be studied are (LSLLLSL)₃NH₂, known as LS2, and (LSSLLSL)₃NH₂, known as LS3 (Lear et al., 1988). These systems are long enough (21 residues each) to span the hydrophobic region of a typical lipid bilayer and are able to aggregate to form channels because of their amphiphilicity. Leucine (L) has hydrophobic properties and a high propensity for helix formation, while serine (S) is hydrophilic with no net charge. Experimentally, the leucine-serine channels exhibit a voltage-dependent gating mechanism (Kienker et al., 1994; Lear et al., 1997). However, even in the absence of an electric field, stable open pore structures are observed for periods of at least a millisecond (Kienker et al., 1994). Because molecular dynamics (MD) simulations of such systems run on the order of nanoseconds, any open pore, once formed, should be stable for the duration of the simulation.

These systems were chosen for the present study because of their simple structures and monovalent charge selectivities. In the current model, the individual helices are packed together in a parallel fashion, packing all of their C termini on the same side of the bilayer. Previously, leucine-serine systems have been studied by MD in vacuo (Mitton and Sansom, 1996; Lear et al., 1988; DeGrado and Lear, 1990) and with an octane slab mimicking a bilayer membrane (Zhong et al., 1998a,b). Tetrameric LS2 and hexameric LS3 systems were found to most closely reproduce experimental observations (Åkerfeldt et al., 1993; Mitton and Sansom, 1996; Dieckmann et al., 1999). In this paper, these systems are solvated for the first time in a palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer with water at each cap and in the pores, to more realistically represent the environment of the channels. This study extends from simulations of OmpF (Tieleman and Berendsen, 1998) in a POPE bilayer, of gramicidin A (Woolf and Roux, 1994), and of alamethicin in a POPC bilayer (Tieleman et al., 1999).

The LS3 peptide forms monovalent cation-selective ion channels with a single-channel conductance (G = 70 pS in 0.5 M KCl) (Lear et al., 1988, 1994; DeGrado and Lear, 1990; Åkerfeldt et al., 1992, 1993). This is consistent with hexameric and possibly even pentameric channels and is represented by the LS3 hexamer in this study, shown in Fig. 1. The LS2 peptide prefers bundles of four helices, is proton selective, has dimensions similar to those of the influenza A M2 protein, and is studied here with the LS2 tetramer.

View larger version (128K): [in this window] [in a new window]   FIGURE 1   A snapshot of the LS3 system after 2 ns of simulation. The protein is shown in ribbon format, and lipid carbonyl oxygen atoms are shown in spacefill mode. The extents of each region are marked by the arrows. w, water; p, phospholipid headgroup; h, hydrophobic region.

With our computer simulations, we seek to show whether the two channels are stable in this realistic membrane environment and to find a detailed explanation of the basis of their different ion selectivities. Thus we have examined the bundle structure and the behavior of the pore-lining residues and pore waters during the simulation.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Models of four helix bundles were generated by restrained in vacuo MD, using a simulated annealing (SA) protocol as previously described (e.g., by Kerr et al., 1994) and were reported by Mitton and Sansom (1996). Briefly, this involves generation of idealized -helical templates, which were made into tetramers or hexamers while ensuring that serine side chains were pore-lining. Helices were capped at both termini by NH₂ groups. These were used in a two-stage SA-MD method incorporating distance restraints to maintain -helical backbone conformations and to hold the helices in four- and six-helix bundle conformations. Each run of this procedure yielded an ensemble of 25 structures, from which structures with high four- and sixfold symmetry were selected as the starting point of extended MD simulations in a bilayer/water environment. In the case of the LS2 system, noncrystallographic symmetry restraints were applied to generate highly symmetrical models.

The LS helix models were then embedded in a preequilibrated lipid bilayer consisting of 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphatidylcholine (POPC) as used in previous simulations (Forrest and Sansom, 1999; Tieleman and Berendsen, 1998; Tieleman et al., 1999). T_c for POPC is 5°C (i.e., 268 K). Thus at 300 K, POPC is in the liquid-crystal phase as opposed to the gel phase, making it a suitable representation of a bilayer environment (Seelig and Waespe-Sarcevic, 1978). Cylindrical holes were made in the center of the bilayer by removing lipids whose phosphorus atoms fell within 1.55 nm of the central axis of the cylinder, followed by running short MD simulations with a radially acting repulsive force to drive atoms out of the cylinder into the bilayer (Tieleman et al., 1999). The helix bundle was then placed within the hole and then energy minimized. The entire system was solvated with a minimum of 30 SPC waters per lipid. Each system was once more energy minimized, followed by a 25-ps MD equilibration stage. LS2 consists of 3517 water molecules (none initially in the pore), 103 lipids, and four helices, for a total of 16,655 atoms. The LS3 system has 4443 water molecules (75 waters initially in the pore region), 102 lipids, and six helices, for a total of 19,737 atoms. MD was run for a full 4 ns on each system. In both systems, waters entered the pore to form a well-defined column after 100 ps.

The MD simulations were carried out using periodic boundary and constant NPT conditions. The simulation box measured x = 5.88 nm (± 0.038 over the course of the simulation), y = 5.53 nm (± 0.05), and z = 7.48 nm (± 0.11), where z is parallel to the bilayer normal. A constant pressure of 1 bar was applied independently in all three directions, using a coupling constant of _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 at 300 K, using a coupling constant _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.8 nm for electrostatic interactions. Although some arguments can be made for using noncutoff methods for electrostatic forces (i.e., Ewald summation), previous work by Tieleman and Berendsen (1998), Tieleman et al. (1999), and Forrest and Sansom (1999) with this protocol has agreed well with experiment. The time step was 2 fs, using LINCS (Hess et al., 1997) to constrain bond lengths.

The MD simulations were carried out on 2-processor and 10-processor, 195 MHz R10000 Origin 2000s and took ~8 days per processor per nanosecond simulation. Simulations and analysis were carried out using the GROMACS (Berendsen et al., 1995) suite (http://rugmd0.chem.rug.nl/^~gmx/gmx.html) with GROMOS87 parameters. Only polar H-atoms were represented explicitly. Initial models were generated using Xplor (Brunger, 1992). Structures were examined using Quanta (Biosym/MSI) and Rasmol, and diagrams drawn using MolScript (Kraulis, 1991).

Pore radius profiles were measured using HOLE (Smart et al., 1997), and average profiles over the whole simulation were used to predict the approximate ionic conductance (G) of the pore of LS3:
where s is the width of the slab at z, A is the area of the cross section, and r' is the effective pore radius. The values z_low and z_high are the z coordinates at the ends of the channel, according to the average of the terminal residues (Smart et al., 1999). Thus the last two terms describe the access resistances of each mouth of the pore. The resistivity (z) is measured using
where D_bulk = 5 × 10⁹ m² s¹ and is given by the maximum value in Fig. 14 B and D(z) is the diffusion constant given by the simulation. The values for G were calculated using _bulk = 0.13 m (Mitton and Sansom, 1996), equivalent to the resistivity of 0.5 M KCl.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Progress of the simulations

Stability of protein structures over the time course of simulations can be evaluated by analyzing the deviation of the protein structure with respect to initial structure. Fig. 2 shows the C- root mean square deviation (RMSD) for each bundle over the duration of each simulation. LS3 exhibits an increase in its RMSD to ~0.2 nm, where it remains for the last 1000 ps, while LS2 fluctuates around 0.11 nm before rising to 0.17 nm, where it remains for the last 1000 ps of the simulation. The fluctuations of the RMSD for both systems are of comparable height and occur at similar frequencies. The overall RMSD is less than 0.2 nm, which is comparable to those seen for simulations starting from the x-ray structures of membrane proteins such as OmpF (Tieleman and Berendsen, 1998), of KcsA (Shrivastava and Sansom, 1999), and of FhuA (Sansom, 1999, unpublished results). Thus by this, albeit relatively crude, criterion the helix bundles had a stability comparable to that of crystallographically determined membrane protein structures.

View larger version (28K): [in this window] [in a new window]   FIGURE 2   Root mean squared deviation (RMSD) of the C atoms of the helix bundles during the simulation.

The magnitudes of deviations from a time-averaged structure can also provide useful information about internal fluctuations in structure. The RMS fluctuations of the C-'s as a function of residue number are plotted in Fig. 3. The termini fluctuate much more than the middle of the peptides. This is expected because the termini are not connected to the lipid bilayer, except through hydrogen bonding. The fluctuations for most of the length of the bundles are low (~0.07 nm). However, the fluctuations for LS3 (Fig. 3 B) may be slightly higher than those for LS2 (Fig. 3 A), partly because of the larger number of water molecules penetrating the pore (see below). Neither the RMSD nor the RMS fluctuations are significantly greater than those observed during simulation of single helices within a lipid bilayer (Forrest and Sansom, 1999) or during simulation of single helices in aqueous solution (Bodkin and Goodfellow, 1995), thus again indicating the stability of these bundles.

View larger version (28K): [in this window] [in a new window]   FIGURE 3   Root mean squared fluctuations per residue for (top) LS2 and (bottom) LS3. Vertical bars indicate the extents of the helices.

Bundle fluctuations

Visual examination of the C- backbones of LS2, superimposed every 250 ps (Fig. 4), provides the reader with a picture of the true stability of the bundle. Interhelix distances plotted as a function of time (Fig. 5) show how the pores become somewhat wider. For LS2 (Fig. 5 A), two interhelix distances (between helices 1-2 and 3-4) are significantly larger than the other two. This may give some indication of a "dimer of dimers" structure, as seen previously by Zhong et al. (1998a) (see below). Some LS3 interhelix distances also increased (Fig. 5 B), but over a larger range and with no distinctive pattern emerging.

View larger version (64K): [in this window] [in a new window]   FIGURE 4   Superimposed C traces taken at intervals of 250 ps for (A) LS2 and (B) LS3 viewed down the z axis and for (C) LS2 and (D) LS3 viewed in the x-y plane, with the N termini at the top. H1-H6 denote the helix number.

View larger version (35K): [in this window] [in a new window]   FIGURE 5   Distances between adjacent helices as a function of time for (A) LS2 and (B) LS3.

Fig. 6 shows the crossing angles for adjacent helix pairs in each of the systems, as described by Chothia et al. (1981). The helix crossing angles for the tetramer (Fig. 6 A) are much larger (10-30°) than for the hexamer (Fig. 6 B) (0-20°), indicating a more pronounced coiled-coil nature, as might be expected, and as observed previously (Mitton and Sansom, 1996). Helix tilt values with respect to the bilayer normal show similar behavior, with the values for LS2 (Fig. 7 A) being slightly greater and spread over a wider range, and the values for LS3 stabilizing to ~10° by the end of the simulation (Fig. 7 B).

View larger version (30K): [in this window] [in a new window]   FIGURE 6   Helix crossing angles () for adjacent helix pairs as a function of time for (A) LS2 and (B) LS3. H1-H6 denote the number of the helix as defined in Fig. 4.

View larger version (29K): [in this window] [in a new window]   FIGURE 7   Angle of tilt of helix with respect to bilayer normal as a function of time for (A) LS2 and (B) LS3.

Interactions between helices can also give an insight into bundle stability. To this end, patterns of hydrogen bond formation between pairs of adjacent helices have been analyzed. In considering the results of such analysis, it should be remembered that side-chain motions are relatively slow (Tieleman et al., 1999) and H-bonds are long lasting, so side-chain interactions are incompletely sampled. The criteria used for assigning an H-bond were an H-donor-acceptor angle of <60° and a donor-acceptor distance of <0.35 nm. Fig. 8 shows the fraction of time each residue spends forming a H-bond with each helix. This fraction (f_H) was calculated as follows. If n_H = the total number of time steps during which a given residue was observed to form a H-bond; n_T = the total number of time steps in the simulation; and q = the number of observed types of H-bond formed by the residue in question during the simulation, then f_H = n_H/(n_T · q). Thus, if a given residue makes a single type of H-bond throughout the entire simulation, f_H = 1. If, for example, a serine hydroxyl donates an H-bond to two different water oxygens during the course of a simulation but always remains H-bonded to one or the other oxygen, f_H = 0.5. 

View larger version (26K): [in this window] [in a new window]   FIGURE 8   Fraction of H-bonding, fH (as defined in the text) for donor atoms from every residue in helices 1-4 or 1-6 to acceptor atoms in the other helices. Bars shaded in black indicate H-bonds formed by serine residues, and unshaded bars indicate H-bonds formed by leucine residues. (A) Interhelical H-bonding pattern for LS2. (B) Interhelical H-bonding pattern for LS3.

Close inspection of Fig. 8 reveals a lack of symmetry for the hydrogen bonding between the helices and residues. This is expected, as we did not study H-bonds between each possible residue pair. To get a symmetrical graph, one would have to calculate the H-bonds between each residue and every other residue. This would result in as many plots as there are residues (88 for LS2!), which would be difficult to analyze. In this paper, H-bonds are examined from each helix in the bundle to all other residues in the bundle. Thus Ser2 on helix A may form many H-bonds to Leu1 and Leu3 on helix B; however, Ser2 on helix B may not form reciprocating bonds to Leu1 and Leu3 on helix A, instead forming H-bonds with residues on helix C or intrahelical H-bonds. This information is not available from our plots.

Significant interactions are noted between each helix and its immediate neighbors, particularly involving serine side chains. In fact, nearly all of the interhelix hydrogen bonds occur through the serine residues. In LS2 (Fig. 8 A) the fraction of H-bonding is slightly greater between helices 1 and 4, and again between helices 2 and 3. Coupled with the slightly longer interhelix distances noted above, a tendency toward a "dimer of dimers" structure seems to exist. However, both the distance and the fraction of H-bonding differences are small. LS3 exhibits no such preference in its H-bonding (Fig. 8 B). For both structures, intrahelix H-bonds were omitted.

Pore and water behavior

During the course of the simulation, a well-defined column of water forms in each of the pores. It is well known that the way in which water behaves inside the lumen of pores is very different from its behavior in its bulk arrangement (Jakobsson and Chiu, 1987; Breed et al., 1996; Sansom et al., 1996; Hartnig et al., 1998; Engels et al., 1995). In addition, the arrangement of the water molecules inside the pore affects ion selectivity and permeation rates. To this end, the dynamics and orientation of the intrapore water have been analyzed in detail, as have certain features of the lining of the pore.

Figs. 9 and 10 show snapshots of both LS3 (A) and LS2 (B) after 4 ns of simulation. It is clear that in both cases, the hydroxyl groups of the serines point their H atoms toward the N termini. This creates a field across the bilayer that is influential in determining the direction in which ions flow across the channel, reflecting their behavior as inward rectifiers, i.e., a lower resistance to ion flow in one direction than in the other (Kienker et al., 1994; Woolley et al., 1997; Dieckmann et al., 1999). Such behavior was also observed previously in in vacuo simulations (Mitton and Sansom, 1996). However, in contrast to the results of Zhong et al. (1998a), the H atoms of the Ser residues are not all pointing toward the backbone carbonyls of the next turn of the -helix.

View larger version (32K): [in this window] [in a new window]   FIGURE 9   The (A) LS3 and (B) LS2 helix bundles, where backbone atoms are drawn in ribbon mode and serine side-chain atoms are shown explicitly, viewed down the z axis. The N-terminal group is facing the reader. The hydrogens (in light gray) of the Ser-OH groups are facing the N terminus. Oxygen atoms are colored black.

View larger version (38K): [in this window] [in a new window]   FIGURE 10   The (A) LS3 and (B) LS2 channels, where backbone atoms are drawn as ribbons and serine side-chain atoms are shown explicitly, viewed in the x-y plane. One helix has been removed from each bundle for clarity. The N termini are at the top. The Ser-OH groups can be seen to point toward the N terminus. Oxygen atoms are colored black; hydrogens are light gray.

Pore radius profiles (Fig. 11) were calculated at t = 0 ns and t = 4 ns for both systems, using HOLE (Smart et al., 1997). Essentially, HOLE measures the diameter of the pore by using Monte Carlo simulated annealing to maximize the radius of a sphere at different lengths along the pore. The radius of each sphere is then recorded as the pore radius at that point. The radius of each pore increases ~0.05-0.10 nm during the course of the simulation. At the end of 4 ns, the LS2 pore is wider at the C termini than at the N termini. LS2 (Fig. 11 A) has an average radius of 0.22 ± 0.05 nmroom for about one or two water molecules across the porewith its narrowest point having a radius of 0.15 nm. LS3 (Fig. 11 B) has an average radius of 0.36 ± 0.05 nmor about two or three water molecules across the porewith its narrowest region ~0.2 nm in radius. These results are slightly larger than the restricted in vacuo simulations of Mitton and Sansom (1996), where the pore radii for LS2 and LS3 are 0.2 ± 0.03 nm and 0.34 ± 0.09 nm, respectively.

View larger version (23K): [in this window] [in a new window]   FIGURE 11   Pore profiles at t = 0 ns (solid lines) and t = 4 ns (dashed lines) for (A) LS2 and (B) LS3. The C termini (ca. z = 2.3 nm) open out more than the N termini (ca. z = 5.6 nm) during the simulation. The pore of LS2 appears to lose its periodicity. The extents of the lipid bilayer are marked with dotted lines.

The z-coordinate trajectories of several waters located in the pores at the end of the simulation are plotted in Fig. 12 as a function of time. Most waters that entered the LS2 pore remained in the pore for the length of the simulation, which explains the increase in pore radius during the simulation. The waters that were observed leaving the pore left on the same side on which they entered (Fig. 12 A, Water 409). It is clear that waters inside the pore of LS2 (Fig. 12 A) remained in the same position for long periods, nearly frozen. In particular, Fig. 12 A shows water 719 moving farther into the pore in jumps of ~0.2 nm at a time (see t = 1400 ps and t = 1600). This is similar to the diffusion seen in a five-staved poly-Ala -helix bundle (Smith and Sansom, 1998). Essentially, when this type of rare movement is observed, the average diffusion coefficient approaches zero, as seen at z = 2.2-5.2 nm in Fig. 13the graph of water diffusion coefficients versus the average z coordinate. Simulations of LS2 in a bilayer mimetic environment (Zhong et al., 1998b) showed three types of water diffusionbulk, pore-bound, and pore-mobile. In this work, however, the pore waters only appear to demonstrate the "pore-bound" behavior, whereby diffusion is considerably restricted. This behavior may enable the formation of a "water wire" network, which is important in the transport of protons via the Grotthüs mechanism (Schmitt and Voth, 1998).

View larger version (29K): [in this window] [in a new window]   FIGURE 12   Trajectories of sample water molecules along the z axis over time. (A) Waters 409, 1020, and 719, from top to bottom, from the LS2 simulation. (B) Water 9244 from the LS3 simulation. (right side) Final positions of water molecules 409, 1020, and 719 from LS2. In each case the extents of the lipid bilayer are marked with dashed lines.

View larger version (15K): [in this window] [in a new window]   FIGURE 13   Diffusion coefficients (D) of water calculated along the length of the pore for (A) LS2 and (B) LS3. In LS2, water in the pore is frozen, while for LS3, the diffusion constant is greatly reduced compared to bulk, but is not completely frozen. The extents of the lipid bilayer are marked with dashed lines.

The water in the pore of LS3 (Fig. 12 B) also showed reduced movement along the z axis while in the pore. Although water 9244 moved completely through the LS3 pore during the 4 ns of simulation, it moved much more slowly than it does in bulk water. This observation is confirmed by cross-referencing with the plot of diffusion coefficient versus z coordinate (Fig. 13 B). The diffusion coefficient remains above zero in LS3, at ~1/10 that of bulk water, in contrast to the values of zero observed in the pore region of LS2 (Fig. 12 A). The pore radius and diffusion coefficient results indicate a slightly smaller pore in this model than in previous simulations of LS3 within an octane environment (Zhong et al., 1998b). This, along with observed lipid-to-amino acid H-bonding (not shown), suggests that one effect of the presence of lipid headgroups is to further stabilize the pore.

To gain a better understanding of the interactions between the pore water and the helices, the degree of H-bonding between each residue and water was studied (Fig. 14). The majority of H-bonding for the LS2 channel (Fig. 14 A) clearly resides with the middle serines. This suggests a strong and long-lasting interaction between specific water molecules and the protein and provides further evidence (in addition to the diffusion coefficients) for the water being essentially "frozen" in the pore. The serine side chains at the termini of the helices show smaller H-bonding fractions to water, reflecting interactions with lipid headgroups and greater numbers of different waters. Thus the number of possible of hydrogen bonds (q) has increased, causing a decrease in f_H. The LS3 channel exhibits a larger capacity for H-bonding in its C-terminal serines (Fig. 14 B). However, the H-bonding in the pore is not nearly as strong as that of LS2. This indicates that the water in the LS3 pore is free and able to reorganize itself. The serines in the LS3 pore are thus free to form H-bonds with many more waters, again resulting in an increase in q and decreased f_H values. The H-bonding of Ser3, at the N terminus of helix 6 of LS3, appears to be considerably more prolific than for other residues in the helices. This is due to consistent H-bonding of both hydroxyl and backbone atoms to a single, unmoving water molecule during the majority of the simulation. The low number of possible H-bonds (q) results in a larger H-bonding fraction, as is also seen to a lesser extent for Ser3 of helix 1. 

View larger version (35K): [in this window] [in a new window]   FIGURE 14   Fraction of hydrogen bonds between each residue and solvent water for (A) each helix of H to H4 in LS2 and (B) each helix of H1 to H6 in LS3. For LS2, the central serines do the majority of the H-bonding.

With ~55 water molecules in the middle 3 nm of the pore and a total of 24 polar serine residues in the bundle, LS2 can and does form a large number of hydrogen bonds to the available water. The water molecules become aligned antiparallel to the helix dipole and have an average dipole moment of ~2 Debyes projected along the pore (data not shown). This should be compared to a dipole moment of 2.3 Debyes for a single simple point charge (SPC) water. However, the constriction of the pore (Fig. 11 A) and the hydrogen bonds formed between the water and the serine hydroxyl groups (Fig. 14 A) limit the mobility of these waters, freezing them inside the pore. This corresponds to an existing theory about why LS2 is proton selective (Mitton and Sansom, 1996). The pore size, coupled with the hydrogen bonding between the serines and waters, restricts the waters from passing through the pore. Instead, one or more water "wires" are formed along the interior of the pore, and hydrogens pass along this wire in a Grottüs fashion. In this theory, the proton transfer rate is limited by the rate at which the water can reorient itself (Zhong et al., 1998a; Pomes and Roux, 1998).

LS3 has a very different water ordering, starting with 111 water molecules in the inner 3 nm of the pore. The diffusion within the LS3 pore is not zero, although it is much reduced compared with that of bulk liquid water (Jakobsson and Chiu, 1987; Ghadiri et al., 1994; Breed et al., 1996; Tieleman and Berendsen, 1998; Smart et al., 1997). The larger diffusion rate should be sufficient to enable the diffusive motion of ions through the pore (Mitton and Sansom, 1996). In addition, the larger pore radius (Fig. 11 B) of the LS3 helix bundle suggests that monovalent cations other than protons may be able to pass through the pore. In particular, K⁺ ions, the ionic radii of which are similar to the van der Waals radii of water, would be able to pass through while remaining partially hydrated throughout (Roux and Karplus, 1993; Jakobsson and Chiu, 1987; Smith and Sansom, 1998). If the pore waters were "frozen" in the same manner as in LS2, then they might be expected to present an energetic barrier to ion permeation, despite this large pore radius. However, the diffusion rate is greater than zero, and thus creates no such barrier.

Given the average radius of the pore over time (Fig. 11 B) and the diffusion constants of water (Fig. 14 B), it is possible to predict the conductance of the LS3 channel by the method of Smart et al. (1997, 1999). Assuming a bulk diffusion of 5 × 10⁹ m² s¹ (using the GROMACS value for bulk water) and the resistivity of bulk 0.5 M KCl solution, _bulk = 0.16 m (Smart et al., 1997), we get a conductance of G = 64 pS. This corresponds to the experimental values of 70 pS for equivalent conditions (Lear et al., 1988). Such strong agreement suggests that the model of LS3 is a very good representation of the true system. The absence of diffusion data in the pore region (Fig. 14 A) and the inapplicability of this simple theory to H⁺ conductance mean that similar calculations have not been carried out for the LS2 helix bundle.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Simulations have been carried out on four- and six-helix bundles of the synthetic peptides LS2 and LS3, respectively, embedded in explicit lipid (POPC) bilayers. The results of the simulations are compared with the experimental properties of LS2 and LS3 and with the results from previous simulations of these systems within different environments. Importantly, the stability of the channels on these time scales, even in the absence of an applied field, indicates that these are valid representations of the true structures and dynamics of these synthetic ion channels. Both the LS2 tetramer and the LS3 hexamer showed low RMSD and RMS fluctuation values over 4 ns. In these simulations we use a starting model in which we assume that the helices are packed around a central water-filled pore. It may be that in doing so we are simulating the behavior of the "open" state of the channel rather than its "closed" state. Conceivably, one might model the closed state by packing the helices in a more compact, less symmetrical fashion, e.g., similar to the helix packing observed in bacteriorhodopsin (Edholm et al., 1995).

The behavior of the water in the pores of each helix bundle reflects the observed conductance properties of channels formed by synthetic peptides of the same sequence. In the pore of the six-helix bundle model of LS3, the waters are observed to move more slowly than in bulk water. However, they are sufficiently mobile to allow the diffusion of monovalent cations though the pore, in line with experimental results. Analysis of the bundle structure indicates that the pore dimensions are large enough to incorporate a diffusing ion. Furthermore, predicted conductance of K⁺ ions through the pore results in values in agreement with experimental evidence.

Electrophysiological experiments on the LS2 peptide show that it forms proton-conducting channels. The current simulations agree with this prediction, as the water within the pore is effectively "frozen" and exhibits very little diffusion. Such rigidity suggests the formation of a "water wire" along which a proton might transfer in "Grotthüs" fashion. Furthermore, the dimensions of the pore are small enough to confirm the selectivity of LS2 channels against larger ions.

The current work differs from previous simulations by the more realistic treatment of the environment surrounding the proteins. The need for interhelix restraints as previously incorporated in in vacuo simulations (DeGrado and Lear, 1990) has been eliminated by the use of an atomistic bilayer. Furthermore, in contrast to the work of Zhong et al. in their simulations of LS2 and LS3 in an octane slab (Zhong et al., 1998a,b), the use of POPC lipids in this case includes a representation of the charged headgroups found in membranes. However, these simulations are highly dependent on the starting structurethis is particularly apparent for the four-helix bundle model of LS2. It was noted that a previous in vacuo model of LS2 had weak packing at one of the helix-helix interfaces. Starting from this model resulted in helix unpacking during a bilayer simulation, indicating a sensitivity to poorly constructed models. Therefore, an in vacuo model with better helix-helix packing was generated for the current study.

It is important to remember that these simulations are still approximations of the true systems. Ions typically take 10-100 ns to move through a pore, so that this is unlikely to be reproduced in these 4-ns simulations. Although a recent simulation of the initial stages of peptide folding in water reached the µs time scale (Duan and Kollman, 1998), such simulations are not currently within the capabilities of most laboratories. Another approximation is the use of a cutoff of 1.8 nm during the calculation of long-range forces. This could be improved by using techniques such as Ewald summation or the particle-particle particle-mesh (PPPM) method. However, previous simulations of OmpF (Tieleman and Berendsen, 1998) (which contains considerable numbers of charged residues) and of alamethicin (Tieleman et al., 1999) in lipid bilayers, using our adopted protocol, have agreed with experimental evidence. Furthermore, the contribution of electrostatics to the interactions of the leucine-serine proteins are likely to be small because of the nature of the amino acids. Although these approximations serve to qualify our simulations, they nonetheless appear to reasonably model the system.

This paper has summarized molecular dynamics simulations conducted on synthetic ion channels in an explicit bilayertwo systems of more than 15,000 atoms each. The results agree with observed experimental results for the conductances of the channels. Future studies of these systems will include an explicit quantum mechanical treatment (Schmitt and Voth, 1998) of the proton transport mechanisms in LS2. The study of such synthetic channels will also eventually lead to greater understanding of simple, naturally occurring channels such as M2 from influenza A and of more complex structures such as KcsA.