INTRODUCTION

Top Abstract Introduction Methods Results & Discussion Conclusion References

Specific lipid-protein interactions involved in the anchoring and stabilization of membrane-bound proteins are of central importance in a large number of fundamental processes occurring at the surface of the cell. However, despite the development of powerful techniques such as x-ray crystallography (Deisenhofer and Michel, 1989), electron microscopy (Henderson et al., 1990), and nuclear magnetic resonance (NMR) (Cross and Opella, 1994), the characterization of lipid-protein interactions remains difficult because of the complexity of the bilayer environment. At the present time, even qualitative information gained by performing detailed computer simulations of protein-membrane complexes can be valuable, because only scarce information is available from experiments about the structure and dynamics of these systems. To begin to understand lipid-protein interactions at the microscopic level, we performed molecular dynamics simulations of melittin in a dimyristoyl phosphatidylcholine (DMPC) bilayer.

Melittin is the major protein component of the venom of the honey bee Apis mellifera that is responsible for lysis of the cell membrane (Habermann, 1972; Sessa et al., 1969). The cationic amphiphilic polypeptide consists of 26 amino acids with the sequence (Habermann and Jentsch, 1967) 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²⁶. Its structure has been determined to a high resolution by x-ray crystallography (Terwilliger and Eisenberg, 1982a) and NMR spectroscopy. The NMR structures were determined in detergent micelles (Inagaki, 1989) and in nonpolar solvent (Bazzo et al., 1988), in which the protein is present as a monomeric form, whereas the crystals were grown from concentrated aqueous salt solutions in which melittin molecules form closely associated tetramers through hydrophobic contacts (Terwilliger and Eisenberg, 1982a). Remarkably, these studies indicate that melittin adopts a very similar -helical conformation in those very different environments. The proline residue at position 14 is responsible for a bend separating two segments of the -helical structure: a hydrophobic segment, going from residue Gly¹ to Leu¹³, and an amphiphilic segment, going from residue Ala¹⁵ to Gln²⁶. Because no significant structural variations are observed in the different environments, it is likely that the -helical structure is a good model for the membrane-bound conformation. Because it has been studied extensively, melittin may serve as an interesting prototypical model of membrane-binding amphiphilic polypeptides (Segrest et al., 1990). More generally, amphipathic -helices represent an important structural motif that plays a role in the anchoring of monotopic membrane proteins (Picot et al., 1994).

The association of melittin with bilayers and the mechanism involved in the initiation of membrane lysis have been investigated by various approaches (Brown et al., 1982; Dawson et al., 1978; Dempsey, 1988; Dempsey and Butler, 1992; Frey and Tamm, 1991; Smith et al., 1992; Vogel and Jahnig, 1983; Vogel et al., 1986). The large number of studies shows that the interaction of melittin with membranes depends on the lipid composition, the peptide concentration, the hydration level, and the membrane potential (Dempsey, 1990). The orientation of melittin in phospholipid bilayer was shown to be dependent on pH, suggesting that the protonation state of the N-terminus of melittin influences its interaction with membranes (Bradshaw et al., 1994). The proline in position 14 as well as the polar residues 23-26 at the C-terminus have been shown to be essential for the lysis activity (Otoda et al., 1992; Rivett et al., 1996; Werkmeister et al., 1993). Moreover, bilayers consisting of lipids with longer acyl chains are less affected by the lysis activity of melittin, illustrating the importance of the membrane composition (Bradrick et al., 1995). It was also shown that negatively charged membranes are less affected by melittin, suggesting that the peptide binds strongly at the membrane surface without penetrating the bilayer (Ohki et al., 1994). Moreover, phosphorus NMR spectroscopy showed that melittin penetrates more deeply into the bilayer of zwitterionic lipids, causing membrane lysis (Monette and Lafleur, 1995).

On the basis of those observations, different mechanisms for the lytic activity of melittin have been proposed. One of the earliest hypotheses suggested that melittin increases membrane permeability by partial penetration of the bilayer (Ash et al., 1978). It was estimated that 250 melittin molecules per unilamellar palmitoyloleoyl phosphatidylcholine (POPC) vesicle of 100 nm diameter (i.e., ~200 lipids/monomer) were required to initiate the lysis, suggesting that this process is due to a collective membrane perturbation by monomeric bound peptides (Benachir and Lafleur, 1995). In contrast, some models involving the formation of a canal structure by the aggregation of four transbilayer melittin molecules have also been proposed (Smith et al., 1994; Tosteson and Tosteson, 1981). Others suggested that aggregated melittin is involved in the solubilization of large lipid disks, which would lead to cell lysis by leaving large holes in the membrane (Dufourc et al., 1986). It has also been proposed that the binding of melittin to other membrane proteins is involved in the initiation of the lytic mechanism (Portlock et al., 1990; Werkmeister et al., 1993). Clearly, a characterization at the molecular level of the association of melittin with membranes is necessary for a better understanding of the microscopic factors playing an important role in its lytic activity. Computer simulations of detailed atomic models represent a powerful approach to understanding such complex systems. In recent years it has been used to gain insight into the structure and dynamics of pure lipid membranes (Berger et al., 1997; Chiu et al., 1995; Pastor, 1994; Venable, 1993) as well as their interations with small solutes (Stouch, 1993) and proteins (Edholm et al., 1995; Huang and Loew, 1995; Woolf and Roux, 1994a, 1996; Shen et al., 1997; Tieleman and Berendsen, 1998; see also Merz and Roux, 1996, and references therein). In particular, one dynamical simulation of an amphipathic helical peptide bound to the surface of a dioleoyl phosphatidylcholine (DOPC) bilayer bilayer has been reported (Huang and Loew, 1995).

