INTRODUCTION

TOP ABSTRACT INTRODUCTION MODELS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS AND REMARKS REFERENCES

Although cell lysis by toxic peptides is a well-known phenomenon and has been investigated by many experiments, the molecular mechanism is not well understood. A particular popular toxin is the bee venom peptide melittin, which is known to cause hemolysis (DeGrado et al., 1982; Tosteson et al., 1985; Katsu et al., 1988) and to induce leakage of fluorescent dyes from lipid vesicles (Schwarz et al., 1992; Benachir and Lafleur, 1995; Rex, 1996; Ladokhin et al., 1997; Matsuzaki et al., 1997b; Rex and Schwarz, 1998).

Melittin's secondary structure is well established to be highly -helical in its crystalline state (Terwilliger and Eisenberg, 1982; Dempsey, 1990) and may form some type of tetrameric aggregate (Terwilliger et al., 1982; Dempsey, 1990). The interaction with membranes, however, changes the properties. The currently accepted view is that under certain conditions melittin molecules insert into a lipid bilayer and form multiple aggregated forms that are affected by temperature, pH value, ionic strength, lipid composition, and lipid-to-peptide ratio. A central problem which has been addressed in many experiments is the type of aggregate and its lytic mechanism. A reasonable view at the present time is the formation of cylindrical pores built by transbilayer helices (Vogel and Jähnig, 1986; Schwarz et al., 1992; Rex, 1996; Ladokhin et al., 1997; Matsuzaki et al., 1997b). In most cases the existence of pores has been concluded on the basis of the efflux of fluorescent dye molecules from large unilamellar vesicles.

However, no evidence for spontaneous aggregation has been reported in other experimental studies (Hermetter and Lakowicz, 1986; Schwarz and Beschiaschvili, 1989; John and Jähnig, 1991). It had been proposed (Talbot et al., 1987; John and Jähnig, 1991) that aggregation takes place at high salt concentration and at peptide/lipid molar ratio above 1:200. These conclusions are also in accord with other experimental results obtained by NMR (Stanislawski and Rüterjans, 1987) and ERS (Altenbach and Hubbell, 1988).

It should be noted that other intriguing mechanisms of melittin-induced membrane permeability can be envisaged. For example, leakage may be induced by melittin due to a cooperative perturbation of the bilayer's permeability (Benachir and Lafleur, 1995).

The sizes of the melittin pores, characterized by the inner pore diameter, have been reported to be in the range 1.0-6.0 nm (Rex, 1996), 1.3-2.4 nm (Matsuzaki et al., 1997b), and 2.5-3.0 nm (Ladokhin et al., 1997). The corresponding numbers n of melittin monomers forming the pores were estimated using the formula n  (d_P/d_M + 1), where d_P and d_M  1.2 nm are the pore diameter and the average diameter of the helix of melittin (Terwilliger and Eisenberg, 1982), respectively. Using the above values for d_P, one finds that the pore should consist of 6-19 monomers. Based on other considerations, however, even smaller pore sizes consisting of four monomers have been postulated (Tosteson and Tosteson, 1981; Vogel and Jähnig, 1986), which has not be confirmed by efflux measurements.

The aim of the present work is to investigate the stability of a hypothetical melittin pore consisting of four monomers using molecular dynamics simulations. Similarly to previous experiments (Schwarz et al., 1992; Benachir and Lafleur, 1995; Rex, 1996; Ladokhin et al., 1997), we consider the pore to be embedded in a POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine) bilayer membrane. The pH value is 7. It should be noted that the explicit simulation of lipids is important for the present study. So far, only a few simulations of helix aggregates in fully hydrated lipid membranes have been reported in the time range >1 ns (Tieleman and Berendsen, 1998; Zhong et al., 1998; Tieleman et al., 1999). Because performing a simulation of the complete system consisting of water, lipids, and proteins is very time-consuming, in many studies the membrane has been replaced by an effective medium and the pore or membrane protein structure has been preserved by "restrained methods" (Sansom et al., 1995). This latter procedure, however, is not suitable for the present model system, as will become apparent from the results presented below.

	MODELS AND METHODS

TOP ABSTRACT INTRODUCTION MODELS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS AND REMARKS REFERENCES

Initial membrane configuration

