INTRODUCTION

Top Abstract Introduction Conclusions References

Ion channels regulate the ionic concentration across a membrane. They are responsible for signal transduction and transmission, and are one of the most important constituents of living cells. Physiological processes involving ion channels range from very elementary sensory systems to the massively interconnected information-processing network of the mammalian brain. Hundreds of functionally distinct ion channels have been identified, and the amino acid sequences of ion channel proteins have been established. The commonly accepted structure for a great variety of ion channels involves an aqueous pore. The protein that makes up the pore walls forms an envelope, undergoing favorable interactions that compensate for the partial dehydration experienced by an ion as it traverses the channel. From an electrostatic point of view, the protein may be considered as a medium, which screens the image charge of the ion, because of the low permittivity of the surrounding membrane environment. Unfortunately, despite the obvious curiosity and speculation concerning the nature of their three-dimensional structure and functional mechanisms, important membrane proteins, such as AChR (Imoto et al., 1988; Unwin, 1993; Beroukhim and Unwin, 1995), are too large and complicated to be investigated reliably by computer simulation at this time. Accordingly, the focus to date has mostly been on much simpler systems (Åqvist and Warshel, 1989; Cooper et al., 1988; Jordan, 1987; Pullman, 1991; Roux and Karplus, 1991).

To understand the function of ion channels, Lear et al. (1988) adopted a "minimalist" approach, starting with the synthesis of small peptides. The synthetic peptides they employed consisted of simple repetitive sequences, but with some of the structural features thought to be important in the natural channel-forming proteins. This minimalist approach offers the advantage that molecular models used to explain observed functional properties have fewer structural possibilities than those generated with natural sequences. Moreover, the simpler sequence peptides can be expected to have somewhat simpler functional properties. The more complex functions of the natural proteins might then be built up by introducing modifications and additions designed to test specific hypotheses concerning the role of individual residues. These model systems are also more amenable to the current generation of molecular simulations, which can then be carried out in much the same environment as that used for the laboratory counterparts.

One of the first minimalist channel peptides designed by Lear et al. (1988) had the sequence Ac-(LSLLLSL)₃-CONH₂, which is 21 residues in length, so that it can span the hydrophobic portion of a typical lipid bilayer. Leucine (L) was chosen for the apolar face of the helix because of its high hydrophobicity and helix-forming propensity. Serine (S) was chosen for the hydrophilic face of the helix because it is polar, but without a net charge. This peptide was found to assemble as a tetramer and behave as a proton-conducting channel (Lear et al., 1988; Chung et al., 1992; Akerfeldt et al., 1993; Kienker et al., 1994; Kienker and Lear, 1995). The simplicity of the heptad repeat in the design sequence minimizes issues related to the protein folding problem.

The -helical coiled-coil structural motif was originally proposed by Crick (1953). In the Crick model, the axes of the -helices are inclined toward each other by an angle of ~20°. Furthermore, the -helices deform slightly, so that they remain in contact and wind around each other, like a two-stranded rope. Coiled-coils are often found as a structural motif in native proteins, such as muscle regulatory protein, DNA-binding protein, and cGMP-dependent protein kinases. The secondary structure is entirely -helix, and the helices are packed according to complementary, so-called knob-into-holes packing. The coiled-coil motif provides an ideal model for studying forces responsible for protein folding and protein-protein recognition. It has recently attracted intensive research in both experimental and computational studies (DeLano and Brünger, 1994).

Molecular dynamics (MD) simulation has proved itself to be a useful tool for understanding the properties and functions of complex systems, ranging from self-assembled monolayers (Hautman and Klein, 1989) to micelles (Watanabe and Klein, 1989) and membranes (Berendsen, 1996). The hydrophobic environment is vital to the structure and function of ion channel proteins; these proteins have also been studied via simulation (Åqvist and Warshel, 1989; Jordan, 1987; Woolf and Roux, 1994; Elber et al., 1995; Roux et al., 1995). To date, much of the focus has been on gramicidin, and important results have been obtained. However, in the case of larger channel proteins, simplifications have often been necessary. For example, in some of the work of Sansom et al. (1995), the -helices composing the channel were constrained to reduce the computational overhead. Although such a procedure can give information on the behavior of the pore water, it clearly cannot give insight into the dynamical comportment of the channel.