In the present work we specifically address questions about the interaction between a single melittin monomer and a phospholipid bilayer. This paper reports the results of molecular dynamics trajectories for a fully hydrated DMPC-melittin system. The initial configuration of the system was assembled using a general approach developed previously for constructing the starting configuration for molecular dynamics simulations of membrane-bound proteins (Woolf and Roux, 1994a, 1996). All atoms are explicitly included in the calculations. Melittin was assumed to be in the -helical conformation taken from the x-ray crystallographic structure (Terwilliger and Eisenberg, 1982a). In accord with its amphipathic amino acid sequence, the helix was oriented parallel to the membrane-solution interface such that the apolar residues are facing the hydrophobic core of the membrane and the polar residues are facing the water bulk phase. This orientation, which corresponds to the "wedge" model that has been described in the literature (Dawson et al., 1978; Terwilliger et al., 1982), is consistent with experimental observations of the orientation of melittin in bilayers (Bradshaw et al., 1994; Citra and Axelsen, 1996; Dempsey and Butler, 1992). Two protonation states of the N-terminus of melittin were simulated. For the first 600 ps of simulation, the N-terminus of melittin was unprotonated. At 600 ps, the N-terminus was protonated and the trajectory was pursued for 400 ps. As a control, the trajectory of the unprotonated system was also continued for 400 ps. In the next section, the theoretical methods and the atomic models are described in detail. The results are then given and discussed. The paper concludes with a brief summary of the main results.

	THEORY AND METHODS

Top Abstract Introduction Methods Results & Discussion Conclusion References

Construction of the microscopic model

The microscopic system consists of one melittin monomer (N-terminus unprotonated), 41 DMPC (17 in the upper layer, which contains melittin, and 24 in the lower layer), and 1595 water molecules, for a total of 10,056 atoms. This constitutes the central unit of a periodic system, the dimensions of which are 32 × 48 × 60 ų. The net charge of the protein with unprotonated N-terminus is +5e. The membrane normal is oriented along the z axis, and the center of the bilayer is at z = 0. Periodic rectangular boundary conditions were applied in the xy directions to simulate an infinite planar layer and in the z direction to simulate a multilayer system. In the simulation cell, melittin is deeply inserted into the upper monolayer, with its helical axis roughly parallel to the y axis. To avoid possible difficulties with constant pressure algorithms due to the significant anisotropy and inhomogeneity of the simulation cell, the trajectory was calculated in the microcanonical ensemble with constant number of particles, energy, and volume (NVE). The average temperature of the system was set to 330 K, above the gel-liquid phase transition of DMPC (Gennis, 1989). The potential energy function used for the calculations was the all-hydrogen PARAM 22 force field (MacKerell et al., 1998) of the biomolecular CHARMM program (Brooks et al., 1983), which includes phospholipids (Schlenkrich et al., 1996) and the TIP3P water potential (Jorgensen et al., 1983).

The construction of such a complex system requires careful attention to several details to get a meaningful trajectory. Because the simulations are computationally very intensive, it is desirable to build a starting configuration that is as representative of the solvated protein-membrane system as possible, thereby limiting the required equilibration time. A general protocol developed by Woolf and Roux (1994a, 1996) was used to construct the initial configuration of the protein-membrane system. This method has been used previously to generate configurations for the gramicidin channel (Woolf and Roux, 1994a, 1996) and bacteriophage Pf1 coat protein in lipid bilayers (Roux and Woolf, 1996). The general strategy for creating a representative starting configuration for the system consists of randomly selecting lipids from a preequilibrated and prehydrated set, dispersing them around the protein, and reducing the number of core-core overlaps between heavy atoms through systematic rigid-body rotations (around the z axis) and translations (in the xy plane) of the preequilibrated and prehydrated lipid molecules. The initial melittin configuration was taken from the 2.0-Å resolution x-ray crystallographic structure (Terwilliger and Eisenberg, 1982a). The helix was oriented parallel to the membrane-solution interface, in a position that left the hydrophilic residues in contact with the bulk water and the hydrophobic residues in contact with lipid acyl chains.

The total cross-sectional area for the simulation of the protein-membrane system must be carefully determined because it has an important influence on the state of the bilayer (Chiu et al., 1995; Heller et al., 1993; Woolf and Roux, 1996). The system shown in Fig. 1 has a marked asymmetry because of the amphipathic segment lying parallel to the membrane surface in the upper layer. It is necessary to account for the cross-sectional area of melittin to determine the appropriate number of lipids to include in the upper and lower halves of the bilayer at the microscopic model. Although the present simulations were carried out at constant volume, a reasonable estimate of the cross-sectional area of the system and the asymmetrical number of lipids in the upper and lower halves of the bilayer is also required, even with constant pressure algorithms in which the membrane cross-sectional area is allowed to vary (Berger et al., 1997; Chiu et al., 1995; Feller and Pastor, 1996). The cross-sectional area of melittin in the membrane-bound conformation is 450 Å². To surround the helix by a complete lipid environment, the dimension of the system in the xy plane was set at 32 Å × 48 Å, corresponding to an area of 1536 Å². Because the average cross-sectional area of a single DMPC molecule is 64 Å² (Gennis, 1989; Nagle, 1993; Pastor et al., 1991), this corresponds to a total of 24 DMPC molecules (lower layer), or to one melittin and 17 DMPC molecules (upper layer). To determine the initial position of each lipid, the DMPC polar heads were first represented by large effective spheres with a cross-sectional area of 64 Å². The position of the large spheres was obtained through molecular dynamics and energy minimization with the same periodic boundary conditions as those used in the simulation of the complete system. The spheres were constrained at z = 17 Å and z = 17 Å for the upper and lower layers, respectively. The resulting configuration is shown in Fig. 1. The effective lipid spheres were then substituted by full DMPC molecules randomly chosen from a library of 2000 preequilibrated phospholipids (Venable et al., 1993). In this library, the polar headgroups of the DMPC are prehydrated by ~20 water molecules constructed on the basis of a molecular dynamics simulation of o-phosphorylcholine (o-PC) in bulk solution (Woolf and Roux, 1994b). The xyz coordinates of the different spheres were used to position the center of mass of the phosphorus and nitrogen atoms of the lipid polar heads. Systematic rigid-body rotations of the lipids around the z axis and translations in the xy plane were then performed to minimize the number of unfavorable contacts and atomic overlaps. The remaining bad contacts were removed by energy minimization. The system was then fully hydrated by overlaying a preequilibrated water box of the appropriate dimension in x and y. The resulting configuration is shown in Fig. 2 A. The system was further refined by energy minimization before the dynamical simulation was started.