The bilayer lipid membrane consists of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine) lipids. This membrane has a low gel-to-fluid transition temperature of T_m = 5°C. The AMBER '91 parameter set (Weiner et al., 1984, 1986), OPLS parameters (Jorgensen and Tirando-Rives, 1988), and our own parameters for unsaturated carbons are used to model the POPC lipid. Carbon atoms with hydrogen atoms are modeled as united-atoms, which totally reduces the number of atoms of one lipid molecule from 134 for an all-atom model to 52 for the united-atom model. The partial charges of the headgroup are adapted from previous studies using AMBER (Essmann et al., 1995a), which were calculated for DPPC lipids using the Gaussian 92 program (Frisch et al., 1992) and CHELPG module (Breneman and Wiberg, 1990) at 3-21 G* and 6-31 G* levels. Other values for the partial charges have also been documented in the literature for GROMOS (Egberts et al., 1994) and for CHARMM (Heller et al., 1993; Feller et al., 1997). An initial configuration of 200 equilibrated POPC lipids in liquid-crystalline phase with 5483 TIP3P (Jorgensen et al., 1983) water molecules was also adapted from previous work (Heller et al., 1993) and was equilibrated again for the change of force field parameters by using a simulated annealing technique. Finally, the system was brought to the desired temperature of 300 K. After a 30-ps simulated annealing simulation, the system reached roughly constant box dimensions (75.1 Å × 55.6 Å × 90.6 Å). Particle-Mesh-Ewald method (Essmann et al., 1995b) was then used for a further 170-ps simulation, while in previous simulated annealing simulation a 9 Å cutoff was used for the Coulomb interaction calculation. Inasmuch as the width of the water layer was ~23.4 Å, which was too narrow for putting a peptide like melittin in between, we decided to add more water molecules to our system. The preparation was done by LEaP (Schafmeister et al., 1995) and our own program. Totally, there were 10,951 water molecules in the system with an ~46.7 Å water layer. The new system with 43,253 atoms was run for 50 ps to reach a new equilibrium state. The physical properties of the lipids are in very good agreement with experiments (Lin and Baumgaertner, 2000).

Initial pore configuration

The initial structure of melittin was taken from the Protein Data Bank (PDB code: 2 mlt). Melittin is 26 amino acids long, including the C-terminus and the N-terminus. It has the sequence (H₂N-Gly*-Ile-Gly-Ala-Val-Leu-Lys*-Val-Leu-Thr-Thr-Gly-Leu-Pro-Ala-Leu-Ile-Ser-Trp-Ile- Lys*-Arg*-Lys*-Arg*-Gln-Gln-CONH₂), where charged amino acids are indicated by an asterisk. The N-terminus part (Ile-Gly-Ala-Val-Leu) is more hydrophobic, whereas the anchor sequence (Lys*-Arg*-Lys*-Arg*-Gln-Gln) is strongly hydrophilic.

A four-helix bundle with its symmetry axis parallel to the membrane normal was constructed (Vogel and Jähnig, 1986; Dempsey, 1990; Sansom, 1991; Bechinger, 1997). The bundle was made compact, but van der Waals contacts were kept minimal. In addition, the helices were rotated about their helical axes until their hydrophilic sides were facing each other. This orientation was anticipated because the hydrophilic side chains must be expected to be in contact with the water molecules inside the pore. Because the hydrophilic side chains are facing water molecules, we keep the charged residues in their protonated state, unlike the situation of a single melittin buried completely inside a lipid membrane (Bernéche et al., 1998). Following procedures similar to previous simulations of polyalanine (Shen et al., 1997), melittin (Bernéche et al., 1998), and porin (Tieleman et al., 1999) a hole 35 Å in diameter in the POPC membrane bilayer was made to insert the melittin bundle. Finally, a cylindrical box 9 Å in diameter containing water molecules was inserted into the bundle.