In this work we have combined some of the recently introduced novel simulation methodologies (Martyna et al., 1992, 1996) with state-of-art parallel computing to probe the behavior of a minimalist synthetic proton channel (Lear et al., 1988) in the presence of a membrane-mimetic environment. Specifically, the latter has been taken into account by placing a four-helix peptide bundle in the octane region of a water/octane simulation box. The evolution of the four-helix bundle was then followed during an MD simulation.

Anticipating our results, we will see that once formed, the close-packed coiled-coil channel persisted for the duration of the MD simulation, which spans several nanoseconds. The serine side chains, lining the hydrophilic inner surface of the channel, were found to be bound mostly through pore H₂O. The latter formed a network, through which an excess proton should easily pass via a Grötthus-type mechanism (Sagnella et al., 1996). From an analysis of the MD trajectory, we have identified a slow damped collective mode of the peptide bundle, with a period of ~300 ps. Two relaxation times of pore H₂O were observed: one around 3 ps, and another that is much slower.

The paper is organized as follows. First we present a brief outline of the simulation procedures. Then we discuss the system set-up and the assembly of the model channel. This section is followed by an analysis of structural and dynamical properties of the four-helix bundle. The article ends with a brief conclusion and discussion of possible avenues for future research.

	MULTIPLE-TIME-STEP MOLECULAR DYNAMICS

To get optimal performance in both system size and time scale, we have employed the MD scheme that is based on explicit reversible integrators, combined with the multiple time step method of Tuckerman et al. (1992). The van der Waals and electrostatic interactions that compose the interatomic forces were each divided into short- and long-range components. The short-range interactions were truncated at 7.0 Å for the electrostatic terms and at 2.0 , where  is the usual Lennard-Jones parameter, for van der Waals interactions. The short-range interactions were calculated every 1.5 fs. However, with the chosen cutoffs, they contain only about half of all of the nonbonding interactions. The major part of the nonbonding interaction, the long-range part, fluctuates slowly and was therefore calculated only every 3 fs. The intramolecular interactions, including stretching, bending, and dihedral angles, are calculated every 0.3 fs.

Periodic boundary conditions were used with an overall cutoff of 2.5  for the van der Waals interactions. A standard correction (Allen and Tildesley, 1987) was utilized for the neglected long-range contributions of van der Waals interactions. The Ewald method was employed to take into account the long-range electrostatic interactions (Allen and Tildesley, 1987; Ciccotti et al., 1987), with a 10-Å real-space truncation, 10 Å¹ as the cutoff in reciprocal space, and a value of  = 0.3 for the weight of the Gaussian damping factor.

The simulation was carried out at nominal room temperature, 300 K, with the temperature controlled by a Nosé-Hoover chain thermostat (Martyna et al., 1992). We used a single Nosé-Hoover chain of length 3, and the frequency factor of the chain was chosen to be 2 ps¹. The thermostat equations of motion yield continuous dynamics that generate a canonical distribution (Martyna et al., 1996). In the part of the trajectory where the averaged properties are calculated, no constraints were applied to the system.

The new simulation code we employed was developed at the University of Pennsylvania (Moore and Klein, 1997) and parallelized using the Message Passing Interface (MPI). The bulk of the simulation was done at the IBM T. J. Watson Laboratory, Yorktown Heights. The code takes ~1000 s/ps of trajectory, using 16 processors of an IBM SP1.

	SIMULATION SYSTEM