View larger version (86K): [in this window] [in a new window]   FIGURE 1   Configuration of the large effective lipid spheres on the two layers: (top) view from the side; (bottom) view from above. The rectangle of 32 Å × 48 Å indicates the central unit of the system. Molecules outside the central box represent periodic images of the molecules inside the central unit. The -helical conformation was taken from the crystallographic structure (Terwilliger and Eisenberg, 1982).

View larger version (79K): [in this window] [in a new window]   FIGURE 2   Representation of the atomic system through the equilibration period. The deformation of the membrane in the lower layer is significant (it should be noted that only the central system is shown and that no significant vacuum exists because of the periodic boundary conditions). The system shown at 150 ps corresponds to the beginning of the production trajectory.

Computational details: equilibration and dynamics procedures

The system was equilibrated for 150 ps by molecular dynamics. Fig. 2, A-D shows snapshots of the system through the equilibration period. To converge to an equilibrium state, the system was coupled for the first 125 ps to a heat bath at 330 K by the use of Langevin dynamics. During the last 25 ps of equilibration the velocities were periodically rescaled to stabilize the temperature. The equations of motion were integrated with a time step of 2 fs, and the coordinates were saved every 0.1 ps. The list of nonbonded interactions was truncated at 12 Å by the use of a group-based cutoff. The nonbonded van der Waals and electrostatic interactions were smoothly switched off over a distance of 3.0 Å, the values being maximum for a radius of less than 8 Å and zero at 11 Å. The SHAKE algorithm (Ryckaert et al., 1977) was used to fix the length of all bonds involving hydrogen atoms. A number of energy restraints were used at the beginning of the equilibration period to ensure a smooth relaxation of the system toward an equilibrated configuration. Harmonic potentials were applied to the melittin backbone to prevent large spurious motions, the center of mass of the lipid polar heads was kept around z = ±17 Å by planar harmonic constraints to maintain the planarity of the membrane, and the penetration of water in the bilayer region was prevented within z by the use of planar potentials. All of those restraints were gradually reduced to get a completely free system after 100 ps of equilibration. During the production trajectory, the center of mass of the protein was restrained to the center of the xy plane. For the first 600 ps of simulation, the N-terminus of melittin was unprotonated. At 600 ps, the N-terminus was protonated and, after 500 steps of energy minimization, the trajectory was continued for 400 ps. As a control, the trajectory of the system with an unprotonated N-terminus was also continued for 400 ps. Detailed views of the system after 600 ps (unprotonated N-terminus) and 850 ps (protonated N-terminus) are shown in Fig. 3, A and B.

View larger version (120K): [in this window] [in a new window]   FIGURE 3   (A) Atomic system after 600 ps of dynamics. The hydrogens atoms of the acyl chains are not shown for the sake of clarity. The presence of a small number of water molecules near the unprotonated N-terminus of melittin and the deformation of the bilayer is observed. (B) Atomic melittin system protonated N-terminus after 850 ps of dynamics. A large number of water molecules forming a continuous network across the membrane is observed.

Calculation of solid-state NMR and polarized FTIR properties

Solid-state NMR spectroscopy (Cross and Opella, 1994) and polarized attenuated total internal reflectance Fourier transform infrared spectroscopy (PATIR-FTIR) (Axelsen et al., 1995) are methods of increasing importance for studying the structure and orientation of membrane-bound peptides and proteins. However, a direct structural interpretation of the measurements obtained from those techniques in terms of a peptide conformation and orientation is often not straightforward. For this reason, it is of interest to examine the relationship of the observed data with the average structure and orientation on the basis of a molecular dynamics trajectory. The time scale of solid-state NMR is much slower than that of rapid molecular motions (Cross and Opella, 1985; Seelig and Seelig, 1980). For this reason, observed properties such as the chemical shift and the deuterium quadrupolar splitting (DQS) correspond to a time average over rapidly fluctuating quantities. To obtain the observed properties, an average over instantaneous values must be performed. For example, the chemical shift for a specifically ¹⁵N-labeled site is given by a time average over a large number of configurations of the projection of the instantaneous second-rank shielding tensor (Cross and Opella, 1985; Woolf and Roux, 1994a), i.e.,

(1)

where _ii(t) and ê_i(t) are, respectively, the instantaneous magnitude and direction of the principal tensor components and is a unit vector in the direction of the bilayer normal. Similarly, the DQS order parameter, S_CD, for a specifically deuterated site is (Seelig and Seelig, 1980; Woolf and Roux, 1994a)

(2)

where (t) is the instantaneous angle between the director of the C-D bond and the bilayer normal. The magnitude and orientation of the three components of the ¹⁵N backbone chemical shift tensor of polypeptides have been determined experimentally from powder spectra (Cross and Opella, 1985; Teng and Cross, 1989). Typically, the largest component of the tensor, ₃₃, has a magnitude of 201 ppm and an orientation approximately parallel to the N-H bond in the amide plane. The component ₂₂ (perpendicular to the amide plane) has a magnitude of 55 ppm, and the component ₁₁ (in the amide plane) a magnitude of 28 ppm (Teng and Cross, 1989). To compute the instantaneous chemical shift for a backbone site, the tensor components were built in the local molecular frame on the basis of the atomic configurations taken from the trajectory. The principal axis of ₃₃ (201 ppm) was constructed in the H-N-C plane with an angle of 105° relative to the N-C bond. The principal axis of ₂₂ (55 ppm) was constructed perpendicular to the H-N-C plane, and that of ₁₁ (28 ppm) was obtained from a cross product of the second and third principal axes (i.e., ê₁ = ê₂ × ê₃). This rigid tensor approximation ignores the rapid fluctuations of the tensor component magnitudes caused by variations in the local backbone geometry (Woolf et al., 1995).

PATIR-FTIR measurements can also be used to determine the orientation of membrane-bound helical peptides (Axelsen et al., 1995). In this method, the absorption coefficients parallel and perpendicular to the interface are measured and the orientational helical order is determined from the amide I dichroic ratio. Because the time scale of FTIR is much faster than that of the slow reorientational movements, the observed orientation corresponds to an ensemble average. The observed order parameter S_PATIR results from a superposition of the amide I band of all the residues,

(3)

where _i is the angle between the amide I transition moment and the normal to the interface for the ith residue. A structural interpretation of the experiments requires information about the orientation of the amide I transition moment with respect to the peptide backbone. Although the amide I vibration arises from concerted displacements within the N-H-C-O group (Krimm and Reisdorf, 1994), in the present treatment we assumed that the transition moment is parallel to the C==O backbone carbonyl bond for the sake of simplicity.