In order to bring the initial conformation closer to an equilibrium state, and to avoid some spurious artifacts due to unequilibrated lipids and water conformations, a combined Monte Carlo/molecular dynamics technique was used. A similar approach using Monte Carlo methods applied to the construction of the initial configuration of the gramicidin A channel embedded in a fully hydrated DMPC (dimyristoyl phosphatidylcholine) bilayer has been reported recently by Woolf and Roux (1996). The details of our approach are the following. We used standard Monte Carlo methods and the AMBER force field to relax the positions of the water and lipid molecules, while the melittin pore was kept immobile. The motion of water molecules was first performed randomly in a sphere centered at the chosen water molecule with radius 0.15 Å, and then rotated along some randomly chosen axis for random angles with a maximum value of 45°. For lipid molecules the moves are composed of two parts, first horizontal translation in the xy plane, and then a vertical translation along the z axis. The maximum displacement in the x and y direction was 0.15 Å, and in the z direction, 0.3 Å. One Monte Carlo step is defined as the scan of trials for all water and lipid molecules; totally, 2000 Monte Carlo steps were used to equilibrate the system. The acceptance ratio was ~0.23. The molecular dynamics run then followed the Monte Carlo simulation. After a 50-ps simulation, the system reached roughly constant box dimensions (84.9 Å × 85.9 Å × 71.6 Å). Particle-Mesh-Ewald method was then used for further a 400-ps simulation, and four lipids in the top layer were removed to make the number of lipids in both layers equal. It should be noted that 4 × 6 Cl counterions were solvated among water molecules. The final model system consisted of 4 melittin peptides, 152 POPC lipids, 24 Cl ions, and 11,133 water molecules, which corresponds in total to 43,071 atoms. The production run was then performed for 5.8 ns.

Computational details

The Monte Carlo program developed in our group at MOD of Forschungszentrum Jülich using the AMBER force field was used to prepare the initial conformation for the molecular dynamics simulations. Molecular dynamic simulations were performed by SANDER in AMBER 5.0 (Case et al., 1997) installed on the CRAY T3E at the Forschungszentrum Jülich. Berendsen thermostat and barostat (Berendsen et al., 1984) were used to keep the system in the specified temperature and constant pressure (1 Bar). Solute and solvent atoms were coupled independently to the thermostat with the same coupling constant 0.2 ps, and the center of mass motion was removed at each picosecond, whereby we removed the artifact due to the velocity rescaling scheme (Zhang et al., 1995; Harvey et al., 1998). An isotropic pressure scaling was used, which means the trace of the instantaneous pressure tensor P_iso = (P_xx + P_yy + P_zz) was used to compare with the reference pressure, and which allows uniform expansions and contractions of the system volume. A coupling constant for the barostat _p = 0.1 ps was used. The average fluctuations of the box dimensions are found to be very small, in the order of 0.5 Å. In such case the ratio of the surface area to system volume serves as the implicit constraint (Zhang et al., 1995). During MD simulations the SHAKE algorithm (Ryckaert et al., 1977) was used, so the hydrogen atoms did not have the bond stretching freedom, which allowed us to use the larger 2-fs time step. The order of B-spline interpolation for Particle-Mesh Ewald was set to 4, which implies a cubic spline approximation. The direct sum tolerance was set to be 0.00001. The scale factors for 1-4 electrostatic interactions (SCEE) and for 1-4 van der Waals interactions (SCNB) were both set to 2.0 (Weiner et al., 1984, 1986; Case et al., 1997). Because explicit water molecules were included in the simulation as solvent molecules, no distance-dependent dielectric constant was used. The atomic coordinates were saved every 1 ps and the atomic velocities were saved every 10 ps to reduce the cost for the need to rerun some part of the simulation. It took ~0.16 h for a 1-ps run on CRAY T3E using 32 PEs.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODELS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS AND REMARKS REFERENCES

Conformational properties of melittin

The end-to-end distance of melittin from the crystal structure (Terwilliger et al., 1982) is ~35.87 Å, which shows that melittin is a transmembrane peptide. Melittin is a comparably short protein which would serve perfectly as a single transmembrane-spanning -helix. In fact, melittin is expected to have a predominantly helical structure with a kink approximately at position 14 caused by proline which does not form hydrogen bonds (Heijne, 1991).

In our simulations we found that the conformation of melittin is remarkably stable over the whole time scale of 5.8 ns. The snapshot of the four melittin molecules in their final state is presented in Fig. 1. The snapshot indicates that melittin has a bent or open hairpin-like conformation. The hairpin is formed by two -helices connected to each other by a kink induced by Pro-14. The bending angle (t) enclosed by the two segments of each helix is presented as a function of time t in Fig. 2. One observes that the bending angle is approximately in the range of 90° <  < 160° and seems to converge for all of the four helices to   134 ± 20°.