The Ac-(LSLLLSL)₃-CONH₂ -helical peptide was first set up in an ideal right-handed -helical form using INSIGHT (Biosym Technologies, San Diego, CA). Then the peptide was minimized via an adapted-basis Newton-Raphson algorithm (ABNR), using an optimized version of the CHARMM program and the CHARMM 19 united atom parameter set (Brooks et al., 1983), with all backbone atoms fixed, so that the atoms on the side chain could relax. Next, the minimized helix was duplicated and assembled as a parallel tetramer with fourfold symmetry. The interhelix separation (~10 Å) was chosen to avoid any nonphysical repulsions. For example, the leucine side chains should likely not be intertwined. The energy of the four-helix bundle was then minimized by ABNR, with all backbone atoms fixed, to eliminate any possible bad contacts between the helices. The resulting structure (see Fig. 1) provided the starting configuration for the ion channel simulation.

View larger version (95K): [in this window] [in a new window]   FIGURE 1   The central portion of the initial set-up of the octane-water-peptide simulation system. The four helices, shown in a ribbon representation, have the parallel tetramer arrangement proposed by Lear et al. (1988), which has fourfold symmetry (see text).

In the laboratory, the channel protein is solvated in a lipid bilayer. Unfortunately, at the present time, the lipid bilayer membrane diphytanoylphosphatidylcholine (diPhy-PC) used in laboratory experiments (Lear et al., 1988) is too demanding in CPU time to be included in a long MD simulation. However, one of the most important properties of the lipid bilayer, among all others, is a well-defined hydrophilic/hydrophobic interface. Here liquid hydrocarbon has been used to provide the channel peptide with an appropriate membrane-mimetic environment. For this purpose, we have chosen the hydrocarbon n-octane (C₈H₁₈), because it is a relatively simple liquid that forms a stable interface with H₂O at room temperature (see Fig. 2). The same basic idea has been used before in other ion channel simulations (Roux and Kaplus, 1991; Sagnella and Voth, 1996).

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Linear density profile of octane () and water (·····) along the channel direction, taken from the MD run. Here L(z) has units of m/Å, where m is in atomic mass units. The horizontal line indicates the bulk H2O value. The sharp interface between the hydrophilic and hydrophobic environments was maintained even after the 6-ns-long simulation run (see text).

Individual simulation systems of both H₂O and C₈H₁₈ were prepared using constant volume (NVT) and constant pressure (NPT) MD simulations, at room temperature and atmospheric pressure, each of which ran for over 500 ps. Then the C₈H₁₈ box was cut into a slab (45.45 × 47.45 Å²) with the same height (34.68 Å) as the ideal -helix of Fig. 1. The lateral dimensions of the simulation box were chosen in such a way that the four-helix bundle was reasonably far from any of its images. Next, a hole was cut in the middle of the octane box and the four-helix bundle was introduced. Then two slabs of H₂O molecules, each 6 Å thick, were taken from the H₂O simulation and added to both ends of the peptide bundle. The H₂O layer thickness was chosen based on the known characteristics of the H₂O/C₈H₁₈ interface (see Fig. 2). A string of H₂O containing 15 molecules was cut from the equilibrium box and inserted in the middle of the bundle. The resulting system, with 222 C₈H₁₈ molecules, 957 H₂O molecules, and the four -helices, with dimensions of 45.45 × 47.45 × 46.68 Å³, composed the starting configuration for this study.

To eliminate the tension between different components and to solvate the four-helix bundle, we initially constrained the bundle and carried out an NVE-MD run on the H₂O and C₈H₁₈ subsystem for 200 ps. The well-known TIP3P model (Jorgensen et al., 1983) was used for both the bulk and pore H₂O. However, because of the multiple-time-step integrator used for the equations of motion, the O-H bond length was not constrained. The parameters used for C₈H₁₈ were those recommended by Siepmann et al. (1993) for a fully flexible united atom model of methyl and methylene groups. The topology and parameters for the peptide were taken from CHARMM 19 (Brooks et al., 1983). This version of the parameter set is based on a united atom model for the peptide residues and their interactions, and hence is consistent with the molecular model used for the hydrophobic C₈H₁₈. All hydrogens, except the polar hydrogen of serine, the amide, and water, were absorbed in heavy atoms, which leads to a significant reduction in the number of interaction sites. The resulting simulation system contained, in total, 5415 united atoms.

	FOUR-HELIX BUNDLE