Continuum model for protein-membrane association

A mean-field potential based on a continuum electrostatic approximation was used to investigate the importance of thermodynamic and energetic factors in the membrane association of melittin. A similar approach has been used by Ben-Tal et al. (1996) to examine the association of a polyalanine -helix with a membrane. According to this approximation, the total free energy of solvation G_tot is decomposed into an electrostatic contribution G_elec and a nonpolar cavity formation G_np (Ben-Tal et al., 1996; Gilson and Honig, 1988),

(4)

The term G_elec accounts for the reaction field contribution to the free energy of charging the peptide in the polar medium. The electrostatic contribution to the free energy of transfer from water to the membrane along the z axis was computed with the Poisson equation for macroscopic continuum electrostatics (Honig and Nicholls, 1995; Honig et al., 1993; Warwicker and Watson, 1982),

(5)

where (r) is the position-dependent dielectric constant and ^prot(r) is the protein charge density. The protein was represented at the microscopic level with its associated atomic radii and charges. The atomic radii used to define the protein-solvent dielectric interface were derived from radial distribution functions calculated for the 20 standard amino acids from molecular dynamics simulations and free energy perturbations with explicit water molecules (Nina et al., 1997). The atomic charges were taken from the all-hydrogen parameters PARAM22 (MacKerell et al., 1998). The membrane was represented by a planar slab 25 Å thick corresponding to the width of the hydrocarbon core of the membrane (White and Wiener, 1996). Dielectric constants were assigned according to the polarity of the medium:  = 80 for bulk water,  = 2 for the membrane, and  = 1 for the protein. Because of the uncertainty of the continuum description of the water-lipid interface forming a transitional dielectric region, the slab of low dielectric was defined to represent the hydrocarbon chain region only. The region corresponding to the polar headgroups was assumed to have a dielectric constant of 80. All calculations were performed with a cubic grid of 70 Å with two grid points per Å. The ionic strength was set to zero (no counterions were included). The geometrical center of the membrane was set to 10 Å along the z axis of a three-dimensional cubic grid. The protein was mapped onto the grid; the center of mass of the helix was placed initially at the geometrical center of the membrane. For each position of melittin, the electrostatic contribution to the free energy of transfer from the water to the membrane was calculated by subtracting the electrostatic energy computed in a continuous medium representative of water ( = 80) from the electrostatic energy computed in a membrane ( = 2) immersed in a solvent region ( = 80),

(6)

where q_i is the charge of the ith atom in the protein, ^memb(r_i) is the total electrostatic potential of the protein initially embedded in a membrane, and ^bulk(r_i) is the total electrostatic potential of the protein immersed in the bulk region. The electrostatic potential (r) was obtained by solving the Poisson equation with the finite-difference algorithm of Klapper et al. (1986), implemented in the PBEQ (Beglov and Roux, unpublished) facility of CHARMM (Brooks et al., 1983).

The term G_np accounts for the nonpolar contributions and is assumed to be related to the water-accessible surface of the peptide. A simple sum over the water-accessible surface of all atoms is used to approximate the nonpolar contributions,

(7)

where z_i is the z position of atom i (the membrane normal is oriented along the z axis), S_i^np is the water-accessible surface of the ith atom, and  = 33 cal/mol/Ų is the surface tension coefficient obtained from experimental free energies of transfer of hydrocarbons from the pure liquid alkane to water (Hermann, 1972). The position of the interface z₀ and its width z were set at 10.0 Å and 2.5 Å, respectively, according to experimental data on lipid bilayers (Jacobs and White, 1989; White and Wiener, 1996). The water-accessible surface S_i was calculated by rolling a probe of 1.4 Å on the protein van der Waals surface. Atomic radii were taken from the all-hydrogen PARM22 parameters file (MacKerell et al., 1998).

	RESULTS AND DISCUSSION

Top Abstract Introduction Methods Results & Discussion Conclusion References

Average structure of the system

Fig. 2 shows the transformations taking place in the membrane structure during the equilibration period. The initial bilayer configuration, which is perfectly planar as a result of the construction protocol, is progressively modified by vertical displacements of the lipids along the z axis. In addition, the order of the lipid acyl chains is decreasing in the upper layer while it is increasing for those in the lower layer. The distribution of the carbons of the acyl chains in the bottom layer is shifted to the top layer. As a consequence, the acyl chains are more extended in the lower layer at the end of the equilibration period relative to the initial configuration. The perturbation of the lipid configurations results in a local curvature and a reduction of the thickness of the membrane (see Fig. 2). Because of this structural reorganization, the hydrophobic core of the membrane is reduced by ~30% from its original thickness near the center of the system. However, no significant structural changes in the protein-membrane system are observed after the relatively short equilibration period of 150 ps. Test simulations with other starting lipid configurations resulted in very similar perturbations of the bilayer structure (not shown).

To characterize the dominant structural features of the melittin-membrane system, the average density profile of the main components was calculated. The first 600 ps of the trajectory in which the N-terminus of melittin is unprotonated was used. The result, shown in Fig. 4, indicates that melittin is located roughly between the lipid-water interface and the center of the membrane. The density of water molecules converges to the normal bulk density (0.033 molecules/ų) away from the membrane and decreases in the membrane region. Water molecules go as deep as 5 Å from the middle of the membrane. The phosphate and choline groups are both located around z = ±18 Å, which reflects the fact that the polar heads are oriented nearly parallel to the membrane-water interface. A broadening of the membrane structure was observed in a study of an amphipathic helix bound to the surface of a DOPC membrane (Huang and Loew, 1995). As indicated by the distribution of the headgroups, no significant difference in the width of the membrane-bulk interface is observed in the present simulation. The density of the hydrocarbon chains is slightly reduced around z = 2.5 Å. In a pure membrane the corresponding minimum is observed in the middle of the bilayer at z = 0 Å (Berger et al., 1997; Pastor, 1994; White and Wiener, 1996). The displacement of the local minimum from the center of the membrane toward the upper leaflet is due to the presence of melittin.