View larger version (51K): [in this window] [in a new window]   FIGURE 1   Snapshot of a melittin pore. Prepared by MolScript v2.1 (Kraulis, 1991). The N-terminus and the C-terminus are denoted by N and C, respectively. The four helices are labeled 1-4.

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Bending angle (t) as a function of time. The fluctuations of (t) have been smoothed by averaging over 100 ps.

The -helicity of the peptide, h(t), is determined by the O(i) to H-N(i + 4) distances. According to the AMBER '91 force-field hydrogen bond parameters for these two atom types, we defined the hydrogen bond length as 2.5 Å. Thus a hydrogen bond between residues i and i + 4 is formed when the O(i) to H-N(i + 4) distance is 2.5 Å. Because there are 26 residues in melittin, a perfect -helix of melittin could form 22 hydrogen bonds; therefore, we can define the helical order parameter h of melittin as the number of hydrogen bonds in this peptide divided by 22. The helicity of melittin as a function of time, h(t), is presented in Fig. 3. It is found that the average helicity of melittin is h  0.75, and does not change significantly as a function of time. Because h values for each of the monomers differs considerably, one may conclude either that the typical correlation time of h(t) is larger than our simulation time of 5.8 ns, or that there is some structural correlation we couldn't identify.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Helicity h(t) as a function of time. The fluctuations of h(t) have been smoothed by averaging over 100 ps. The numbers 1-4 correspond to the labels given in Fig. 1.

The orientation of each of the helices with respect to the normal of the membrane surface has been estimated from our MD results. The orientational angle  between the end-to-end vector and the surface normal is shown in Fig. 4 as a function of time. From the data shown in Fig. 4 it can be concluded that the orientation of the helices relax to a stable orientation with respect to the membrane surface, which is of the order of   150 ± 10°. A similar tilt angle has been observed in other helix aggregates, e.g., in the case of the transmembrane four-helix bundle of influenza A M2 protein channel (Kovacs and Cross, 1997).

View larger version (19K): [in this window] [in a new window]   FIGURE 4   The orientation (t) of melittin with respect to the membrane normal as a function of time t.

Configuration of the melittin aggregate

The initial and the final locations of the four -helices in the xy plane after 5.8 ns are depicted in Fig. 5 by the dotted and the full circles, respectively. Each circle represents the projection of an ideal helix of diameter 12 Å located at their center of mass. The box size, as depicted in Fig. 5, has the actual xy dimensions of the periodic cell at the time of 5.8 ns. Comparing the locations of the initial and the final places, one finds that the initial tetrameric state is unstable and has changed after 5.8 ns to a trimeric state, consisting of helices 1-3, plus one monomer, helix 4. The transition from the tetramer to the trimer plus monomer is continuous. This can be concluded from the time evolution of the distances R_ij(t) between pairs (i, j) of helices. This is presented in Fig. 6. The interhelical distances R_ij are defined between the center of mass of the helices in the xy plane. The data clearly indicate that helices 1-3 remain in close contact with each other (Fig. 6, bottom) forming a stable trimeric aggregate where the distances between the pairs (1, 2) and (2, 3) are approximately the same, 20 Å, during the whole simulation. Because the distance between helices 1 and 3 is ~33 Å, the configuration of the trimeric bundle corresponds to an equilateral triangle in the xy plane. In contrast, helix 4 becomes continuously separated from the trimer up to a distance of ~50 Å after 5.8 ns (Fig. 6, top).

View larger version (17K): [in this window] [in a new window]   FIGURE 5   Initial and final configuration of the four melittin molecules in the xy plane.

View larger version (31K): [in this window] [in a new window]   FIGURE 6   Distances Rij(t) between pairs (i, j) of melittin helices versus time.