A 200-ps NVE-MD run was started from the set-up described above. This was followed by a 400-ps run under NVT conditions. During both runs, the velocities were reassigned every 100 steps to hold the temperature at 300 K. After 600 ps of temperature scaling, a NVT-MD run was carried out, which lasted more than 5 ns.

The behavior of the system first appeared to stabilize after 1000 ps, most likely because the initial setup was far from the preferred conformation of the peptide bundle and optimal side-chain packing. During the NVT run, the helices tilted, twisted, and associated to form a channel. The helices first evolved from the initial idealized tetrameric structure proposed by Lear et al. (1988) to a state with two associated dimers and overall twofold symmetry (see Fig. 3). The latter state persisted for over 3 ns. Then the helical bundle evolved back to a tetrameric structure (see Fig. 4). After that, no significant changes were observed over a further 2 ns of simulation. The final tetramer structure is consistent with the open channel structure proposed by Lear et al. (1988), in which the helices form a left-hand coiled-coil (see Fig. 5). Besides this, the helices are close packed and support the necessary pore structure with a hydrophilic inner and a hydrophobic outer surface (see Fig. 6). In the final structure, the inner, polar side chains of serine are solvated by the pore H₂O molecules, and the outer, hydrophobic side chains of the leucine residues are solvated in C₈H₁₈. Very importantly, despite significant excursions of each helix around average positions, the bundle shows no tendency to dissociate over the time of this simulation.

View larger version (120K): [in this window] [in a new window]   FIGURE 3   An instantaneous configuration of the backbone of the four-helix peptide bundle at ~2 ns, as viewed from the N terminus. The 24 oxygen atoms of the serine side chains are drawn (in green) with van der Waals radii. The carbonyl oxygens are drawn in red, and the amide nitrogen and hydrogens are in blue and grey, respectively. All other atoms are omitted for clarity. The bundle has the "dimer-of-dimers" structure (see text).

View larger version (112K): [in this window] [in a new window]   FIGURE 4   An instantaneous configuration of the backbone of the four-helix peptide bundle at ~4 ns, viewed from the N terminus. The 24 oxygen atoms of the serine side chains are drawn (in green) with van der Waals radii. The carbonyl oxygens are drawn in red, and the amide nitrogen and hydrogens are in blue and grey, respectively. All other atoms are omitted for clarity. The bundle has a tetramer structure and a larger pore than the configuration shown in Fig. 3 (see text).

View larger version (58K): [in this window] [in a new window]   FIGURE 5   An instantaneous configuration of the four-helix bundle at ~5 ns viewed from the side, with the N terminus at the top. The helices, which are drawn in a ribbon representation, have formed a left-handed coiled-coil structure.

View larger version (128K): [in this window] [in a new window]   FIGURE 6   An instantaneous configuration illustrating the compact nature of the four-helix peptide bundle at 5 ns, as viewed from the C terminus. The atoms of the pore structure are drawn with van der Waals radii: red, oxygen; black, carbon; blue, nitrogen; white, hydrogen.

Throughout the MD trajectory, the H₂O molecules are pumped into the channel pore and others are pumped out. In the dimer and tetramer states, respectively, there are ~30 and 40 H₂O molecules remaining in the channel. Fig. 7 actually shows the time evolution of the number of pore H₂O molecules in the central portion of the channel only. Figs. 8 and 9 show instantaneous configurations of pore H₂O taken from the simulation.

View larger version (19K): [in this window] [in a new window]   FIGURE 7   Number of pore H2O molecules, N, associated with the "dimer-of-dimers" (- - -) and tetramer () structures as a function of time. Pore H2O molecules were defined as those whose oxygen atoms lie within 7.5 Å of the center of mass of the bundle.