View larger version (21K): [in this window] [in a new window]   FIGURE 4   Average density profile of the main components of the melittin-DMPC system. The density profile of the heavy atoms of melittin, water, and the hydrocarbon chains, the ester oxygens of the glycerol region, and the headgroup phosphate and nitrogen of the DMPC lipid molecules are shown.

Stability of the protein and macroscopic motion

During the trajectory, the center of mass of melittin (with unprotonated N-terminus) drifted toward the center of the membrane, from 13.0 Å to 11.0 Å along the z axis. To characterize the global motions of melittin at the interface, its instantaneous configurations were reoriented relative to an average reference structure by minimizing the coordinate root mean squared (RMS) difference. The results indicate that the helix is undergoing rocking motions with fluctuations on the order of 5-8°. On average, the rotation axis was oriented parallel to the main helical axis. Nevertheless, as shown in Fig. 5, the backbone conformation did not changed significantly. The initial structure remained stable, although there are small deviations from a perfect -helical configuration near the C-terminus. This is in qualitative accord with circular dichroism (CD) data, which indicate that the helical content of membrane-bound melittin is 72% (Vogel, 1987). It is likely that the initial helical conformation near the C-terminus of the peptide is not conserved because its polar and charged residues, Lys²¹-Arg²²-Lys²³-Arg²⁴-Gln²⁵-Gln²⁶, are extensively solvated by water molecules (see below). In contrast, the helical conformation of the nonpolar segment of the peptide, which is embedded in the hydrocarbon core of the membrane, is very stable. The average angle between the two helical segments (defined by the two vectors joining the C of residues Val⁵ to Gly¹², and Leu¹⁶ to Lys²³, respectively) was ~145°, with fluctuations on the order of 5°.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   Superposition of five configurations of the backbone of melittin taken from the trajectory. The configurations are separated by 100 ps.

Experimental techniques, such as solid-state NMR and PATIR-FTIR, provide important information concerning the conformation and orientation of a membrane-bound peptide. However, interpretation of the data in terms of a unique microscopic structure may be difficult. Thus, even though the data are not yet available, it is of interest to examine how the conformation would be characterized by those experimental techniques on the basis of the current theoretical model. Fig. 6 A shows the average solid-state NMR backbone ¹⁵N chemical shifts, computed from the trajectory. There are large variations because of the helical periodicity. The residues of the amphipathic helical segment (15-26) have the lowest chemical shift (between 50 and 80 ppm), in accord with a parallel orientation relative to the membrane surface. The residues of the hydrophobic helical segment (1-14), which is in a diagonal orientation relative to the bilayer normal, have intermediate values (between 80 and 150 ppm). The C-terminal segment is disordered and deviates from the -helical structure. For this reason, the chemical shift of the residues at the C-terminus differs from the value obtained for the other residues of the amphiphilic helix. The calculated deuterium-C bond order parameters are shown in Fig. 6 B. As in the case of the ¹⁵N chemical shift, there are large variations because of the helical periodicity. The very high values observed for Ala⁴, Val⁸, and Ala¹⁵ (the order parameter is 0.8) suggest that they may be well suited for future investigation of the orientation of melittin at the membrane surface by solid-state NMR experiments.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   NMR property calculations. (A) Backbone 15N chemical shift for the 26 residues of the melittin. The values were calculated on the basis of Eq. 1. (B) Order parameter of the C-H bond for the 26 residues of the melittin.

The helical orientational order of melittin corresponding to PATIR-FTIR amide I experiments was examined. The calculated distribution of the instantaneous values of the angle between the carbonyl C==O bond and the membrane normal is shown in Fig. 7. The contributions of the two helical segments, corresponding to residues 1-13 and 14-25, are shown separately. The distribution of the instantaneous angle is very broad, reflecting the intramolecular fluctuations in the orientation of the carbonyl bonds. The maximum in the distribution is around 110°. Overall, the backbone carbonyl bonds are oriented parallel to the membrane surface, although the hydrophobic helix formed by residues 1-13 makes a larger angle with the membrane surface. Experimental values reported for melittin in phospholipid monolayer or bilayer range from 0.19 to 0.36 (Citra and Axelsen, 1996). The average order parameter S_PATIR calculated from the trajectory is 0.19, in excellent agreement with the experimental result. Nevertheless, the significant spread in the distribution suggests that the data should be interpreted with caution in the case of a flexible membrane-bound polypeptide.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   Distribution of the angle of the C==O bond of the melittin backbone with respect to the bilayer normal. The two segments of melittin (residues 1-13 and 14-25) are plotted separately (the COO of residue 26 was not included).

Solvation of melittin

The atomic density profile of the main components of the system shown in Fig. 4 is in qualitative accord with the amphipathic character of melittin, but it is not sufficient to give a detailed picture of the protein solvation. To better characterize the environment of the protein, the average number of water molecules, acyl chain carbons, and lipid headgroups surrounding each side chain was calculated from the radial distribution functions. The results are shown in Fig. 8. On average, the solvation requirements imposed by the amphipathicity of the protein are satisfied, because polar side chains are exposed to water and the nonpolar side chains are exposed to hydrocarbon chains. A similar solvation pattern was observed in a previous simulation of an amphipathic helix at a membrane surface (Huang and Loew, 1995). However, some of the residues are exposed to a more complex environment. For example, the nonpolar part of the side chain of Lys⁷ is surrounded by hydrocarbon chains, whereas its positively charged nitrogen is in contact with water molecules. Fig. 9 A shows the hydrogen bond complex formed by the interaction of Lys⁷ with interfacial water molecules. The environment of Lys⁷ is in accord with the "Snorkel model" proposed by Segrest et al. (1990) to account for the amphipathic character of lysine side chains. The interactions of Trp¹⁹ with both the hydrocarbon chains and water molecules also reflect the amphipathic character of the indole side chain. Fig. 9 B shows one configuration of Trp¹⁹ at the membrane-solution interface. The observed environment of Trp¹⁹ is in accord with data which suggest that it is limited in its motion by forming hydrogen bonds with water molecules or lipid carbonyl groups (Chattopadhyay and Rukmini, 1993). The amphiphilic nature of the indole side chain is particularly well suited to stabilization of the protein at the membrane-water interface. The charged and polar residues Lys²³-Arg²⁴-Gln²⁵ at the C-terminus contribute to the association of melittin with the membrane by interacting directly with the polar headgroups (see Fig. 8). Electrostatic interactions are important for the stability of the protein at the membrane interface: experimental studies showed that charged lipids prevent cell lysis and that melittin has more affinity for charged lipids than zwitterionic lipids (Beschiaschvili and Seelig, 1990; Monette and Lafleur, 1995). The present simulation indicates that electrostatic interactions with the polar headgroups are important, even if the bilayer is not constituted of charged lipids.