One interesting question is concerned with the internal motion of the trimer. As these three helices are in proximity to each other, in particular the consecutive pairs (1, 2) and (2, 3), some type of cooperative motion can be expected. The center-of-mass motions, as shown in Fig. 6, do not exhibit some indications related to cooperativity. However, some kind of collective rotation is still possible. One simple way to detect such an effect is to analyze, instead of the center-of-mass motion, the relative displacements of the C-termini and the N-termini of the helices. If R_ij^C and R_ij^N are the distances between the C-termini and the N-termini of helices i and j, respectively, then R_ij = R_ij^C  R_ij^N gives the difference between the separations of the two C-termini and the two N-termini. The quantity R_ij can be considered to measure the relative orientation of helices i and j with respect to each other. The "collective orientation" R_ij(t) as a function of time is shown in Fig. 7. As far as we can conclude from a simulation of 5.8 ns, it seems conceivable that the trimer performs a collective motion between two states, probably of an oscillating form. In one state the helices are parallel and hence R_ij very small, whereas in the other state the helices are tilted with respected to each other and R_ij is large. According to Fig. 7 (top) the tilted configuration is very specific and probably of a form as depicted in the cartoon in the bottom part of Fig. 7: in the case of pair (1, 2) the distance between the N-termini is larger than the distance between the C-termini, hence R_ij < 0; in the case of pair (2, 3) it is just the opposite, and therefore R_ij > 0. Of course, cooperative tilting of the three helices must not necessarily occur in one plane, as suggested by the two-dimensional cartoon in Fig. 7 (bottom), but rotation with respect to the xy plane must be expected. The correlation time between the two states seems to be of the order of 2 ns, which may be considered as a lower bound. It should be noted that during the whole simulation the residue Lys-7 remains directed toward the aqueous part of the pore.

View larger version (34K): [in this window] [in a new window]   FIGURE 7   Collective orientation Rij(t) between pairs (i, j) of the melittin trimer versus time.

Pore formation supported by lipids

A typical snapshot of the aqueous part of the pore in its final stage (5.8 ns) is presented in Fig. 8. From these four pictures, taken from four different directions in the xy plane, one can estimate the effective pore diameter d_P, which is in the range of 28 < d_P < 35 Å. This is in reasonable agreement with estimates obtained by efflux experiments of fluorescent molecules (Rex, 1996; Matsuzaki et al., 1997; Ladokhin et al., 1997). Because the radius of the pore increases with time, the number of water molecules inside the pore must increase accordingly. We have estimated the number N_W(t) of water molecules inside the pore as a function of time. In order to avoid end effects near the lipid head regions the length of the pore was restricted to 17 Å. The result for N_W(t) is shown in Fig. 9. Some important observations, which may provide a basis for an explanation of the expansion of the pore and the concomitant transition of the tetrameric melittin configuration to a trimeric one, are related to the behavior of the lipids as a function of time.

View larger version (119K): [in this window] [in a new window]   FIGURE 8   Snapshot of the distribution of water. For the sake of clarity we have removed lipids and melittin.

View larger version (16K): [in this window] [in a new window]   FIGURE 9   Number of water molecules, NW(t), inside the pore as a function of time.

Because the diameter of the aqueous pore increases with time, its cylindrical surface increases accordingly. Also, as the distance between the trimeric and the monomeric melittin increases with time, the initial shielding of water from direct contact with the acyl chains by the initial tetrameric melittin becomes less and less effective as the distance between the trimer and the monomer increases. Because the hydrophobic parts of the lipids are exposed to water, the interesting question arises how this conflicting situation is resolved by the system.

First we have identified those lipids which are close to the water molecules of the pore. The number of these pore-forming lipids, N_L(t), at the rim of the pore is shown in Fig. 10 and is in agreement with the idea of an expanding pore: N_L(t) increases with time.

View larger version (21K): [in this window] [in a new window]   FIGURE 10   Number of lipid molecules, NL(t), participating in pore formation as a function of time.

An analysis of these pore-forming lipids had revealed that during the simulations some lipids had changed the orientation of their acyl chains with respect to the membrane normal from a parallel to a perpendicular orientation. This reorientation is accompanied by a translocation of the lipid heads from the membrane surface to the surface of the pore. This kind of translocation of some lipid heads is shown in Fig. 11, where we have plotted the z coordinates of all 152 PO₄ groups of the lipid heads as a function of time. Initially, all coordinates are localized at z  52 ± 3 Å and z  25 ± 3 Å, which is expected for a perfect bilayer membrane. After ~1 ns, however, some lipid heads (denoted by the solid lines in Fig. 11) start to move to the central core of the bilayer. Finally, at 5.8 ns we have identified eight lipid heads (solid lines, Fig. 11) in the range 45 < z < 33 Å. Most importantly, these lipid heads are intercalated between the melittin trimer and the monomer, which is shown in Fig. 12 by a top view and a side view at 5.8 ns, including some water molecules, some lipid molecules, and the melittin molecules. The flipped lipids are depicted separately. No lipid head intercalated between the melittins of the trimeric aggregate have been found during the whole simulation. It should be noted that this type of lipid reorientation and translocation is a phenomenon well-known from the structures of "hydrophilic" membrane pores purely formed by lipids as, e.g., in electroporation phenomena (Neumann et al., 1989; Chang et al., 1992).