View larger version (74K): [in this window] [in a new window]   FIGURE 8   An instantaneous configuration of the pore H2O molecules in the "dimer-of-dimers" state of channel. The pore H2O molecules, oxygen (yellow), and hydrogen (grey) are drawn with van der Waals radii.

View larger version (76K): [in this window] [in a new window]   FIGURE 9   An instantaneous configuration of the pore H2O molecules in the tetramer state of the channel. The color code is the same as in Fig. 8. The chains of hydrogen-bonded H2O molecules should be noted.

	STRUCTURAL AND DYNAMICAL PROPERTIES

During the simulation, each peptide is rotating, stretching, and bending, but the -helices remained intact, with the hydrogen bonds along each back-bone well maintained. In addition, the bundle undergoes certain modes of collective motion, namely breathing, twisting, and stretching.

An ideal -helix has 3.6 residues per turn. On the other hand, the channel helix has a periodicity of 7. Thus it has 0.2 residues less than required for two complete turns. The helix therefore acquires a twist of ~10° to have all of the polar side chains pointing to the pore. The net twist of the bundle was calculated by monitoring the angle between the total channel dipole and the dipole of each helix and taking the average over the MD trajectory. The result, 17.5 ± 4.8°, is consistent with the ideal value of 20° proposed by Crick (1953).

To explore further the structural properties, we have calculated the dipole moment of each helix as well as that of the whole bundle (see Fig. 10). At equilibrium, each peptide has a dipole moment of ~15.0 e × Å, where e is the charge of the electron, which is essentially equivalent to a charge of ±0.50e located at the N and C termini of the -helix, respectively. This is consistent with the conventional picture of the -helix and supports the assumption underlying the electrodiffusion model of Kienker et al. (1994).

View larger version (34K): [in this window] [in a new window]   FIGURE 10   The total dipole moment of the peptide bundle, M(t), as a function of time for the tetramer structure. The unit of dipole moment is electron charge (e) times Ångstrom. The horizontal line represents the total value of the dipole moment when each helix has a charge of ± 0.50e at the N and C termini, respectively, with a helix length of 30.0 Å.

Further analysis of the dynamics was carried out by calculating the inertia tensor of the bundle. The inertia tensor has two roughly equal components in the plane perpendicular to the channel direction (see Fig. 11). These components and the one along the channel direction show periodic oscillations (see Fig. 11). We have also calculated the power spectrum of each component of the inertia tensor. There is a broad peak around 0.003 ps¹, a frequency corresponding to breathing of the helix bundle.

View larger version (25K): [in this window] [in a new window]   FIGURE 11   Evolution of the inertia tensor, I(t), of the four-helix bundle for the tetramer structure projected along the principal axes in units of m × Å2, where m is in atomic mass units. The two components in the plane perpendicular to the channel (upper two panels) are very similar. All components, including those along the channel axis (lower panel), show evidence of the slow collective motions of the bundle.

The serine side chains, which form the hydrophilic surface of the channel peptide, lie mostly in the plane perpendicular to the channel direction. The protons of the serine hydroxyl groups are generally pointed toward the backbone carbonyl at the next turn of the -helix. This arrangement helps to stabilize the helical structure and exposes the serine hydroxyl oxygen atoms to the pore. The serine hydroxyl oxygen atoms on neighboring helices are linked together, mostly through pore H₂O, via hydrogen bonds (see Figs. 8 and Fig. 9). In this way, the helices form a stable channel, with the serine polar side chains adopting a layered structure that acts as a bottleneck in the channel. However, the collective motion of the helices allows the pore H₂O to be transferred between the layers of serines and eventually pass through the channel.

We have calculated the diffusion coefficient of the H₂O molecules in the system. Three kinds of behavior of H₂O have been observed (see Fig. 12). The pore H₂O molecules have a diffusion coefficient that is about one-third of the bulk value, whereas the bound H₂O molecules stay linked to the serine side chains and hold the helical bundle together. However, each H₂O changes its role among the three kinds of behavior during the course of the MD trajectory.