View larger version (41K): [in this window] [in a new window]   FIGURE 8   Contribution from the main components of the membrane system (water oxygens, acyl chain carbons, phosphate and choline polar headgroups) to the solvation of the side chains of melittin. The average solvation number was calculated by counting the number of nearest neighbors within a distance of 4.5 Å around each side chain. The number of neighbors was averaged by normalizing with respect to the total number of atoms in the side chains.

View larger version (43K): [in this window] [in a new window]   FIGURE 9   Hydrogen bonds formed by the side chain of (A) Lys7 and (B) Trp19 with water molecules during the simulation.

Lipid conformations and dynamics

The association of melittin with membranes is thought to depend partly on its affinity for the lipid hydrocarbon chains (Terwilliger et al., 1982). Fig. 10 A shows a top view of the protein-membrane system (water molecules are not shown for the sake of clarity). The hydrophobic segment near the N-terminus is embedded deep in the membrane hydrocarbon. Direct contacts between hydrophobic residues and the acyl chains are observed (e.g., Leu¹³, Leu¹⁶, and Ile¹⁷). The amphiphilic moiety near the C-terminus occupies a large area of the membrane surface. The lipids adopt particular conformations to avoid leaving a large cavity under the helix. The marked asymmetry of the system is reflected in the carbon-deuterium order parameters of the acyl chains. For the upper leaflet, the average order parameters are 0.119 and 0.179 for the Sn-1 and Sn-2 chains, respectively. For the lower leaflet, the average order parameters are 0.151 and 0.224 for the Sn-1 and Sn-2 chains, respectively. The order of the two layers differs by ~20%, the lipids in the same layer as the melittin being less ordered. Although it is not possible to distinguish the upper and lower leaflets experimentally, the overall decrease in lipid order indicated by the calculations is in qualitative agreement with the values observed by Dufourc et al. (1986) in the presence of melittin in DPPC membranes.

View larger version (122K): [in this window] [in a new window]   FIGURE 10   Packing of hydrocarbon chains around the melittin. (A) Top view without water molecules. (B) Side view of a few lipids surrounding melittin.

The disorder in the lipids of the top layer induced by the presence of the amphipathic helix is illustrated in Fig. 10 B. The organization of the hydrocarbon chains around the protein shown in the figure is strikingly reminiscent of the model proposed by Terwilliger et al. (1982). The disorder in the top layer leads to the creation of vacant space in the middle of the bilayer that the chains in the lower layer fill by adopting more extended configurations. Nevertheless, no empty space is left under the amphipathic helix, as indicated by the distribution of free volume shown in Fig. 11. The largest density of free volume accessible to a particle the size of a water molecule is in the center of the membrane. This is in accord with the results obtained for pure membrane simulations (Marrink and Berendsen, 1994). The presence of melittin in the upper layer disturbs the bilayer by breaking the original symmetry of the membrane, and the maxima are shifted. The distribution of free empty volume shows that the membrane is well packed in the polar head region (around ±17 Å), because there is less vacuum space. In contrast, the free volume fraction in the hydrocarbon region is higher because it is less densely packed. Nevertheless, the fraction of free volume is very similar to that observed in pure membranes (Marrink and Berendsen, 1994). Similarly, no significant number of defects appeared in the packing of the hydrocarbon chains in a previous simulation of an amphipathic helix at the surface of a DOPC bilayer (Huang and Loew, 1995). These results suggest that phospholipid bilayers possess a considerable plasticity, allowing them to adapt to the perturbation of an amphipathic helix. It is observed experimentally that bilayers formed by lipids with short acyl chains are more susceptible to the perturbation of melittin (Bradrick et al., 1995). This suggests that the short lipid chains may not have sufficient plasticity to adapt to the large perturbation due to the presence of melittin.

View larger version (14K): [in this window] [in a new window]   FIGURE 11   Average distribution of cavities large enough to accommodate a spherical particle the size of a water molecule along the z axis. The results with all atoms present in the system () and without the water molecules (·····) are plotted. - - -, The fraction of empty free volume with all atoms present in the system.

The dynamics of the hydrocarbon chains is affected by the presence of melittin. Fig. 12 illustrates lipid conformations at different times along the trajectory. The superposition of configurations shows that the movements of the acyl chains close to the protein have a smaller amplitude than others on the opposite layer or farther from the protein. Bradrick et al. (1995) observed that the presence of melittin causes a decrease in fluidity of the surrounding lipids over a large distance (~50 Å) at a temperature below the melting transition temperature of the lipids; the fluidity of the membrane did not seem to be affected at higher temperatures. Even though the present simulation was performed at a temperature slightly above the melting transition, we note that melittin has an influence on the dynamics of the lipids. No significant rotational or lateral translational motions of the lipids is observed, although there are important displacements along the axis normal to the bilayer. As shown in Fig. 4, the distribution of the polar heads extends over a region from z = 13 Å to z = 21 Å, i.e., a distance 4 Å from their initial position (at z = 17 Å). The spread was larger in the bottom layer, from z = 12 Å to z = 22 Å. Different studies have shown that the influence of melittin on the polar headgroups is lipid dependent (Dempsey and Watts, 1987; Beschiaschvili and Seelig, 1990). It has been suggested that the presence of melittin in a PC bilayer alters the average orientation of the headgroup dipoles by shifting the N+ end of the PC headgroups toward the bulk water away from the membrane surface (Kuchinka and Seelig, 1989). In the present simulation, the PC headgroups remained roughly parallel to the membrane surface, and their orientation was not perturbed. Interestingly, it is observed that the dipoles of the PC headgroups of the lipids surrounding the amphiphilic helix are mostly oriented parallel to the protein surface, whereas no preferential orientation is observed on the opposite layer.