View larger version (70K): [in this window] [in a new window]   FIGURE 11   Trajectories of the z-coordinates of all lipid heads as a function of time. For the sake of clarity, only those trajectories of lipid heads exhibiting a transition from the membrane surface to the inner part of the pore have been depicted by solid lines; trajectories of other lipid heads are given by dotted lines.

View larger version (67K): [in this window] [in a new window]   FIGURE 12   Snapshot of the pore from top (A) and from side (B). Depicted are the four melittins, the water inside the pore, and lipids, including the eight flipped lipids. The phosphorus atoms of the flipped lipids are shown by the space-filling model (orange). Figures were prepared with RasMol 2.6 (R. Sayle).

Similar conjectures related to the intercalation of lipids between helices of a pore have been reported for magainin (Matsuzaki et al., 1996, 1997a; Ludtke et al., 1996). Such a type of pore has been termed a "toroidal" (wormhole) model (Ludtke et al., 1996).

The important and interesting point, however, is the intercalation. It is conceivable that this spontaneous intercalation, once it has happened, leads to a destruction of a stable n-meric aggregate. As a necessary prerequisite for this event to take place, we must postulate that sufficiently large fluctuations of the interhelical distances should occur. This hypothesis would imply that, ultimately, even the trimer should decay into monomers after a sufficiently long time, much longer than 5.8 ns.

If the intercalation of some lipids is a secondary effect, then the question arises why the trimeric state had remained stable during the whole simulation. A second question is concerned with the increasing separation between trimer and monomer, which is probably related to the electrostatics of the interhelical interactions. These two questions are discussed in the next two sections.

Stability of the trimer

It is conceivable that the stability of the trimer is related to the so-called hydrophobic effect. More precisely, the unfavorable contact between water and the acyl chains may induce an attractive interaction between the helices if the helices could expose their hydrophilic side chains to the water molecules inside the pore.

The proof of the validity of this hypothesis is quite delicate and difficult. As an attempt to support this hypothesis we have calculated the "shielding" probability P(t) of water at the boundary of the pore. We have calculated this quantity locally, i.e., limited to the neighborhood of the pairs of helices.

The "shielding" probability P(t) is defined as follows (compare also the cartoon in the bottom part of Fig. 13). Consider in the xy plane the midpoint r_i,i+1 = [r_i^cm + r_i+1^cm]/2 between the two center-of-mass coordinates, r_i^cm and r_i+1^cm, of one consecutive pair of helices, i and i + 1. This midpoint is located in the center of a cylinder along the z axis. The radius of the cylinder is r =  r_i^cm + r_i+1^cm /2 and its length is 11 Å. Counting the number of water molecules inside such a cylinder, this would provide a measure of the local concentration of water in that regime between the helices i and i + 1. In addition, if the cylinder is subdivided along the z axis into two semi-cylinders where the separation plane is defined according to the cartoon as depicted in the lower part of Fig. 13, then one can define the ratio P_o(t)/P_i(t) between the numbers of water molecules in the upper ("out") and the lower ("in") compartments. The results for the four consecutive pairs of helices are shown in the top part of Fig. 13. The data have been calculated at intervals of 100 ps. The "shielding" probabilities related to the trimer (pairs (1, 2) and (2, 3)) are much smaller than the curves related to helix 4 (pairs (3, 4) and (1, 4)). This indicates that the trimer can successfully shield the water from contact with lipids. Perfect shielding would require P_o(t)/P_i(t) = 0. The non-zero values for the trimer, however, can most probably be attributed to our inaccurate definition of locality defined by a rigid cylinder. For example, effects from orientational and conformational fluctuations of the helices are not taken into account. In conclusion, it is very likely that amphiphatic helices may form stable pores, at least for a certain lifetime.

