- Email: [email protected]
Molecular dynamics study of structure and gating of low molecular weight ion channels
Parallel Computing 26 (2000) 965±976
www.elsevier.com/locate/parco
Molecular dynamics study of structure and gating of low molecular weight ion channels Dennis M. Newns a,*, Qingfeng Zhong b, Preston B. Moore b, T. Husslein b, Pratap Pattnaik a, Michael L. Klein b a
b
IBM, Thomas J. Watson Research Center, Yorktown Heights, NY 10598, USA Center for Molecular Modeling and Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323, USA Received 12 July 1999; accepted 26 October 1999
Abstract We implement molecular dynamics (MD) simulations on low molecular weight alpha helixbased functional synthetic and native ion channels. The synthetic channels are the LS2 proton channel and the LS3 voltage-gated channel. The simulation manifests key features of the channels such as the coiled-coil structure of the alpha-helix bundle and the continuous aqueous pore. By implementing simulations with and without an applied voltage, we develop a hypothesis as to the voltage-gating mechanism. The native channel is the M2 proton channel in the in¯uenza A virus, which plays an essential role in the infection process. This channel is pH gated via protonation of one or more imidazole rings in the H37 residues. Simulation of the neutral channel reveals a coiled-coil structure whose pore is penetrated by water, but not threaded by a water column. By means of simulations with dierent numbers of charged H37 residues, we demonstrate a possible gating action via opening up the channel to a continuous water column, and also provide support for the alternative proton relay gating mechanism. Ó 2000 Published by Elsevier Science B.V. All rights reserved. Keywords: Ion channel; M2; Proton channel; Molecular dynamics
1. Introduction Molecular dynamics (MD) simulation has reached a stage of development, where it can be realistically applied in environments as complex as protein/water and *
Corresponding author. E-mail address: [email protected] (D.M. Newns).
0167-8191/00/$ - see front matter Ó 2000 Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 1 9 1 ( 0 0 ) 0 0 0 2 1 - 1
966
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
protein/lipid systems. Current computer power places limits on the size of systems and/or the length of time for which systems can be studied. The protein folding problem, involving real times of order ms, poses a major challenge in terms of number of time steps required. However, there are problems of major signi®cance in biochemistry which have time constants of order ns and can be readily addressed with presently available computer capabilities, especially exploiting parallelism. Such problems can also serve as a test of the reliability of MD simulations of large biomolecules. Transmembrane proteins constitute 30% of a typical genome, and among these ion channels are recognized as having extensive biological importance; the central role, e.g., of the K channel in nerve conduction is well known, but ion channels are increasingly recognized as playing key roles in many other areas. The ion channel is a transmembrane protein whose role is, with some degree of selectivity, to transport an ion across the impermeable lipid membrane when subject to a stimulus, which may be, e.g., a voltage across the membrane, or a change in pH. The action of this stimulus in opening the channel to transport is termed `gating'. A typical ion channel structure is a peptide bundle (or a continuous protein) which in its transmembrane section is organized into a bundle of alpha helices. The center of the bundle constitutes the pore for ion transport; to what extent a water column pervades the pore is an issue in the ®eld which is highly relevant to the nature of ion channel action. The peptide sections external to the membrane, in the aqueous medium, may not be coiled. The earliest ion channel studies by MD [1±7] by Roux and others involved a slightly atypical channel, gramicidin, which is a beta barrel rather than an alphahelix bundle. Sansom's [8] group studied seven-helix bundle ion channels, and the behavior of water within various models of ion channels [9]. Sansom's work on the in¯uenza A virus M2 channel [10] provided an important basis of comparison for part of the present work. This group has also carried out an extensive study of the alamethicin channel [11,12]. In the work to be described below, Zhong et al. and Husslein and coworkers [13±17] have studied both functional synthetic channels and the M2 channel, with special emphasis on understanding of the gating mechanism. 2. Objectives The MD technique is intrinsically capable of locating the stable state of the macromolecular system provided (a) sucient simulation time is allowed and (b) no deep local free energy minima are encountered. We have found that in several nanoseconds ion channel conformations such hydrated, coiled-coil bundles can be stably established. However, a result such as bundle dissociation, which requires a longer timescale than nanoseconds to occur fully, can only be established in a more speculative manner within the framework of nanosecond simulations. In this paper, we shall describe MD simulations of two very dierent types of channels, selected for their relatively low molecular weight and the availability of
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
967
extensive experimental data on their properties, including information on their gating behavior. The ®rst systems are a pair of functional synthetic ion channels invented by Lear, Wasserman and DeGrado [18], aided by extensive modeling by Dieckmann et al. [19]. The synthetic channels, both built from 21 amino acid peptide sequences, are termed LS2 and LS3. The LS2 peptide consists of the sequence Ac± LSLLLSL3 ±CONH2 , and LS3 has the sequence Ac± LSSLLSL3 ±CONH2 (where L leucine and S serine). In the presence of a lipid layer, the LSN peptides are believed [18] to self-assemble into a transmembrane bundle of alpha helices, which experimentally [18±23] acts as a functional ion channel. The architecture [18] of the sequence is designed so that, viewed down its axis, each helix has an angular sector h which is hydrophilic, consisting of serine residues, and which is expected to orient towards the pore, while the remaining, larger, angular sector 2p ÿ h is hydrophopic, consisting of leucine residues. Thus, when the hydrophilic sectors face the pore, the bundle is expected to have a hydrophobic coating, enabling it to embed into the lipid bilayer, and a hydrophilic pore, capable of maintaining a transmembrane water column. In a naive model of the channel bundle as a symmetric ring of n touching cylinders with hydrophilic sector inwards, then n
2p 2 L S ; pÿh LÿS
1
where L and S are the number of leucines and serines, respectively. This is only a semiquantitative model, but it predicts that the larger the fraction of serine, the larger is the number n of peptides in the channel. For the LS2 peptide, Eq. (1) estimates n 4:7, and LS2 is believed to form [18] a tetrameric channel. Consistent with the small size of the pore in this case, LS2 is proton selective. For serine-rich LS3, Eq. (1) predicts the large value n 14; but in fact this channel is believed to form a hexameric or heptameric bundle. Experimentally it can transport larger ions such as Clÿ , and is voltage-gated [23]. The excessive hydrophilic sector h in LS3 relative to pore perimeter raises the possibility of signi®cant interhelix hydrogen bonding due to the overlap of OH groups from dierent chains. Molecular dynamics can be used to investigate a number of issues in the LSN channel area. Key issues are, whether the hypothesized bundles are stable, whether their architecture indeed includes the expected feature of pore-lining OH groups, whether there is a water column and how large is it, whether there is interhelix hydrogen bonding in LS3, and how does the voltage-gating of the LS3 channel work. Our second set of studies is of a native channel, the M2 channel of the in¯uenza A virus [24]. At two stages in the infection process of a cell by the virus, the M2 channel plays the role of a proton pump and is crucial to successful infection. In fact, there is an in¯uenza drug, amantadine, which is believed to operate by blocking the M2 channel [25,26]. The gating of the channel is caused by lowering the pH below, approximately pH 5.8. A key simpli®cation in the study of the M2 channel is the discovery by Du and Ashley [25] of a truncated 25 amino acid M2 sequence SSDPLVVAASIIGILHLILWILDRL, consisting of a transmembrane section with small extensions, which nevertheless gates like a complete M2 channel. The
968
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
functional truncated Du±Ashley sequence is ideal for MD simulation. Like the LS2 sequence, the Du±Ashley M2 peptide is believed to self-assemble into a transmembrane alpha-helix bundle in the presence of the lipid bilayer, the assembly being tetrameric [27]. In many ways, however, the M2 contrasts with the synthetic LSN type of system. First, in the M2, in contrast with LSN, there are few hydrophilic groups lining the pore. Secondly, the M2 is believed to contain a well-de®ned gating structure consisting of four histidine residues, one from each peptide in the bundle. The four histidine residues form a channel block towards the C-terminus of the peptides. There are two main hypotheses as to how the pH-controlled gating of the channel may occur. Both hypotheses are based on the gating pH 5.8 lying close to the pKa for protonation of the histidine imidazole ring [28]. In the multiple ionization hypothesis, two or possibly more histidines are protonated below pH 5.8, and the resulting inter-peptide chain electrostatic repulsion (which is largely unscreened due to the absence of a water column at this point) forces apart the bundle at the pore blocking point, leading to the formation of a water column and the consequent opening of the channel to proton transport [29]. The initial version of this multiple ionization hypothesis, due to Sansom [10], assumed ionization of all four histidine groups. However, according to the results of our simulations which are unconstrained [15], the quadruply ionized channel is unstable, while we ®nd that double ionization is sucient to open the channel, without endangering its stability. A second hypothesis for M2 gating is the proton relay mechanism of Pinto et al. [27]. This hypothesis assumes only single ionization, but introduces the idea that the proton picked up by the imidazole ring from the upstream (N-terminus) water concentration can be ionized o again on the downstream side of the channel block, following a spontaneous rotation of the imidazole ring. This mechanism does not require a continuous water column to be present in the conducting state of the channel. The key issues which we address regarding the M2 simulation are then: is there a stable peptide bundle, to what extent does a water column penetrate the pore despite the absence of polar pore-lining groups, and how is the water column aected by single and multiple imidazole protonation, thus throwing light on the gating mechanism. 3. Methodology Recent methodological developments in MD have been aimed at larger systems and longer timescales as well as a more accurate description of the interactions. Ecient Ewald summation of the electrostatic interactions and multiple time-step algorithms are two of the important developments that have enabled us to implement the series of ion channel simulations to be described. As regards the Zhong et al. work [13±17], described in this article, simulations have been implemented within one signi®cant approximation, the replacement [4,29] of the lipid bilayer in which the ion channel is embedded by an octane [30]
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
969
lipid-mimic. The motivation for this approximation is twofold. Our interest lies primarily in aspects of the channel relating to ion transport and thus with the pore. The lipid/channel interface is furthest from the pore, and thus is anticipated to play a secondary role, a role primarily concerned with the thermodynamics of overall channel stability, suggesting that this is an aspect of the problem open to approximation. We believe that the success of the approximation is justi®ed by results. The second motivation is that the lipid bilayer is found to be a relatively slow equilibrating structure, a disadvantage absent from the octane lipid-mimic. Hence, the productivity of the project, which is a computation-limited one, has been substantially enhanced by the octane substitution. The setup procedure and the simulation methods we employ have been described in detail elsewhere [13]. We shall describe the procedure used in the LS3 simulations [14], which can be taken as typical, apart from slight variants, of our approach. We use INSIGHT (BIOSYM Technologies, San Diego, CA) to set up our systems, starting with a bundle of ideal alpha helices arranged as parallel rods in a symmetric manner. The bundle is then minimized and equilibrated with backbone atoms ®xed using an optimized version of the CHARMM program [31], so that the side chain atoms can relax. The bundle is then immersed in a preequilibrated box of octane C8 H18 cut to the same height as the bundle. The octane slab provides the membrane mimetic hydrophobic environment. Layers of water H2 O are added at the two ends of the helices. The composite system is then equilibrated with ®xed helices by carrying out an NVE MD run on the H2 O and C8 H18 to eliminate any possible tension between the dierent components. The whole system is equilibrated for about 200 ps. with the velocities reassigned every 0.1 ps. In order to get optimal performance in both system size and timescale, we use the technique of explicit reversible integrators combined with a multiple time-step scheme [32]. The van der Waals and electrostatic interactions are each divided into for electrostatic intershort and long range parts. The short range cut-o is 7.0 A actions and 2.0r, where r is the usual Lennard±Jones parameter, for van der Waals interactions, respectively. The bonding interactions, including stretching, bending and dihedral angles (both proper and improper), are calculated every 0:3 fs. The short range interactions have to be calculated every 1:5 fs, however they only contain a small portion of the non-bonding interactions. The major part of the non-bonding interaction is the long range part, which is calculated every 3 fs. for the elecPeriodic boundary conditions are introduced with a cut-o of 10 A trostatic interactions and 2.5r for the van der Waals interactions, respectively. Beyond the cut-o, the van der Waals interaction is taken into account by a standard long range correction on an average basis [33,34]. The Ewald method is used to take into account the long range contribution from the electrostatic interaction [33,34]. We ÿ1 , ÿ1 as the cut-o in k-space and the calculation is converged at a 0:3 A use 10 A where a is the Gaussian damping factor. We use a Nose-Hoover chain of length 3 to control the temperature. The frequency of the chain is chosen as 2 psÿ1 . The resulting equations of motion give continuous dynamics that generates a canonical distribution [35,36].
970
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
The well-known TIP3P model [37] is used for both the bulk and pore H2 O. However, since we are using the multiple time-step integrator, we do not introduce constraints on the O±H bond length. The parameters used for C8 H18 are for a fully ¯exible molecule, which consists of methyl and methylene groups [30]. The topology and parameter set for the peptide are based on the CHARMM parameter set (CHARMM 19) [31], which is consistent with the molecular model used for C8 H18 . All hydrogens, except the polar hydrogen and the amide hydrogen, are absorbed into heavy atoms so that the simulation can be carried out more eciently. The simulations are all carried out at 300 K, room temperature, on an unconstrained system. The MD code was developed at the University of Pennsylvania [38], and parallelized using MPI. Computations were carried out on the SP1 and SP2 supercomputers at the T.J. Watson Research Center, at the Maui High Performance Computing Center, at the Cornell Supercomputing Center, and at the Center for Molecular Modelling at the University of Pennsylvania. 4. Results In Fig. 1, we illustrate an axial view, normal to the lipid plane of an equilibrated (4 ns) LS2 channel (Ref. [13]). The peptide chains remain formed into alpha helices (maintaining their setup structure) in this and in all the following ion channel simulations. The alpha helices, however, deviate from their setup orientation (parallel to the pore axis), being now locked into a Crick [39] left-handed coiled-coil structure (see Fig. 2), with a mean tilt angle of 17:5 . In the coiled-coil structure, the hydrophilic OH groups line the pore, as expected from the design philosophy for these LSn channels [18]. The pore contains a water column of about 40 H2 O molecules, compared to 15 molecules in the setup; the strong hydration which has occurred is attributable to the hydrophilic lining of the pore. The physical properties of the pore water have been investigated [13]. The number of hydrogen bonds per water is found to be about 70% of that in TIP3P water [37]. The H2 O self-diusion coecient is about 1/3 of that in bulk water. There are also some anomalies in the dipole relaxation rate relative to bulk H2 O [13]. The pore water properties seem close enough to those of bulk water to support a conventional proton transport mechanism through these synthetic channels [29,40]. Simulation of the LS3 channel (Ref. [14]) has been implemented with both six and seven peptide monomers. The heptamer was found to be highly unstable. In Fig. 3, left panel, we illustrate the equilibrated LS3 hexamer. Although the hexamer also does not appear to form a stable bundle in the simulation, robust hydrogen-bonded peptide dimers are found and persist as a structural motif. The dimer consists of two alpha helices shifted with respect to each other by about 0.75 of a turn, with a twist angle of about 7 . In this way, the leucine hydrophobic side chains of the neighboring helices pack into a ``knob-to-hole'' pattern and the serine hydrophilic side chains form interhelix hydrogen bonds, which plausibly accounts for the stability of the dimers. Note that the arguments based on Eq. (1) suggested that interhelix hydrogen
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
971
Fig. 1. Instantaneous con®guration of LS2 channel after 4 ns, axial view from N-terminus. Atoms drawn with van der Waals radii. Oxygen of serine OH in green, carbonyl oxygens in red, amide nitrogen in blue, amide hydrogen in gray. All other atoms omitted for clarity. Fig. 2. Instantaneous con®guration of LS2 channel after 5 ns, side view with N-terminus at the top. The helices which are drawn in ribbon representation have formed a left-handed coiled-coil structure. Fig. 3. Instantaneous con®guration of LS3 channel, after 1 ns. Left: in absence of electrostatic ®eld channel forms three associated dimers. Right: in presence of electrostatic ®eld, the channel adopts a left-handed coiledcoil structure and the stability of the channel has changed dramatically. Atoms drawn with van der Waals radii, color code is: red (oxygen), black (carbon), blue (nitrogen) and white (hydrogen). Water and octane omitted for clarity.
972
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
bonding could be signi®cant. It is thought that on a more extended timescale the hexamer channel, if formed, may be unstable with respect to dissociation into the dimers. To simulate the gating voltage [14] in this voltage-gated channel [23], we introduced an electric ®eld into the simulation as a smoothed saw-tooth potential. In the presence of the ®eld, the MD simulation shows (Fig. 3, right panel) that the helices adopt a coiled-coil structure [39]. The bundle, which consists of three dimers, has roughly 3-fold symmetry. The tendency towards dissociation along the channel axis is inhibited by the presence of the external ®eld (see Fig. 3). Although we cannot simulate the actual process of gating with the available techniques and resources, the observed change of stability caused by the ®eld is consistent with the idea of a voltage-gated function. In the coiled-coil helical bundle, the axis of each alpha helix is inclined towards the pore by an angle of about 20°. Calculations of the dihedral angles / and w show that the peptides remain predominantly alpha-helical, despite thermal ¯uctuations. Analysis of the structure shows that the helical bundle supports a pore with which is about 0.5 A larger at the N-termian average inner radius of about 5 A, nus. The number of H2 O molecules in the pore is quite stable at 95 10. The physical properties of the pore H2 O are similar to those reported for the LS2 channel. We now turn to the M2 channel (Refs. [15±17]). Since the protonation of the H37 residues is believed to be the key to gating and stability of the M2 channel, we studied neutral [15], singly charged [16,17], doubly (on diagonally opposite H37's) charged [16,17], and quadruply charged [15] M2 channels. Charges were built up incrementally so as to maintain stability in an attempt at mimesis of the slow charge motion involved in a proton migration event. Fig. 4 shows the equilibrated neutral channel after 3 ns. It is seen that the alphahelices form a coiled-coil structure qualitatively similar to that of the LSN channels, and qualitatively similar to that proposed earlier [10]. Despite the absence of a hydrophilically lined pore, there is substantial H2 O penetration into both ends of the channel, although there is no continuous water column through the channel due to the presence of the channel blocking histidine groups, which are highlighted in Fig. 5. Comparison of Figs. 4 and 5 shows that the break in the water column occurs exactly at the location where the histidine H37 groups block the channel. The channel blocking group, and its interruption of the water column, are responsible for the absence of proton transport when the pH is above the ionization threshold of the imidazole ring. A characteristic reminiscent of LS3 is the presence of interhelix hydrogen bonding at the C-terminus, speci®cally between D44 and the guanadino of R45 of adjacent helices. At the N-terminus, it is suggested that the hydration of the S23 and D24 residues, which lie at the octane/water interface, plays a role in the formation of the funnel-like structure at the entrance to the channel. The four solvated D24 residues may even serve to attract extracellular protons into the funnel, whence they can migrate to reach the H37 groups, by a mechanism very similar to that proposed recently for the K channel selector [41,42].
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
973
In a more acidic environment pH 6 5:8, one or more of the imidazole rings ionize [28] by acquisition of a proton. MD simulation of the singly protonated channel yields a picture similar to that of the neutral channel. There is a coiled-coil structure, whose angle of tilt, 28°, is similar to that obtained in recent NMR measurements 33° [43]. But the H37 residues remain close packed so that there is no continuous H2 O channel, although H2 O molecules can penetrate the channel from either direction, as in the neutral M2 channel (Fig. 4). Although without quantum simulation we cannot conclude that the singly protonated channel is a conducting channel, there are features of the simulation which are consistent with the proton relay mechanism of Pinto et al. [27]. First of all, with a liquid H2 O region penetrating far into the channel from the N-terminus, a proton could ®nd its way to the H37 channel block. Secondly, we observe several imidazole ring ¯ips during the simulation (3 ns) [16,17], ¯ips essential to the proton relay mechanism [27]. In the doubly protonated channel, a continuous water column appears (see Fig. 6), a dramatic change in channel con®guration attributable to the electrostatic repulsion between the opposite sides of the channel acting as a channel opening force. In the doubly charged state, the secondary structure and the coiled-coil con®guration of the M2 channel remained stable. We have calculated the inertia tensor of the peptide bundle, which shows that the bundle ¯uctuates around its equilibrium position without any sign of dissociation or collapse. The number of water molecules in the channel pore has been calculated 45 15, and again shows that the whole structure is stable without large ¯uctuations. The opening of the H37 channel block occurs via the rotation of one of the charged H37 imidazole rings towards the Nterminus as in a ¯ap valve. In this `open' structure, which forms a continuous water network, proton transport through the channel could occur in a conventional way [29]. Finally, simulation of the quadruply ionized channel [15] originally postulated by Sansom [10], within a constrained simulation technique, as a possible candidate for the `open' channel state, was carried out. Disappointingly, our unconstrained simulation showed clear signs of instability [15], attributable to the overstrong electrostatic repulsion between the alpha helices in the bundle. 5. Conclusion Our MD studies have revealed a degree of structural universality, but also quite a variety of individual behavior, in the ion channels we have studied. The universal features are the stability of the alpha-helix structure in the transmembrane sections of the multi-helix bundles, and their self-assembly into coiled-coil structures with a hydrated pore. In the case of the synthetic LSN channels, the pore is, as anticipated, lined with hydrophilic serine OH groups, and the water column is large and forms a continuous path through the ion channel. The physical properties of the pore water have been studied in the LS2 case, revealing it to be close to, but not identical with, bulk water
974
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
Fig. 4. Instantaneous con®guration of the M2 channel at 3 ns. Side view, N-terminus at top. The lefthanded supercoil is evident, the interior of the bundle contains water, but channel water column is blocked by the H37 residues. Peptides rendered in ribbon representation, channel waters drawn with van der Waals radii. Octane and non-pore water omitted for clarity. Fig. 5. The same con®guration as in the previous ®gure, showing the H37 channel blocking residues. Fig. 6. Instantaneous con®guration of the M2 channel with a pair of opposite H37 residues protonated. Side view, N-terminus at the top. Now, there is a continuous water column permeating the pore. Channel waters drawn with van der Waals radii. Carbonyl oxygen in red, amide nitrogen and hydrogen in blue and gray, respectively. Octane and non-pore water omitted for clarity.
in its properties. In the LS3 hexamer case, a tendency towards break up into three hydrogen-bonded dimers was revealed, which may indicate instability of the hexamer LS3 channel on a more extended timescale. However, in the presence of the
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
975
gating voltage, the tightly bound coiled-coil structure reasserted itself. It is somewhat speculative to identify the voltage stabilization of the coiled-coil state of the LS3 channel with the gating action, but it is a plausible suggestion whose veri®cation must be left to future work. The M2 channel lacks many pore-lining hydrophilic groups, but nevertheless water penetrates into the pore from both termini, without the neutral channel being able to form a continuous water column due to the presence of the channel blocking histidine residues. When the pH drops below 5.8, one or more of the imidazole rings in the channel blocking region becomes ionized due to protonation. In the singly ionized case, the water column is still not continuous, and the only possibility for the channel to transport protons would be via the proton relay mechanism. However, double ionization of the channel causes the channel blocking region to open due to the electrostatic repulsion between the two charged histidine residues. In this case, there is a continuous water column through the channel, which would allow proton transport via a conventional process of proton diusion through the water matrix. Experimental distinction between the proton relay and multiple ionization mechanisms should be possible, for instance a successful demonstration of pH gating for another small ion capable of transport through the channel would support the multiple ionization hypothesis. Note added in proof Some recent support for the double protonation model of M2 gating exists in the form of measurements by W.F. DeGrado, J. Lear, et al. (private commun.), of the pKa for single ionization of M2. This is found to be about 1.5 points higher than the gating pH, suggesting that double ionization could indeed be the process occuring at the higher proton concentration of the gating pH.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
P.C. Jordan, J. Phys. Chem. 91 (1987) 6582. K.E. Cooper, P.Y. Gates, R.S. Eisenberg, J. Membr. Biol. 106 (1988) 95. J. Aqvist, A. Warshel, Biophys. J. 56 (1989) 171. B. Roux, M. Karplus, Biophys. J. 59 (1991) 961. A. Pullman, Chem. Rev. 91 (1991) 793. B. Roux, B. Prodhom, M. Karplus, Biophys. J. 68 (1995) 876. T.B. Woolf, B. Roux, Proc. Natl. Acad. Sci. USA 91 (1994) 11631. M.S.P. Sansom, H.S.S.R. Sankararamakrishnan, I.D. Kerr, J. Breed, Biophys. J. 68 (1995) 1295. M.S.P. Sansom, H.S.S.R. Sankararamakrishnan, I.D. Kerr, J. Breed, Biophys. J. 70 (1996) 1643. M.S.P. Sansom, I.D. Kerr, G.R. Smith, H.S. Son, Virology 233 (1997) 163. D.P. Tielemann, M.S.P. Sansom, H.J.C. Berendsen, Biophys. J. 76 (1999) 40. D.P. Tielemann, H.J.C. Berendsen, M.S.P. Sansom, Biophys. J. 76 (1999) 1757. Q.F. Zhong, Q. Jiang, P.B. Moore, D.M. Newns, M.L. Klein, Biophys. J. 74 (1998) 3. Q.F. Zhong, P.B. Moore, D.M. Newns, M.L. Klein, FEBS Lett. 427 (1998) 267. Q.F. Zhong, T. Husslein, P.B. Moore, D.M. Newns, P.C. Pattnaik, M.L. Klein, FEBS Lett. 434 (1998) 265.
976 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976 Q.F. Zhong, D.M. Newns, P.C. Pattnaik, J.D. Lear, M.L. Klein, FEBS Lett. 473 (2000) 195±198. T. Husslein, P.B. Moore, Q. Zhong, D.M. Newns, P.C. Pattnaik, M.L. Klein, preprint. J.D. Lear, Z.R. Wasserman, W.F. DeGrado, Science 240 (1988) 1177. G.R. Dieckmann, J.D. Lear, Q. Zhong, M.L. Klein, W.F. DeGrado, K.A. Sharp, Biophys. J. 76 (1999) 618. A.A. Chung, J.D. Lear, W.F. DeGrado, Biochemistry 31 (1992) 6608. K.S. Akerfeldt, J.D. Lear, Z.R. Wasserman, L.A. Chung, W.F. DeGrado, Acc. Chem. Res. 26 (1993) 191. P.K. Kienker, W.F. DeGrado, J.D. Lear, Proc. Natl. Acad. Sci. USA 91 (1994) 4859. P. Kienker, J.D. Lear, Biophys. J. 68 (1995) 1347. R.J. Sugrue, A.J. Hay, Virology 180 (1991) 617. K.C. Du, R.H. Ashley, Virology 190 (1992) 485. C. Wang, K. Takeuchi, H.P. Pinto, R.A. Lamb, J. Virol. 67 (1993) 5585. L.H. Pinto, G.R. Dieckmann, C.S. Gandhi, C.G. Papworth, J. Braman, M.A. Shaughnessy, J.D. Lear, R.A. Lamb, W.F. DeGrado, Proc. Natl. Acad. Sci. USA 94 (1997) 11301. C. Wang, R.A. Lamb, L.H. Pinto, Biophys. J. 69 (1995) 1363. D.E. Sagnella, K. Laasonen, M.L. Klein, Biophys. J. 71 (1996) 1172. J.I. Siepmann, S. Karaborni, B. Smit, Nature 365 (1993) 330. B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan, M. Karplus, J. Comp. Chem. 4 (1983) 187. M. Tuckerman, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. G. Ciccotti, D. Frenkel, I.R. McDonald, Simulation of Liquids and Solids: Molecular Dynamics and Monte Carlo Methods in Statistical Mechanics, Elsevier, Amsterdam, 1987. G.J. Martyna, M.L. Klein, M. Tuckerman, J. Chem. Phys. 97 (1992) 2635. G.J. Martyna, M. Tuckerman, D.J. Tobias, M.L. Klein, Mol. Phys. 87 (1996) 1117. W.L. Jorgensen, J. Chandresekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. P.B. Moore, M.L. Klein, Technical Report, University of Pennsylvania, 1997. F.H.C. Crick, Acta Crystallogr. 6 (1953) 689. R. Pomes, B. Roux, Biophys. J. 71 (1996) 19. D.A. Doyle, J.M. Cabral, R.A. Pfuetzner, A. Kuo, J.M. Gulbis, S.L. Cohen, B.T. Chait, R. MacKinnon, Science 280 (1998) 69. A. Kreusch, P.J. Pfanger, C.F. Stevens, S. Choe, Nature 392 (1998) 945. F.A. Kovacs, T.A. Cross, Biophys. J. 73 (1997) 2511.
www.elsevier.com/locate/parco
Molecular dynamics study of structure and gating of low molecular weight ion channels Dennis M. Newns a,*, Qingfeng Zhong b, Preston B. Moore b, T. Husslein b, Pratap Pattnaik a, Michael L. Klein b a
b
IBM, Thomas J. Watson Research Center, Yorktown Heights, NY 10598, USA Center for Molecular Modeling and Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323, USA Received 12 July 1999; accepted 26 October 1999
Abstract We implement molecular dynamics (MD) simulations on low molecular weight alpha helixbased functional synthetic and native ion channels. The synthetic channels are the LS2 proton channel and the LS3 voltage-gated channel. The simulation manifests key features of the channels such as the coiled-coil structure of the alpha-helix bundle and the continuous aqueous pore. By implementing simulations with and without an applied voltage, we develop a hypothesis as to the voltage-gating mechanism. The native channel is the M2 proton channel in the in¯uenza A virus, which plays an essential role in the infection process. This channel is pH gated via protonation of one or more imidazole rings in the H37 residues. Simulation of the neutral channel reveals a coiled-coil structure whose pore is penetrated by water, but not threaded by a water column. By means of simulations with dierent numbers of charged H37 residues, we demonstrate a possible gating action via opening up the channel to a continuous water column, and also provide support for the alternative proton relay gating mechanism. Ó 2000 Published by Elsevier Science B.V. All rights reserved. Keywords: Ion channel; M2; Proton channel; Molecular dynamics
1. Introduction Molecular dynamics (MD) simulation has reached a stage of development, where it can be realistically applied in environments as complex as protein/water and *
Corresponding author. E-mail address: [email protected] (D.M. Newns).
0167-8191/00/$ - see front matter Ó 2000 Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 1 9 1 ( 0 0 ) 0 0 0 2 1 - 1
966
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
protein/lipid systems. Current computer power places limits on the size of systems and/or the length of time for which systems can be studied. The protein folding problem, involving real times of order ms, poses a major challenge in terms of number of time steps required. However, there are problems of major signi®cance in biochemistry which have time constants of order ns and can be readily addressed with presently available computer capabilities, especially exploiting parallelism. Such problems can also serve as a test of the reliability of MD simulations of large biomolecules. Transmembrane proteins constitute 30% of a typical genome, and among these ion channels are recognized as having extensive biological importance; the central role, e.g., of the K channel in nerve conduction is well known, but ion channels are increasingly recognized as playing key roles in many other areas. The ion channel is a transmembrane protein whose role is, with some degree of selectivity, to transport an ion across the impermeable lipid membrane when subject to a stimulus, which may be, e.g., a voltage across the membrane, or a change in pH. The action of this stimulus in opening the channel to transport is termed `gating'. A typical ion channel structure is a peptide bundle (or a continuous protein) which in its transmembrane section is organized into a bundle of alpha helices. The center of the bundle constitutes the pore for ion transport; to what extent a water column pervades the pore is an issue in the ®eld which is highly relevant to the nature of ion channel action. The peptide sections external to the membrane, in the aqueous medium, may not be coiled. The earliest ion channel studies by MD [1±7] by Roux and others involved a slightly atypical channel, gramicidin, which is a beta barrel rather than an alphahelix bundle. Sansom's [8] group studied seven-helix bundle ion channels, and the behavior of water within various models of ion channels [9]. Sansom's work on the in¯uenza A virus M2 channel [10] provided an important basis of comparison for part of the present work. This group has also carried out an extensive study of the alamethicin channel [11,12]. In the work to be described below, Zhong et al. and Husslein and coworkers [13±17] have studied both functional synthetic channels and the M2 channel, with special emphasis on understanding of the gating mechanism. 2. Objectives The MD technique is intrinsically capable of locating the stable state of the macromolecular system provided (a) sucient simulation time is allowed and (b) no deep local free energy minima are encountered. We have found that in several nanoseconds ion channel conformations such hydrated, coiled-coil bundles can be stably established. However, a result such as bundle dissociation, which requires a longer timescale than nanoseconds to occur fully, can only be established in a more speculative manner within the framework of nanosecond simulations. In this paper, we shall describe MD simulations of two very dierent types of channels, selected for their relatively low molecular weight and the availability of
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
967
extensive experimental data on their properties, including information on their gating behavior. The ®rst systems are a pair of functional synthetic ion channels invented by Lear, Wasserman and DeGrado [18], aided by extensive modeling by Dieckmann et al. [19]. The synthetic channels, both built from 21 amino acid peptide sequences, are termed LS2 and LS3. The LS2 peptide consists of the sequence Ac± LSLLLSL3 ±CONH2 , and LS3 has the sequence Ac± LSSLLSL3 ±CONH2 (where L leucine and S serine). In the presence of a lipid layer, the LSN peptides are believed [18] to self-assemble into a transmembrane bundle of alpha helices, which experimentally [18±23] acts as a functional ion channel. The architecture [18] of the sequence is designed so that, viewed down its axis, each helix has an angular sector h which is hydrophilic, consisting of serine residues, and which is expected to orient towards the pore, while the remaining, larger, angular sector 2p ÿ h is hydrophopic, consisting of leucine residues. Thus, when the hydrophilic sectors face the pore, the bundle is expected to have a hydrophobic coating, enabling it to embed into the lipid bilayer, and a hydrophilic pore, capable of maintaining a transmembrane water column. In a naive model of the channel bundle as a symmetric ring of n touching cylinders with hydrophilic sector inwards, then n
2p 2 L S ; pÿh LÿS
1
where L and S are the number of leucines and serines, respectively. This is only a semiquantitative model, but it predicts that the larger the fraction of serine, the larger is the number n of peptides in the channel. For the LS2 peptide, Eq. (1) estimates n 4:7, and LS2 is believed to form [18] a tetrameric channel. Consistent with the small size of the pore in this case, LS2 is proton selective. For serine-rich LS3, Eq. (1) predicts the large value n 14; but in fact this channel is believed to form a hexameric or heptameric bundle. Experimentally it can transport larger ions such as Clÿ , and is voltage-gated [23]. The excessive hydrophilic sector h in LS3 relative to pore perimeter raises the possibility of signi®cant interhelix hydrogen bonding due to the overlap of OH groups from dierent chains. Molecular dynamics can be used to investigate a number of issues in the LSN channel area. Key issues are, whether the hypothesized bundles are stable, whether their architecture indeed includes the expected feature of pore-lining OH groups, whether there is a water column and how large is it, whether there is interhelix hydrogen bonding in LS3, and how does the voltage-gating of the LS3 channel work. Our second set of studies is of a native channel, the M2 channel of the in¯uenza A virus [24]. At two stages in the infection process of a cell by the virus, the M2 channel plays the role of a proton pump and is crucial to successful infection. In fact, there is an in¯uenza drug, amantadine, which is believed to operate by blocking the M2 channel [25,26]. The gating of the channel is caused by lowering the pH below, approximately pH 5.8. A key simpli®cation in the study of the M2 channel is the discovery by Du and Ashley [25] of a truncated 25 amino acid M2 sequence SSDPLVVAASIIGILHLILWILDRL, consisting of a transmembrane section with small extensions, which nevertheless gates like a complete M2 channel. The
968
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
functional truncated Du±Ashley sequence is ideal for MD simulation. Like the LS2 sequence, the Du±Ashley M2 peptide is believed to self-assemble into a transmembrane alpha-helix bundle in the presence of the lipid bilayer, the assembly being tetrameric [27]. In many ways, however, the M2 contrasts with the synthetic LSN type of system. First, in the M2, in contrast with LSN, there are few hydrophilic groups lining the pore. Secondly, the M2 is believed to contain a well-de®ned gating structure consisting of four histidine residues, one from each peptide in the bundle. The four histidine residues form a channel block towards the C-terminus of the peptides. There are two main hypotheses as to how the pH-controlled gating of the channel may occur. Both hypotheses are based on the gating pH 5.8 lying close to the pKa for protonation of the histidine imidazole ring [28]. In the multiple ionization hypothesis, two or possibly more histidines are protonated below pH 5.8, and the resulting inter-peptide chain electrostatic repulsion (which is largely unscreened due to the absence of a water column at this point) forces apart the bundle at the pore blocking point, leading to the formation of a water column and the consequent opening of the channel to proton transport [29]. The initial version of this multiple ionization hypothesis, due to Sansom [10], assumed ionization of all four histidine groups. However, according to the results of our simulations which are unconstrained [15], the quadruply ionized channel is unstable, while we ®nd that double ionization is sucient to open the channel, without endangering its stability. A second hypothesis for M2 gating is the proton relay mechanism of Pinto et al. [27]. This hypothesis assumes only single ionization, but introduces the idea that the proton picked up by the imidazole ring from the upstream (N-terminus) water concentration can be ionized o again on the downstream side of the channel block, following a spontaneous rotation of the imidazole ring. This mechanism does not require a continuous water column to be present in the conducting state of the channel. The key issues which we address regarding the M2 simulation are then: is there a stable peptide bundle, to what extent does a water column penetrate the pore despite the absence of polar pore-lining groups, and how is the water column aected by single and multiple imidazole protonation, thus throwing light on the gating mechanism. 3. Methodology Recent methodological developments in MD have been aimed at larger systems and longer timescales as well as a more accurate description of the interactions. Ecient Ewald summation of the electrostatic interactions and multiple time-step algorithms are two of the important developments that have enabled us to implement the series of ion channel simulations to be described. As regards the Zhong et al. work [13±17], described in this article, simulations have been implemented within one signi®cant approximation, the replacement [4,29] of the lipid bilayer in which the ion channel is embedded by an octane [30]
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
969
lipid-mimic. The motivation for this approximation is twofold. Our interest lies primarily in aspects of the channel relating to ion transport and thus with the pore. The lipid/channel interface is furthest from the pore, and thus is anticipated to play a secondary role, a role primarily concerned with the thermodynamics of overall channel stability, suggesting that this is an aspect of the problem open to approximation. We believe that the success of the approximation is justi®ed by results. The second motivation is that the lipid bilayer is found to be a relatively slow equilibrating structure, a disadvantage absent from the octane lipid-mimic. Hence, the productivity of the project, which is a computation-limited one, has been substantially enhanced by the octane substitution. The setup procedure and the simulation methods we employ have been described in detail elsewhere [13]. We shall describe the procedure used in the LS3 simulations [14], which can be taken as typical, apart from slight variants, of our approach. We use INSIGHT (BIOSYM Technologies, San Diego, CA) to set up our systems, starting with a bundle of ideal alpha helices arranged as parallel rods in a symmetric manner. The bundle is then minimized and equilibrated with backbone atoms ®xed using an optimized version of the CHARMM program [31], so that the side chain atoms can relax. The bundle is then immersed in a preequilibrated box of octane C8 H18 cut to the same height as the bundle. The octane slab provides the membrane mimetic hydrophobic environment. Layers of water H2 O are added at the two ends of the helices. The composite system is then equilibrated with ®xed helices by carrying out an NVE MD run on the H2 O and C8 H18 to eliminate any possible tension between the dierent components. The whole system is equilibrated for about 200 ps. with the velocities reassigned every 0.1 ps. In order to get optimal performance in both system size and timescale, we use the technique of explicit reversible integrators combined with a multiple time-step scheme [32]. The van der Waals and electrostatic interactions are each divided into for electrostatic intershort and long range parts. The short range cut-o is 7.0 A actions and 2.0r, where r is the usual Lennard±Jones parameter, for van der Waals interactions, respectively. The bonding interactions, including stretching, bending and dihedral angles (both proper and improper), are calculated every 0:3 fs. The short range interactions have to be calculated every 1:5 fs, however they only contain a small portion of the non-bonding interactions. The major part of the non-bonding interaction is the long range part, which is calculated every 3 fs. for the elecPeriodic boundary conditions are introduced with a cut-o of 10 A trostatic interactions and 2.5r for the van der Waals interactions, respectively. Beyond the cut-o, the van der Waals interaction is taken into account by a standard long range correction on an average basis [33,34]. The Ewald method is used to take into account the long range contribution from the electrostatic interaction [33,34]. We ÿ1 , ÿ1 as the cut-o in k-space and the calculation is converged at a 0:3 A use 10 A where a is the Gaussian damping factor. We use a Nose-Hoover chain of length 3 to control the temperature. The frequency of the chain is chosen as 2 psÿ1 . The resulting equations of motion give continuous dynamics that generates a canonical distribution [35,36].
970
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
The well-known TIP3P model [37] is used for both the bulk and pore H2 O. However, since we are using the multiple time-step integrator, we do not introduce constraints on the O±H bond length. The parameters used for C8 H18 are for a fully ¯exible molecule, which consists of methyl and methylene groups [30]. The topology and parameter set for the peptide are based on the CHARMM parameter set (CHARMM 19) [31], which is consistent with the molecular model used for C8 H18 . All hydrogens, except the polar hydrogen and the amide hydrogen, are absorbed into heavy atoms so that the simulation can be carried out more eciently. The simulations are all carried out at 300 K, room temperature, on an unconstrained system. The MD code was developed at the University of Pennsylvania [38], and parallelized using MPI. Computations were carried out on the SP1 and SP2 supercomputers at the T.J. Watson Research Center, at the Maui High Performance Computing Center, at the Cornell Supercomputing Center, and at the Center for Molecular Modelling at the University of Pennsylvania. 4. Results In Fig. 1, we illustrate an axial view, normal to the lipid plane of an equilibrated (4 ns) LS2 channel (Ref. [13]). The peptide chains remain formed into alpha helices (maintaining their setup structure) in this and in all the following ion channel simulations. The alpha helices, however, deviate from their setup orientation (parallel to the pore axis), being now locked into a Crick [39] left-handed coiled-coil structure (see Fig. 2), with a mean tilt angle of 17:5 . In the coiled-coil structure, the hydrophilic OH groups line the pore, as expected from the design philosophy for these LSn channels [18]. The pore contains a water column of about 40 H2 O molecules, compared to 15 molecules in the setup; the strong hydration which has occurred is attributable to the hydrophilic lining of the pore. The physical properties of the pore water have been investigated [13]. The number of hydrogen bonds per water is found to be about 70% of that in TIP3P water [37]. The H2 O self-diusion coecient is about 1/3 of that in bulk water. There are also some anomalies in the dipole relaxation rate relative to bulk H2 O [13]. The pore water properties seem close enough to those of bulk water to support a conventional proton transport mechanism through these synthetic channels [29,40]. Simulation of the LS3 channel (Ref. [14]) has been implemented with both six and seven peptide monomers. The heptamer was found to be highly unstable. In Fig. 3, left panel, we illustrate the equilibrated LS3 hexamer. Although the hexamer also does not appear to form a stable bundle in the simulation, robust hydrogen-bonded peptide dimers are found and persist as a structural motif. The dimer consists of two alpha helices shifted with respect to each other by about 0.75 of a turn, with a twist angle of about 7 . In this way, the leucine hydrophobic side chains of the neighboring helices pack into a ``knob-to-hole'' pattern and the serine hydrophilic side chains form interhelix hydrogen bonds, which plausibly accounts for the stability of the dimers. Note that the arguments based on Eq. (1) suggested that interhelix hydrogen
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
971
Fig. 1. Instantaneous con®guration of LS2 channel after 4 ns, axial view from N-terminus. Atoms drawn with van der Waals radii. Oxygen of serine OH in green, carbonyl oxygens in red, amide nitrogen in blue, amide hydrogen in gray. All other atoms omitted for clarity. Fig. 2. Instantaneous con®guration of LS2 channel after 5 ns, side view with N-terminus at the top. The helices which are drawn in ribbon representation have formed a left-handed coiled-coil structure. Fig. 3. Instantaneous con®guration of LS3 channel, after 1 ns. Left: in absence of electrostatic ®eld channel forms three associated dimers. Right: in presence of electrostatic ®eld, the channel adopts a left-handed coiledcoil structure and the stability of the channel has changed dramatically. Atoms drawn with van der Waals radii, color code is: red (oxygen), black (carbon), blue (nitrogen) and white (hydrogen). Water and octane omitted for clarity.
972
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
bonding could be signi®cant. It is thought that on a more extended timescale the hexamer channel, if formed, may be unstable with respect to dissociation into the dimers. To simulate the gating voltage [14] in this voltage-gated channel [23], we introduced an electric ®eld into the simulation as a smoothed saw-tooth potential. In the presence of the ®eld, the MD simulation shows (Fig. 3, right panel) that the helices adopt a coiled-coil structure [39]. The bundle, which consists of three dimers, has roughly 3-fold symmetry. The tendency towards dissociation along the channel axis is inhibited by the presence of the external ®eld (see Fig. 3). Although we cannot simulate the actual process of gating with the available techniques and resources, the observed change of stability caused by the ®eld is consistent with the idea of a voltage-gated function. In the coiled-coil helical bundle, the axis of each alpha helix is inclined towards the pore by an angle of about 20°. Calculations of the dihedral angles / and w show that the peptides remain predominantly alpha-helical, despite thermal ¯uctuations. Analysis of the structure shows that the helical bundle supports a pore with which is about 0.5 A larger at the N-termian average inner radius of about 5 A, nus. The number of H2 O molecules in the pore is quite stable at 95 10. The physical properties of the pore H2 O are similar to those reported for the LS2 channel. We now turn to the M2 channel (Refs. [15±17]). Since the protonation of the H37 residues is believed to be the key to gating and stability of the M2 channel, we studied neutral [15], singly charged [16,17], doubly (on diagonally opposite H37's) charged [16,17], and quadruply charged [15] M2 channels. Charges were built up incrementally so as to maintain stability in an attempt at mimesis of the slow charge motion involved in a proton migration event. Fig. 4 shows the equilibrated neutral channel after 3 ns. It is seen that the alphahelices form a coiled-coil structure qualitatively similar to that of the LSN channels, and qualitatively similar to that proposed earlier [10]. Despite the absence of a hydrophilically lined pore, there is substantial H2 O penetration into both ends of the channel, although there is no continuous water column through the channel due to the presence of the channel blocking histidine groups, which are highlighted in Fig. 5. Comparison of Figs. 4 and 5 shows that the break in the water column occurs exactly at the location where the histidine H37 groups block the channel. The channel blocking group, and its interruption of the water column, are responsible for the absence of proton transport when the pH is above the ionization threshold of the imidazole ring. A characteristic reminiscent of LS3 is the presence of interhelix hydrogen bonding at the C-terminus, speci®cally between D44 and the guanadino of R45 of adjacent helices. At the N-terminus, it is suggested that the hydration of the S23 and D24 residues, which lie at the octane/water interface, plays a role in the formation of the funnel-like structure at the entrance to the channel. The four solvated D24 residues may even serve to attract extracellular protons into the funnel, whence they can migrate to reach the H37 groups, by a mechanism very similar to that proposed recently for the K channel selector [41,42].
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
973
In a more acidic environment pH 6 5:8, one or more of the imidazole rings ionize [28] by acquisition of a proton. MD simulation of the singly protonated channel yields a picture similar to that of the neutral channel. There is a coiled-coil structure, whose angle of tilt, 28°, is similar to that obtained in recent NMR measurements 33° [43]. But the H37 residues remain close packed so that there is no continuous H2 O channel, although H2 O molecules can penetrate the channel from either direction, as in the neutral M2 channel (Fig. 4). Although without quantum simulation we cannot conclude that the singly protonated channel is a conducting channel, there are features of the simulation which are consistent with the proton relay mechanism of Pinto et al. [27]. First of all, with a liquid H2 O region penetrating far into the channel from the N-terminus, a proton could ®nd its way to the H37 channel block. Secondly, we observe several imidazole ring ¯ips during the simulation (3 ns) [16,17], ¯ips essential to the proton relay mechanism [27]. In the doubly protonated channel, a continuous water column appears (see Fig. 6), a dramatic change in channel con®guration attributable to the electrostatic repulsion between the opposite sides of the channel acting as a channel opening force. In the doubly charged state, the secondary structure and the coiled-coil con®guration of the M2 channel remained stable. We have calculated the inertia tensor of the peptide bundle, which shows that the bundle ¯uctuates around its equilibrium position without any sign of dissociation or collapse. The number of water molecules in the channel pore has been calculated 45 15, and again shows that the whole structure is stable without large ¯uctuations. The opening of the H37 channel block occurs via the rotation of one of the charged H37 imidazole rings towards the Nterminus as in a ¯ap valve. In this `open' structure, which forms a continuous water network, proton transport through the channel could occur in a conventional way [29]. Finally, simulation of the quadruply ionized channel [15] originally postulated by Sansom [10], within a constrained simulation technique, as a possible candidate for the `open' channel state, was carried out. Disappointingly, our unconstrained simulation showed clear signs of instability [15], attributable to the overstrong electrostatic repulsion between the alpha helices in the bundle. 5. Conclusion Our MD studies have revealed a degree of structural universality, but also quite a variety of individual behavior, in the ion channels we have studied. The universal features are the stability of the alpha-helix structure in the transmembrane sections of the multi-helix bundles, and their self-assembly into coiled-coil structures with a hydrated pore. In the case of the synthetic LSN channels, the pore is, as anticipated, lined with hydrophilic serine OH groups, and the water column is large and forms a continuous path through the ion channel. The physical properties of the pore water have been studied in the LS2 case, revealing it to be close to, but not identical with, bulk water
974
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
Fig. 4. Instantaneous con®guration of the M2 channel at 3 ns. Side view, N-terminus at top. The lefthanded supercoil is evident, the interior of the bundle contains water, but channel water column is blocked by the H37 residues. Peptides rendered in ribbon representation, channel waters drawn with van der Waals radii. Octane and non-pore water omitted for clarity. Fig. 5. The same con®guration as in the previous ®gure, showing the H37 channel blocking residues. Fig. 6. Instantaneous con®guration of the M2 channel with a pair of opposite H37 residues protonated. Side view, N-terminus at the top. Now, there is a continuous water column permeating the pore. Channel waters drawn with van der Waals radii. Carbonyl oxygen in red, amide nitrogen and hydrogen in blue and gray, respectively. Octane and non-pore water omitted for clarity.
in its properties. In the LS3 hexamer case, a tendency towards break up into three hydrogen-bonded dimers was revealed, which may indicate instability of the hexamer LS3 channel on a more extended timescale. However, in the presence of the
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976
975
gating voltage, the tightly bound coiled-coil structure reasserted itself. It is somewhat speculative to identify the voltage stabilization of the coiled-coil state of the LS3 channel with the gating action, but it is a plausible suggestion whose veri®cation must be left to future work. The M2 channel lacks many pore-lining hydrophilic groups, but nevertheless water penetrates into the pore from both termini, without the neutral channel being able to form a continuous water column due to the presence of the channel blocking histidine residues. When the pH drops below 5.8, one or more of the imidazole rings in the channel blocking region becomes ionized due to protonation. In the singly ionized case, the water column is still not continuous, and the only possibility for the channel to transport protons would be via the proton relay mechanism. However, double ionization of the channel causes the channel blocking region to open due to the electrostatic repulsion between the two charged histidine residues. In this case, there is a continuous water column through the channel, which would allow proton transport via a conventional process of proton diusion through the water matrix. Experimental distinction between the proton relay and multiple ionization mechanisms should be possible, for instance a successful demonstration of pH gating for another small ion capable of transport through the channel would support the multiple ionization hypothesis. Note added in proof Some recent support for the double protonation model of M2 gating exists in the form of measurements by W.F. DeGrado, J. Lear, et al. (private commun.), of the pKa for single ionization of M2. This is found to be about 1.5 points higher than the gating pH, suggesting that double ionization could indeed be the process occuring at the higher proton concentration of the gating pH.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
P.C. Jordan, J. Phys. Chem. 91 (1987) 6582. K.E. Cooper, P.Y. Gates, R.S. Eisenberg, J. Membr. Biol. 106 (1988) 95. J. Aqvist, A. Warshel, Biophys. J. 56 (1989) 171. B. Roux, M. Karplus, Biophys. J. 59 (1991) 961. A. Pullman, Chem. Rev. 91 (1991) 793. B. Roux, B. Prodhom, M. Karplus, Biophys. J. 68 (1995) 876. T.B. Woolf, B. Roux, Proc. Natl. Acad. Sci. USA 91 (1994) 11631. M.S.P. Sansom, H.S.S.R. Sankararamakrishnan, I.D. Kerr, J. Breed, Biophys. J. 68 (1995) 1295. M.S.P. Sansom, H.S.S.R. Sankararamakrishnan, I.D. Kerr, J. Breed, Biophys. J. 70 (1996) 1643. M.S.P. Sansom, I.D. Kerr, G.R. Smith, H.S. Son, Virology 233 (1997) 163. D.P. Tielemann, M.S.P. Sansom, H.J.C. Berendsen, Biophys. J. 76 (1999) 40. D.P. Tielemann, H.J.C. Berendsen, M.S.P. Sansom, Biophys. J. 76 (1999) 1757. Q.F. Zhong, Q. Jiang, P.B. Moore, D.M. Newns, M.L. Klein, Biophys. J. 74 (1998) 3. Q.F. Zhong, P.B. Moore, D.M. Newns, M.L. Klein, FEBS Lett. 427 (1998) 267. Q.F. Zhong, T. Husslein, P.B. Moore, D.M. Newns, P.C. Pattnaik, M.L. Klein, FEBS Lett. 434 (1998) 265.
976 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
D.M. Newns et al. / Parallel Computing 26 (2000) 965±976 Q.F. Zhong, D.M. Newns, P.C. Pattnaik, J.D. Lear, M.L. Klein, FEBS Lett. 473 (2000) 195±198. T. Husslein, P.B. Moore, Q. Zhong, D.M. Newns, P.C. Pattnaik, M.L. Klein, preprint. J.D. Lear, Z.R. Wasserman, W.F. DeGrado, Science 240 (1988) 1177. G.R. Dieckmann, J.D. Lear, Q. Zhong, M.L. Klein, W.F. DeGrado, K.A. Sharp, Biophys. J. 76 (1999) 618. A.A. Chung, J.D. Lear, W.F. DeGrado, Biochemistry 31 (1992) 6608. K.S. Akerfeldt, J.D. Lear, Z.R. Wasserman, L.A. Chung, W.F. DeGrado, Acc. Chem. Res. 26 (1993) 191. P.K. Kienker, W.F. DeGrado, J.D. Lear, Proc. Natl. Acad. Sci. USA 91 (1994) 4859. P. Kienker, J.D. Lear, Biophys. J. 68 (1995) 1347. R.J. Sugrue, A.J. Hay, Virology 180 (1991) 617. K.C. Du, R.H. Ashley, Virology 190 (1992) 485. C. Wang, K. Takeuchi, H.P. Pinto, R.A. Lamb, J. Virol. 67 (1993) 5585. L.H. Pinto, G.R. Dieckmann, C.S. Gandhi, C.G. Papworth, J. Braman, M.A. Shaughnessy, J.D. Lear, R.A. Lamb, W.F. DeGrado, Proc. Natl. Acad. Sci. USA 94 (1997) 11301. C. Wang, R.A. Lamb, L.H. Pinto, Biophys. J. 69 (1995) 1363. D.E. Sagnella, K. Laasonen, M.L. Klein, Biophys. J. 71 (1996) 1172. J.I. Siepmann, S. Karaborni, B. Smit, Nature 365 (1993) 330. B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan, M. Karplus, J. Comp. Chem. 4 (1983) 187. M. Tuckerman, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. G. Ciccotti, D. Frenkel, I.R. McDonald, Simulation of Liquids and Solids: Molecular Dynamics and Monte Carlo Methods in Statistical Mechanics, Elsevier, Amsterdam, 1987. G.J. Martyna, M.L. Klein, M. Tuckerman, J. Chem. Phys. 97 (1992) 2635. G.J. Martyna, M. Tuckerman, D.J. Tobias, M.L. Klein, Mol. Phys. 87 (1996) 1117. W.L. Jorgensen, J. Chandresekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. P.B. Moore, M.L. Klein, Technical Report, University of Pennsylvania, 1997. F.H.C. Crick, Acta Crystallogr. 6 (1953) 689. R. Pomes, B. Roux, Biophys. J. 71 (1996) 19. D.A. Doyle, J.M. Cabral, R.A. Pfuetzner, A. Kuo, J.M. Gulbis, S.L. Cohen, B.T. Chait, R. MacKinnon, Science 280 (1998) 69. A. Kreusch, P.J. Pfanger, C.F. Stevens, S. Choe, Nature 392 (1998) 945. F.A. Kovacs, T.A. Cross, Biophys. J. 73 (1997) 2511.