View larger version (45K): [in this window] [in a new window]   FIGURE 12   Superposition of five configurations of lipids surrounding melittin. The configurations taken from the trajectory are separated by 100 ps.

Continuum mean-field model

The molecular dynamics simulation of the atomic model provides a detailed view of the melittin-membrane system. Nevertheless, no large movements of the peptide are observed during the trajectory because of the relatively short time scale that can be simulated with current computer resources (i.e., a few nanoseconds). In addition, information concerning the energetic and thermodynamic factors involved in the association of melittin with a membrane surface cannot be directly extracted from the trajectory. For instance, the current dynamics cannot be used to assess the global stability of the starting configuration: an amphipathic -helix roughly parallel to the membrane interface with the unprotonated N-terminus buried in the hydrocarbon core. An important question is whether the starting configuration is consistent with the electrostatic solvation energy of the backbone (the so-called helix dipole).

To better characterize the melittin-membrane system, an implicit mean-field model was used in which the membrane and the surrounding solvent were represented as structureless continuum dielectric material. The model, which is described by Eqs. 4, 6, and 7, is based on numerical solution of the Poisson equation and is essentially the same as that used previously by Ben-Tal et al. (1996). Although such a mean-field model provides a simplified view of the membrane environment, it is nevertheless useful for characterizing the free energy of melittin at the interface between two media of different polarities.

The solvation free energy surface G_tot was explored along two different rigid-body movements of the helix: a translation of melittin along the z axis, and a rocking motion around the helical axis. The calculations were performed using a single configuration of melittin taken from the trajectory at t = 600 ps. From Fig. 13, it is observed that a free energy minimum exists at the interfacial region, at which point the amphipathic helix is partially inserted into the nonpolar region (the interface is located at 12.5 Å in the model). The results of the mean-field calculations differ slightly with the average position and orientation of the helix observed from the molecular dynamics simulation. The average position of the center of mass during the trajectory was around z = 12 Å (the center of mass shifted from 13 to 11 Å during the trajectory). The calculations suggest that melittin is stabilized by about 18 kcal/mol when its center of mass is located around 12-13 Å along the z axis. The optimal position corresponds to a conformation in which the hydrophobic residues are in the low dielectric region, whereas the hydrophilic residues are in the high dielectric region. This is in accord with general ideas about the association of amphipathic helices at the membrane-water interface (Terwilliger and Eisenberg, 1982b; White and Wiener, 1996).

View larger version (16K): [in this window] [in a new window]   FIGURE 13   Solvation free energy of melittin (A) along the z axis and (B) in rotation around an axis going through the center of mass of the protein and parallel to the x axis.

According to the mean-field model, the free energy well corresponding to the membrane-associated state is very broad, and there is no barrier opposing the association of the helix with the interface. A translation of the center of mass from z = 13 to z = 15 and a global rocking of the helix on the order of ±10° correspond to small variations in the free energy of association relative to the thermal energy k_BT. This suggests that spontaneous thermal fluctuations can affect the position and orientation of melittin at the membrane-solution interface. Hints of such global movements are observed during the trajectory; the center of mass of the helix fluctuates from 13 to 11 Å along the z axis, and the magnitude of the rocking motions is on the order of 5-10°. This suggests that the association of melittin with the membrane surface involves significant fluctuations, in terms of both translation and orientation. Such fluctuations contribute to the association constant of the helix, with the membrane as part of the translational and rotational entropy (Ben-Tal et al., 1996). The mean-field calculations indicate that the association of the amphipathic helix with the membrane results from opposite variations in the electrostatic and nonpolar free energy contributions. The unfavorable electrostatic contribution increases as the helix is inserted deeply into the membrane slab, reflecting the energy cost of transferring the polar backbone from the high-dielectric aqueous phase to the nonpolar hydrocarbon membrane region. Complete insertion of the helix (z = 0) is unfavorable because of the large positive electrostatic contribution. Similar trends have been observed previously in continuum electrostatic calculations by Ben-Tal et al. (1996) in the free energy of association of a polyalanine helix with a membrane.

Movements of water molecules and lytic activity

The lytic property of melittin is usually associated with an increase in membrane permeability to water and other small molecules (Dempsey, 1990). It is thus important to analyze the movements of water molecules in the system. Analysis of the solvation of the peptide as discussed above indicates that local interaction with specific residues is responsible, in part, for the increased penetration of the membrane by water molecules. As shown in Fig. 8, a high penetration of water molecules is generally observed in the neighborhood of all of the hydrophilic residues near the C-terminus. This is in accord with experimental data showing that water molecules more readily penetrate the hydrocarbon core of membranes in the presence of membrane-bound proteins (Ho and Stubbs, 1992) or small peptides (Jacobs and White, 1989). More specifically, residues Ile¹⁷-Ser¹⁸-Trp¹⁹ form a particularly favorable region for the penetration of water molecules into the hydrophobic core of the membrane. Interestingly, Trp¹⁹ plays an important role in the lytic activity of melittin; its omission results in a 124-fold decrease in activity as compared to the wild type (Blondelle and Houghten, 1991). In addition to local interactions, the large-scale reorganization of the membrane induced by the presence of the amphipathic helix (i.e., the reduction of membrane thickness and the increased curvature) appears to be correlated with the penetration of the bilayer by water molecules. As observed in Fig. 3 A, a larger number of water molecules are observed in the region where the bilayer is the most distorted. The membrane reorganization involves the upward movement of some lipids of the lower layer and the solvation of the polar heads of water molecules. From a dynamical point of view, it is observed during the trajectory that water molecules left the bulk region to make short transient excursions near the unprotonated N-terminal group of the peptide in the middle of the membrane. The permeating water molecules entered the bilayer from the region where the membrane deformation is the greatest. Marrink and Berendsen (1994) showed that the rate of translocation of water through the bilayer is limited by the interfacial region in which the glycerol backbone and the acyl chains adjacent to the polar headgroups are closely packed. According to their analysis, the presence of defects and vacuum at the interface contributes to the increase in the diffusion of water molecules across the lipid bilayer. In the present simulation, the membrane curvature induced by the amphipathic helix in both the upper and lower layers results in a disordered and less densely packed interface, thus facilitating the penetration of water molecules.