View larger version (24K): [in this window] [in a new window]   FIGURE 13   Shielding probability ratio Po(t)/Pi(t) versus time (top). Construction of the local volume (bottom).

Interhelical electrostatic energies

Previous theoretical considerations on the helix-helix interactions in lipid bilayers (Ben-Tal and Honig, 1996) have indicated that the electrostatic interactions between two -helices can be quite strong in an alkane phase, because desolvation effects are essentially nonexistent and helix-helix interactions are not well screened. Therefore, antiparallel helix orientation in a helical bundle is strongly favored over the parallel orientation. The helix-helix interaction is typically of the order of a few kcal/mol.

In the present case, however, the helix-helix interaction is much larger due to the six charges on each melittin monomer. The electrostatic energies E_ij(t) between pairs (i, j) of melittin molecules versus time are presented in Fig. 14. The separation of the trimer from the monomer is accompanied by a decrease of the electrostatic energies E_ij(t) between the melittin molecules, which is shown in Fig. 14. All energies are given in units of kcal/mol. From the top part of Fig. 14 one finds that, in agreement with the stability of the trimeric state, the electrostatic energies between pairs of helices 1-3 fluctuate about some constant values, at least after 2 ns. This is in contrast to the electrostatic energies between the monomeric melittin, helix 4, and the helices of the trimer. In this case, which is shown in the bottom part of Fig. 14, the energies gradually decay with increasing time.

View larger version (41K): [in this window] [in a new window]   FIGURE 14   Electrostatic energy Eij(t) between pairs (i, j) of melittin helices versus time.

Of course, the separation of the trimer from the monomer cannot continue to a large extent because of the finite size of the system and, even more important, because of the repulsion between the helices and their images in the neighboring cells due to the periodic boundary conditions applied to the present system. Therefore, the separation must be expected to come to a rest at a certain size. This might then be an artifact due to the simulation technique.

	CONCLUSIONS AND REMARKS

TOP ABSTRACT INTRODUCTION MODELS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS AND REMARKS REFERENCES

One implication from the present study is concerned with the stability of aqueous helix bundles and related ion channels. Based on the present observation that the onset of pore growth and its concomitant reorganization become evident only after 2 ns, future simulations on related objects and phenomena have to address the question of stability with great care. Similarly, from the observed stability of the trimeric melittin aggregate, it is not possible to make any prediction of an average lifetime, but the implication is that at least for a few nanoseconds amphiphatic helices exist as stable aggregates, eventually even as water-filled helix bundles.

An interesting point concerning the driving force of the pore growth may be related to the competition between the surface tension and a hypothetical line tension created by the pore-lining lipids and the four melittin helices. The competition between surface and line tensions has been postulated to be responsible for the growth of pure lipid-formed ("hydrophobic") membrane pores which are observed, e.g., during electroporation. The mechanical contribution to a pore free energy is (Litster, 1975; Taupin et al., 1975; Weaver and Barnett, 1992)
where  is the surface energy density of the membrane-water interface and  is the pore edge energy. In the case of a planar membrane, growth is believed to occur if the pore achieves a radius r > r_c, where the critical radius for expansion is r_c = /. For radii smaller than r_c, the pore shrinks. Based on this concept and with regard to the present situation of a protein-lipid-composed pore, one may hypothesize that melittin initiates pore formation and creates, eventually together with lipids, a sufficiently large line tension leading to the onset of pore growth. The strength of the effective line tension induced by amphiphatic helices may differ in general from one protein to another. In the particular case of melittin, it is likely that the electrostatic repulsion between the anchor sequences (Lys*-Arg*-Lys*-Arg*-Gln-Gln) plays some role and eventually favors intercalation of lipids and water. Once this has happened the regular mechanism of pore growth may take over.

We have not examined in detail the properties of water inside the pore, as has been done in previous studies (Roux and Karplus, 1994; Breed et al., 1996). Of course, it would be of interest to compare our study with previous results. However, the present situation is not related to a stable equilibrated pore configuration, but rather to an expanding pore where the influx of water and the continuously changing boundary may spoil a reasonable analysis of our data.