View larger version (17K): [in this window] [in a new window]   FIGURE 12   Calculated mean square displacements, [z(t)  z(0)]2, of three typical H2O molecules. The diffusion coefficient of the mobile pore H2O (- - -) is about one-third of the bulk H2O () value. The pore-bound H2O (·····) is immobile on the time scale shown.

The number of hydrogen bonds made by each oxygen of the pore H₂O fluctuates around the value 1.4, with a mean square deviation of ±0.2, which is ~30% lower than the value of 1.9, associated with bulk TIP3P water (Jorgensen et al., 1983). We also calculate the dipole moment of the pore H₂O. The averaged dipole moment of a pore H₂O molecule is zero perpendicular to the channel direction, but 0.35e × Å in the channel direction. This suggests that most of the pore H₂O molecules are arranged in such a way that one of their hydrogens mostly points along the channel direction opposing the dipole of the bundle.

We have calculated the dipole moment auto correlation function for the H₂O molecules in the system, because this quantity is accessible to NMR experiment. For the pore H₂O, contrary to a previous report (Breed et al., 1996), we found that the dipole relaxation cannot be fitted to a single exponential decay. However, we found that a two-exponential regression is able to describe the relaxation of the H₂O dipole (see Fig. 13):

(1)

The fitted values reveal that the first term has essentially the same relaxation time as the bulk H₂O. However, the second term has a relaxation time ~100 times larger. The same calculation has been repeated for the H₂O in the cap area, and a similar behavior was observed. The difference between the cap and pore H₂O is that the latter is dominated by the second, slower relaxation term, whereas the former is dominated by the first term. An intuitive picture that may explain the double-exponential regression is that the pore H₂O can rotate like bulk H₂O about the channel axis direction, but rotation about the axis perpendicular to the channel direction is constrained. In the pore, the large field in the channel direction resulting from the dipole of helices gives rise to the large value of B. In the bulk H₂O, there is no preferred direction, and B = 0, i.e., the second term vanishes. The behavior of the cap H₂O molecules is intermediate between those of the bulk and the pore.

View larger version (15K): [in this window] [in a new window]   FIGURE 13   Rotational autocorrelation function, µ(t) · µ(0), of cap () and pore H2O (- - -) as a function of time. The correlation functions can be fitted by two exponential functions, one of which has the typical bulk H2O relaxation time, and the other of which has a relaxation time about 100 times longer. The circles and squares are a double-exponential fit to the simulation data (see text).

	CONCLUSIONS

Top Abstract Introduction Conclusions References

We have performed a constant-volume MD simulation, spanning more than 5 ns, on a synthetic ion channel consisting of a four-helix peptide bundle. The necessary hydrophobic/hydrophilic environment has been taken into account so that the channel can be followed dynamically during the simulation. We have shown that the four -helices assembled in the presence of a membrane-like environment form a close-packed pore structure. The present synthetic channel first spontaneously adopted an associated dimer conformation, which then evolved to a left-handed coiled-coil tetramer structure. The assembled four-helix bundle undergoes certain well-defined oscillatory low-frequency motions during the MD simulation. The two states observed may possibly correspond to the two conducting states of the channel noted by Lear et al. (1988). The hydrophilic serine side chains are solvated by the pore H₂O. The latter form a network through which a proton should easily be able to pass (Sagnella et al., 1996; Pomes and Roux, 1996). The hydrophobic leucine side chains are solvated in C₈H₁₈. The simulation also shows that the pore H₂O has a solvation structure different from that of the bulk H₂O. Specifically, we found evidence for two distinct reorientation times of pore H₂O. Our results also suggest that the CHARMM 19 parameter set provides at least a reasonable initial description of this synthetic ion channel. A more detailed comparison of the associated dimer state and the tetramer state, and dynamical information, including the nature and the damping of the collective mode, the effect of introducing a proton, and the dynamical properties of pore H₂O, along with a parallel method utilizing a different potential parameter set, are currently under investigation.