On the basis of neutron scattering measurements with melittin in which the side chains of Ala⁴ and Ala¹⁵ were specifically deuterated, Bradshaw et al. (1994) proposed that the pH of the bulk solution affects the orientation of the peptide in the membrane through the protonation of its N-terminus. According to these observations, melittin with unprotonated N-terminus binds parallel to the membrane surface, whereas melittin with protonated N-terminus binds in a transbilayer way. Interconversion between the two binding modes appears to be possible under equilibrium conditions. These observations have implications concerning the mechanism of membrane lysis. The current simulation shows that water molecules can penetrate transiently into the membrane interior near the unprotonated N-terminus. This suggests that rare excursions of a hydrated hydronium ion, which would provide a mechanism for the protonation of the N-terminus of melittin, may also be possible. Movements of ions near the membrane surface, although not favorable, are energetically possible. For example, potential of mean force calculations have shown that a hydrated potassium ion (similar in size to hydronium) can penetrate deep into the membrane surface (Gambu and Roux, 1997). The interconversion between the binding modes could be the pathway for a number of microscopic processes that would affect the membrane integrity significantly and could subsequently cause membrane lysis (e.g., a large increase in water permeability, formation of channels by association of transmembrane-bound melittin). This view is consistent with a mechanism proposed by Weaver et al. (1992) in which a partial translocation of melittin bound parallel to the membrane surface would expose the hydrophilic residues to the hydrophobic core of the membrane. The extensive water penetration after the protonation of the N-terminus would ultimately disrupt the membrane. Recent works have also suggested a similar hypothesis (Matsuzaki et al., 1997). According to this proposed mechanism, melittin would first bind parallel to the membrane surface with unprotonated N-terminus. The association step is followed by the protonation of the N-terminus, which leads to a conversion into a transbilayer orientation. To assess this proposed mechanism, it is of interest to better characterize the microscopic processes after the protonation of the N-terminus of melittin.

The influence of the protonation state of the N-terminus of melittin was examined using molecular dynamics simulations based on the current atomic model. At t = 600 ps, the N-terminus of melittin was protonated, and the trajectory was continued for 400 ps (note that an isolated water molecule was transiently making an excursion near the unprotonated N-terminus at t = 600 ps). For comparison, the trajectory of the unprotonated system was also continued. According to the simulations, the behaviors of the two systems are markedly different. Fig. 14 shows the number of water molecules within a radius of 4.5 Å from the N-terminus (protonated or unprotonated) of melittin along the two trajectories. Small peaks are observed around t = 70, t = 580-600, and t = 900 ps, corresponding to rare and transient contacts of the unprotonated N-terminus of melittin with water molecules. After t = 600 ps, the penetration of water molecules increased rapidly in the neighborhood of the protonated N-terminus. In contrast, the number of water molecules does not increase during the corresponding period for the system with an unprotonated N-terminus. At the end of the simulation period, water molecules are observed all across the membrane near the protonated N-terminal group, forming a continuous cluster of water. A configuration of the system with protonated N-terminus is shown in Fig. 3 B. Approximatively 15-25 water molecules were in the hydrophobic core of the bilayer from t = 700 ps to the end of the simulation. This is in accord with the observation of Bradshaw et al. (1994), who estimated that the number of water molecules or hydrated protons in the hydrophobic region of the membrane for each melittin was ~20-30. As a consequence of the penetration of water molecules near the protonated N-terminus, the residues Gly¹ to Val⁵ forming the N-terminal segment lost the -helical structure. The N-terminal group moved toward the opposite side of the bilayer (from z = 2 Å at t = 600 ps to z = 6 Å at t = 1000 ps). The computer experiment shows that the protonation of the N-terminus of melittin can have profound effects on the membrane structure.

View larger version (23K): [in this window] [in a new window]   FIGURE 14   Number of water molecules near the N-terminus of melittin. Values were computed by counting the number of waters at a radius of 4.5 Å from the nitrogen of the N-terminus group.

	CONCLUSION

Top Abstract Introduction Methods Results & Discussion Conclusion References

The interaction of melittin with a fully hydrated DMPC bilayer was examined by molecular dynamics simulations. The initial configuration of the system was constructed from preequilibrated and prehydrated phospholipid molecules. The approach was previously shown to be general and applicable to different proteins of arbitrary shape and size (Woolf and Roux, 1994a, 1996; Roux and Woolf, 1996). The initial configuration of the system was constructed with melittin in an -helical conformation bound parallel to the membrane-solution interface. This corresponds to the wedge model that has been described previously (Dawson et al., 1978; Terwilliger et al., 1982). Because melittin perturbs the bilayer significantly, the current application represents a challenging problem for the methodology. In strong support of the present approach, the initial configuration converged rapidly to a stable simulation, despite the marked asymmetry of the system caused by the presence of melittin.

The simulation provides a detailed view of a membrane-bound amphipathic helix at the atomic level. Observable order parameters have been calculated on the basis of solid-state NMR and PATIR-FTIR spectroscopy techniques. Such analysis may be helpful in the interpretation of experimental assays of melittin and other amphipathic peptides. The order of the lipid acyl chains is smaller in the upper layer, whereas it is larger for those in the lower layer. The perturbation of the bilayer results in a local curvature and a reduction of the thickness of the membrane. Because of the presence of melittin, the hydrophobic core of the membrane is reduced by ~30% from its original thickness near the center of the system. However, the acyl chains of the lipids adopt particular conformations to avoid leaving a large cavity under the amphipathic helix. The organization of the hydrocarbon chains in the neighborhood of melittin is very similar to the picture proposed by Terwilliger et al. (1982) nearly 15 years ago. The analysis of the membrane structure shows that the phospholipid bilayer possesses a remakable plasticity and is able to adapt to the significant perturbation caused by the presence of melittin.

To characterize the thermodynamic and energetics of membrane association, a simplified mean-field free energy model, based on continuum electrostatics, was used. The results of the mean-field model suggest that the free energy well corresponding to the membrane-associated state is very broad and that the association of melittin with a membrane involves significant fluctuations. Although such a mean-field model provides only a simplified view of the membrane environment, it is very useful for characterizing the energetics of melittin at the interface between two media of different polarities. The molecular dynamics trajectory provides more atomic detail, but it is limited by the time scale that can be simulated with current resources. The mean-field model is computationally inexpensive and does not suffer from the same limitations.