INTRODUCTION

TOP ABSTRACT INTRODUCTION THE SIMULATION MODELS THE SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

Membrane peptides (and proteins) mediate many important cell processes, including signal transduction, viral infection, cell aggregation, membrane fusion, and membrane rupture (lysing). The peptide's function is dependent on its position and conformation within the membrane, which is controlled by a large number of system variables. Typical conformations include one or more -helical segments that can span the lipid bilayer or adsorb to the bilayer surface. Although many different membrane/peptide systems have been investigated, both experimentally and theoretically, our understanding of the microscopic relationships underlying the observed behaviors is still far from complete. A better understanding of the way in which system parameters affect peptide adsorption, insertion, conformation, and membrane fusion and lysis, would be invaluable in the optimization of these behaviors for biotechnological and pharmaceutical applications.

Experimental studies provide most of our information about membrane/peptide systems, including peptide orientation and depth of penetration (i.e., does the peptide insert), and changes to the bilayer, such as surface area expansion, membrane leakage, rupture and fusion (for a few examples, see Griffith et al., 1974; McIntosh and Simon, 1993; Ho et al., 1995; Longo et al., 1997, 1998; Channareddy and Janes, 1999; Zhou et al., 2000, Nicol et al., 2000). Although some molecular-level interactions can be inferred from experiment, theoretical work is required in conjunction with these results to develop and test detailed molecular-level theories.

The behavior of a membrane/peptide system is controlled by the properties of the peptide/solvent phase, properties of the membrane, and global properties such as temperature and pH. Several theoretical studies have concentrated on the effects of membrane properties on peptide insertion behavior, representing the peptides as simplified solid inclusions of varying geometries. Although they require a number of approximations and generalizations to be tractable, and give no mechanistic information, these models show that lipid composition (which encompasses membrane thickness, density, spontaneous curvature, and bending stiffness), inclusion/membrane thickness mismatches, and the presence of membrane soluble solutes, as well as peptide properties such as concentration in solution, inclusion conformation (geometry), and inclusion hydrophobicity, can all affect the thermodynamic equilibrium state of a membrane/peptide system (Dan and Safran, 1995, 1998; Dan, 1996; Cantor, 1999; Chou et al., 2001). In contrast, peptide folding and the nature of the hydrophobic and polar interactions of peptides with the head and tail regions of the bilayer are emphasized in other theoretical studies of membrane/peptide systems. These include purely theoretical studies (Engelman et al., 1986), statistical mechanical studies (Milik and Skolnick, 1992; Ben-Shaul et al., 1996), Monte Carlo (MC) simulation studies (Milik and Skolnick, 1993, 1995; Baumgaertner, 1996; Sintes and Baumgaertner, 1998; Efremov et al., 1999a,b;), and molecular dynamics (MD) simulation studies (Berneche et al., 1998; Lin and Baumgaertner, 2000). These researchers propose several different models, each of which predicts peptide chain behavior (insertion, adsorption, orientation, etc.) in agreement with experimental results for the same systems. These studies show the relationship between peptide amino-acid composition and the preferred equilibrium state, as well as providing information about mechanistic pathways.

Computer simulations of biological systems can be divided into two general classes; a full-atom model, containing as much detail as possible, and a coarse-grained approximation, containing only those details relevant to the properties of interest. MD simulations are generally used with a full-atom model (although MC techniques can also be applied). In an MD simulation, the position and momentum of every particle, the forces between them, and all external fields must be known at all times. This information is used to simultaneously solve Newton's equations of motion for every particle in the system, to determine its new position and momentum after a very small time step. This procedure is repeated for several million such time steps, usually spanning just a few nanoseconds in real time. The series of particle configurations generated in this way represents a chain of states through which the system passes as it evolves over time. If a suitable set of interaction potentials is available, MD simulations of this type can be extremely useful in providing detailed information about short timescale properties of the system, even though large-scale movements of the peptide are not observed, and great care must therefore be taken in choosing a starting configuration for the system. To study longer timescale behaviors of biological systems (e.g., protein folding, protein insertion into bilayers), the atomic detail must be simplified, an approach known as coarse graining.

Coarse graining recognizes the fact that many of the details in a full-atom model play little or no part in the behavior of interest, or can be replaced by much simpler terms (although such details may be important for a different property of the system). In this way, a much larger system can be simulated for a considerably longer time. We present one such coarse-grained technique that fits somewhere between analytical theoretical methods, and a full-atom simulation. Unlike these methods, our approach allows us to investigate equilibrium states that are kinetically allowed and thermodynamically stable, and consider longer-timescale mechanistic pathways.

There are also many examples of composite models, which include aspects of both full-atom and coarse-grained techniques. In particular, several models have used coarse-grained (mean-field) bilayers, similar to the one presented here, with a full-atom peptide representation (Biggin and Sansom, 1999).

The conditions that facilitate protein immobilization at the surface, self-assembly into a helix and subsequent insertion into the membrane, are of particular interest in the fields of targeted drug delivery, gene therapy, and disease control. Because these are the particular behaviors we are attempting to understand and predict, it is important for our model peptide to be more sophisticated than the simple hydrophobic block used in many theoretical studies, both in the way that it moves and folds and in the details of its interaction with the lipid bilayer. As a result, following Milik and Skolnick (1993), our model emphasizes the complex hydrophobic, polar, and hydrogen-bonding interactions between the membrane and the different amino-acid residues of the peptide, as well as the flexibility of the peptide chain. To avoid overcomplicating the model, this added complexity is offset by a simplification of some other characteristics of the membrane. As a result, we have neglected many of the physical membrane effects mentioned earlier (e.g., spontaneous curvature), and any steric effects.

We apply the MC simulation technique to a theoretical model of a peptide/lipid bilayer system, to investigate the insertion behavior of two peptides, Magainin2 and M2. Magainin2 is an antimicrobial peptide from the skin of the Xenopus frog, that binds to cell membranes (Zasloff, 1987), but does not lyse red blood cells (Bechinger et al., 1991). It is predominantly helical in detergent micelles, but unfolded in water (Milik and Skolnick, 1993), and solid-state NMR spectra suggest that it adsorbs to the surface of a bilayer but does not insert (Bechinger et al., 1991). In contrast, M2 readily inserts into lipid bilayer membranes, forming a trans-bilayer helix, and lysing red blood cells (Bechinger et al., 1991; Kersch et al., 1989).

Initial attempts were made to simulate the insertion of these peptides using the models and methods of Milik and Skolnick (1993) and Baumgaertner (1996), but without success. As a result, our model includes the "chain of spheres" peptide proposed by Milik and Skolnick (1993) and a variation of their effective medium bilayer approximation, along with a new set of amino-acid transfer energies derived from the work of Roseman (1988a,b) and Jacobs and White (1989). In recognition of the growing evidence that the tail region has some polar content and is not simply a homogeneous nonpolar phase, we can vary the tail polarity in our model from alkane-like to octanol-like via a single parameter, f_q, the polarity factor. In addition, we take into account the diminished internal (intramolecular) hydrogen bonding at each end of the helical peptide (Engelman et al., 1986; Roseman, 1988a), which has been neglected in previous coarse-grained membrane/peptide models despite the great effect it can have on insertion behavior. A key feature of our model is that the functions and parameters used to define the membrane (water content, polarity, and hydrophobicity) and its interaction with the peptide (amino-acid transfer energies) were not chosen simply to reproduce the experimental insertion behavior of Magainin2 and M2. They were instead based on several experimental and theoretical studies of phospholipid membrane structure (Jacobs and White, 1989; Milik and Skolnick, 1993; Griffith et al., 1974: Engelman et al., 1986) and an experimental study of the interaction of a lipid bilayer with a series of tripeptides that solubilized in the head region (Jacobs and White, 1989). By setting up our membrane model in this way, it is not specific to the Magainin2 and M2 studies reported here, but should be generally applicable to a wide range of other peptides, allowing us to reproduce and predict their insertion behaviors as well.

The peptides Magainin2 and M2, are similar in many ways, but interact quite differently with a phospholipid membrane. A sophisticated theoretical model is therefore required to reproduce their respective behaviors. However, such a model can do far more than simply reproduce the experimentally observed behavior of each system. It also allows us to follow positional and energetic trajectories of individual amino-acid residues during dynamical processes (e.g., insertion), and collect positional, energetic, and orientational probability distribution data for equilibrated configurations, including metastable states. Therefore, in addition to presenting results that show that our model accurately reproduces the macroscopic behavior of both peptide/membrane systems, we also characterize the stable equilibrium configurations of both peptides (and a metastable equilibrium configuration of M2), and briefly describe two distinct insertion mechanisms for M2 and a translocation mechanism for Magainin2. We also present results for the interaction of M2 with bilayers of different thicknesses, highlighting the dramatic differences in system behavior that can result from a small change in membrane thickness.

We believe that this model simulates peptide behavior during adsorption at the surface, self-assembly into an -helix, and subsequent insertion (under appropriate conditions) more accurately than previous models. Such a detailed understanding of the membrane/peptide system is the first step toward controlling and improving current processes and materials, and designing totally new ones with specific properties. An equally important aspect of this work is that it greatly enhances our understanding of sequence-structure-function relationships.

	THE SIMULATION MODELS

TOP ABSTRACT INTRODUCTION THE SIMULATION MODELS THE SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

The peptide model

The model peptide is essentially the same as that described by Milik and Skolnick (1993), and Baumgaertner (1996). It consists of a linked chain of hard spheres, each representing a peptide residue containing backbone and side-chain elements (-NH-CHR-COO-). The -carbon (C) is located at the center of the sphere, and spheres are linked by 3.8-Å-long virtual bonds. The virtual bonds ensure that the C atoms are always at their experimentally observed equilibrium separation (the mean distance for proteins in the Brookhaven Protein Database), but have no energetic component, acting only as a constraint on the spatial configuration of the peptide. Van Gunsteren and Karplus (1982) have shown that fixing the virtual bond lengths in this way does not adversely affect the simulation dynamics.

In the model of Gregoret and Cohen (1990), the excluded volume of a residue depends on the size of its side-chain. Represented by a sphere, the ideal radius varies from 2.02 ± 0.25 Å for Gly to 2.65 ± 0.45 Å for Lys, with a second sphere (His, Phe, Tyr) and third sphere (Trp) added for the aromatic residues. Such a model requires the primary sphere to be centered at the center of mass of the residue, which is actually somewhat removed from the backbone and the C atom, where our hard spheres are centered. As a result, we are unable to use the residue radii given by Gregoret and Cohen, because they would, in all cases, overlap with adjacent residues along the 3.8-Å virtual bonds. However, because the detailed packing of a peptide is not the focus of our work, we have used the 3.0-Å diameter (1.5-Å radius) spheres proposed by Baumgartner (1996) to more efficiently study the insertion phenomenon. This choice prohibits unrealistic peptide folding, which could occur with smaller hard spheres or no hard spheres at all, while maintaining the simplicity of the model. Our main concern is that the model is able to reproduce the helical folding of a peptide from a random, extended conformation, under the appropriate conditions. This model has such a capability, due to the internal hydrogen-bonding interactions built into the energy calculations, which we describe in the lipid-bilayer model section.

The total energy of the peptide chain, U, is the sum of six energy terms,

(1)

These can be grouped into two different classes; environment-independent terms and environment-dependent terms. There are three environment-independent terms: the angle energy, U_A; the torsional energy, U_T; and the steric energy, U_S. Although indirectly affected by the local environment of the peptide, the calculation of these terms requires only the peptide chain conformation (i.e., the relative position of each residue with respect to all the other residues). The calculation is thereafter quite straightforward using Eqs. 2. There are also three environment-dependent terms: the total hydrogen bonding energy, U_H; the hydrophobic energy, U_B; and the polar (hydrophilic) energy, U_Q. The calculation of these terms requires the position of each residue with respect to the model lipid and surrounding aqueous phase, as well as the overall conformation of the peptide in the case of U_H. Because the calculation of these terms is dependent on the details of our modified lipid-bilayer model, they will be discussed in the lipid-bilayer model section.

The simplest of the environment independent terms is U_S, which is zero when no residues overlap, and positive infinite when any two or more residues overlap:

(2a)

N is the number of residues in the peptide chain, r_ij is the distance between residues i and j, and  = 3.0 Å is the residue hard-sphere diameter. Although listed as an energetic term, U_S acts more as a spatial constraint on the peptide, ensuring that no configuration contains an overlap. U_A is a function of angle  (Fig. 1 a), the angle between adjacent virtual bonds:

(2b)

where  = 2.0 kcal/mol and ₀ = 89.5°. U_T is a function of angle  (Fig. 1 b), the angle between adjacent planes, each of which is defined by two adjacent virtual bonds,

(2c)

where  = 1.5 kcal/mol and ₀ = 52.1°. The angle and torsional energy functions have shallow minima at bond and torsional angles of 89.5° and 52.1° respectively, which correspond to the lowest energy (least strained) conformation of the peptide.

View larger version (13K): [in this window] [in a new window]   FIGURE 1   (a) The angle between adjacent virtual bonds, and (b) the angle between adjacent planes in the peptide model.

The lipid bilayer model

A lipid bilayer is composed of two opposing lipid monolayers. In an aqueous solution, the hydrophobic acyl chain tails face inward toward the center of the bilayer, and the amphiphilic lipid heads face outward, in contact with the aqueous phase. Our standard bilayer model follows Milik and Skolnick (1993), having a 4.5-Å-wide head region at each side and a 27.0-Å-wide tail region in the center (i.e., 2z₀, where z₀ is the length of the acyl chain tail of a single lipid). When a bilayer of a different thickness is modeled, z₀ changes, but the head region remains the same size. Different model bilayers are therefore composed of lipids with similarly sized head groups, but acyl chain tails of different lengths.

Three functions are used to model the head and tail regions of the bilayer, and the surrounding aqueous phase. These functions represent the fractional water content, w(z), the polarity, p(z), and the hydrophobicity, y(z), and are shown in Fig. 2. Each one is a function of the z coordinate of the system only, where the z axis is defined perpendicular to the bilayer surface. w(z) is defined by

(3)

where w defines the shape of the boundary region (where the water content changes from one to zero), z₀ is the length of the acyl chain tail of a single lipid, and w_shift is the displacement of the boundary center from z = z₀. w(z) = 1 in the aqueous phase, on either side of the bilayer region, but drops through the head region, reaching zero a short way into the tail region. Although w = 2 is used throughout our work (fixing the shape and width of the boundary), w_shift = 1.65 Å was chosen by fitting simulation data to the experimental data of Jacobs and White (1989) for the adsorption of a series of tripeptides to the surface of a bilayer membrane (see Appendix B). The fitted w(z) profile matches the observed distribution of water in a lipid bilayer quite well; the amphiphilic head region supports a significant water content, which decreases toward the center of the bilayer, and the outer edges of the hydrophobic tail region contain a small amount of water in the presence of peptide residues (Jacobs and White, 1989). w(z) is used in conjunction with the peptide conformation to calculate U_H, the total hydrogen-bonding energy.

View larger version (24K): [in this window] [in a new window]   FIGURE 2   The model bilayer functions (symmetrical about z = 0). The vertical dashed lines denote the head region.

Eq. 4 gives the phenomenological helicity factor for a residue pair, V_H(r_n), proposed by Milik and Skolnick (1993). It is a function of the measured separation between a residue and its nth neighbor, r_n.

(4)

where a_n is the optimal separation in an -helix, and  is the decay length of the potential. Only the third and fourth neighbors along the peptide chain are considered, because these are the main contributors to the internal hydrogen bonding of a given residue (Baumgaertner, 1996), and the parameters a_n and  are fixed throughout this work (a_|3| = 5.04 Å, a_|4| = 6.3 Å; Barlow and Thornton, 1988;  = 0.1a_n). V_H(r_n) is almost a square well potential, with a maximum value of one around the optimal separation decaying rapidly to zero outside this range. V_H(r_n) is used to calculate H_int(i), the internal hydrogen-bonding energy of each residue in the peptide chain,

(5a)

where H₀ is the maximum hydrogen-bonding energy per residue, including internal (intramolecular) bonding along the peptide backbone and external (intermolecular) bonding with surrounding water molecules. H₀ = 6.12 kcal/mol, which corresponds to the transfer energy of an unbonded peptide group (N---H, O==C), from CCl₄ to water (Roseman, 1988b). This includes the hydrogen-bonding energy (5.50 kcal/mol), and the polar interaction energy of the peptide group (0.62 kcal/mol). We have included the polar interaction energy with the hydrogen-bonding energy because the extent to which both terms contribute to the transfer energy is dependent on the degree of local helicity of the peptide group. A helically folded peptide effectively screens the polar groups of the backbone from the surrounding environment, so there is no energy change on transferring the backbone groups between the aqueous phase and the bilayer. It is therefore appropriate to include the backbone polar interactions in the hydrogen-bonding calculation, rather than with the side chain polar interactions, which are independent of peptide conformation. Our H₀ value is the same one used by Jacobs and White (1989), and Baumgaertner (1996), but different from that used by Milik and Skolnick (1993). Milik and Skolnick used a smaller value (4.1 kcal/mol), which appears to originate from the work of Klotz and Farnham (1968). However, calorimetric studies later showed the thermodynamic cycle from which this value was derived to be incorrect (Kresheck and Klotz, 1969). As a result, the corrected transfer energy for an unbonded peptide group, calculated by Roseman (1988b), is 6.12 kcal/mol. For most residues within the peptide chain, the minimum value of H_int(i) is H₀, but the last four residues on each end of the chain only have a third and fourth nearest neighbor in one direction along the peptide chain, increasing their minimum H_int(i) values to 4.59 kcal/mol for the third and fourth residues from the end, and 3.06 kcal/mol for the end two residues. The external hydrogen-bonding energy of each residue, H_ext(i), is given by

(5b)

This is the energy resulting from hydrogen bonding between the backbone of a residue and the surrounding water molecules. w(z) accounts for the diminished probability that hydrogen bonds will be formed with water molecules when the surrounding phase has a lower water content. The total hydrogen-bonding energy for a peptide chain, U_H, is given by

(6)

The polar energy of a residue, Q(i), is the interaction energy of the charged and partially charged functional groups of the side chain with the solvent environment,

(7)

where q₀(i) is the polar energy of the residue i side chain in water (Table 1), and p(z) is the polarity function,

(8)

where f_q is a polarity factor used to select the relative polarity of the tail region, p defines the shape of the boundary region (p = w = 2), and p_shift is the displacement of the boundary center from z = z₀ (p_shift = 1.5 Å). The value of p_shift is chosen so that p(z) = 1 throughout the aqueous phase and amphiphilic bilayer head region, where Jacobs and White (1989) suggest that all polar interactions are satisfied. p(z) only begins to drop when more than half of the residue has passed into the tail region, reaching its minimum value several Angstroms further into the region, past the small amount of water that gives the tail region some polarity at its outer edges. The total polar energy of the peptide chain, U_Q, is given by

(9)

The hydrophobicity function, y(z), is given by

(10)

where y defines the shape of the boundary region (y = p = w = 2), y_shift1 and y_shift2 are the displacements of the two boundary centers from z = z₀ (y_shift1 = 1.5 Å, y_shift2 = 3.0 Å), and f_y is the fraction of the total hydrophobic energy that is realized when a residue is adsorbed into the head region (f_y = 0.56). The values of y_shift2, and f_y were chosen to best fit the data of Jacobs and White (1989), whose experimental studies of tripeptide adsorption suggest that ~56% of the available hydrophobic energy is associated with adsorption into the head region, and that adsorbed (but not inserted) residues reside primarily at the headgroup midplane (see Appendix B). y_shift1 was chosen to yield the remaining hydrophobic energy as the residue passes into the tail region. y(z) = 0 throughout the aqueous phase, and y(z) = 1 for most of the tail region.

y(z) is used to calculate B(z_i), the hydrophobic energy of a peptide residue at any point along the z axis:

(11a)

where b₀(i) is the water-to-alkane hydrophobic Gibbs transfer energy of residue i (Table 1), and b_hlx(i) is the side-chain hydrophobic energy change of residue i due to the peptide conformation in the aqueous phase. b₀(i) is calculated from A_a(i), the stochastic accessible surface area of a residue (Table 4 of Jacobs and White, 1989; Rose et al., 1985), and C_s, the solvation coefficient for an alkane solvent (C_s = 22 cal/mol Ų; Richards, 1977) as shown in

(11b)

(Reynolds et al., 1974). b_hlx(i) derives from the loss of accessible surface area when the peptide chain has some degree of local helicity,

(11c)

where b₁(i) is the maximum reduction in hydrophobic energy due to helical folding (Table 4 of Roseman, 1988a) (Table 1), and V_H(r_n) is the helicity factor of Eq. 5a. U_B is the total hydrophobic energy of the peptide chain,

(12)

It should be noted that, in general, entropy terms are omitted from energy summations used in computer simulations because the simulation method itself incorporates entropic effectsless likely configurations of the system will be sampled less often. Only the enthalpy (or internal energy, depending on the ensemble) needs to be considered explicitly. However, the energy change that accompanies the transfer of a hydrophobic group from water to a nonpolar solvent (the hydrophobic effect) is due almost entirely to an entropy change, the enthalpy change being negligible (Jacobs and White, 1989; Jain et al., 1985). Because our model has a simplified peptide representation and no water molecules, the local ordering of water around hydrophobic groups, which causes the hydrophobic effect, cannot occur. The simulation will therefore not implicitly account for the entropy associated with water structuring, which must be included in the hydrophobic energy term (this energy term thus becomes a Gibbs energy term), as in any analytical solution. Here, we have avoided the problem of calculating the transfer enthalpies and entropies separately by using empirical Gibbs transfer energies (via C_s and A_a). Other entropic effects, related to the peptide conformation and localization in the bilayer, do not need to be considered explicitly, because they are accounted for by the simulation procedure, even for our simplified peptide model.

The model parameter set

Because there are almost as many residue/bilayer parameter sets as there are models, each one working to the satisfaction of its creators but not widely used elsewhere, there was no clear choice for our work (for a good review, see Engelman et al., 1986). As a result, we have created a new parameter set from reliable primary data taken from experimental and theoretical sources. This parameter set was generated independent of the simulation results, with no parameter fitting. Because it is not tied to any specific peptide/bilayer system, it can be applied equally well to a wide range of systems.

The parameter set (Table 1) consists of three residue-specific values, the water-to-alkane Gibbs transfer energy (hydrophobic energy), b₀(i); the side-chain polar energy in the aqueous phase, q₀(i); and the hydrophobic energy change due to local helicity, b₁(i). The b₀(i) values are calculated via Eq. 11b, and include the hydrophobic energy of the residue backbone, b₀(bb_i), and residue side-chain, b₀(sc_i). Because all residues have the same backbone structure, and glycine has no side-chain, the backbone hydrophobic energy of every residue is equal to the total hydrophobic energy of glycine, b₀(bb_i) = b₀(glycine) = 1.94 kcal/mol. Therefore, the side-chain hydrophobic energy for each residue can easily be found:

(13)

The side-chain hydropathies of Table 3 of Roseman (1988a), corrected for self-solvation, are sums of the polar and hydrophobic energies for water-to-alkane transfer of the residue side-chains at pH 7.0, q₀(sc_i) + b₀(sc_i). Subtraction of b₀(sc_i), calculated in Eq. 13, yields q₀(sc_i). Because the polar interaction energy of the residue backbone, q₀(bb_i) (i.e., the peptide group described earlier), is small and included in the hydrogen bonding energy, q₀(i) is just the side-chain polar energy, q₀(sc_i), which must either be zero for nonpolar groups, or negative for polar groups. A positive value would erroneously infer that a polar group prefers a nonpolar environment to water. However, the method outlined above produces a slightly positive polar energy (ranging from 0.01 to 0.98 kcal/mol) for six residues, and, so, a small correction is applied to b₀(sc_i) for these residues, which reduces q₀(i) to zero. These corrections to b₀(i) are noted in Table 1. The q₀(i) values for some of the polar residues may not appear to be negative enough at first glance. Indeed, simple theoretical estimates and partitioning studies conducted using individual amino acids both suggest that the polar transfer energies should be more positive (q₀(i) more negative) than those presented here. However, Roseman (1988a) shows that such estimates are invalid for amino acids within a peptide chain, and that a variety of effects, including self-solvation, reduce the effective charges on the side chains. As a result, the insertion of polar residues within a peptide chain is actually considerably easier than early studies suggested. Our q₀(i) transfer energies take these charge-reducing effects into account, and are therefore the appropriate values to use for the transfer of residues within a peptide chain. Finally, the b₁(i) values in Table 1 come directly from Table 4 of Roseman (1988a).

Although the q₀(i) parameters given in Table 1 are independent of the bilayer structure, the polarity function, p(z), is dependent on the polarity of the bilayer tail region. The tail region is not simply a homogeneous solvent, but rather a complex tangle of acyl chains, cavities, and a small number of water molecules, giving it characteristics of both octanol and hexadecane-type solvents (Roseman, 1988a). The electron spin resonance studies of Griffith et al. (1974) even suggest a gradient of polarity from the center of the tail region, which is hexane-like, to the edge of the tail region, which is more octanol-like (the arguments made by a number of investigators in favor of an octanol-like or alkane-like tail region are reviewed by Stein (1986)). Our polarity function, Eq. 8, therefore includes a polarity factor, f_q, which allows us to vary the polarity of the bilayer tail region to model all types of bilayer tail polarity. The polarity factor modifies the polar energy of each residue, Q(i), as it enters the tail region, via p(z) in Eq. 7. If the bilayer tail region is assumed to be alkane-like, f_q = 1.00, and if it is octanol-like, f_q = 0.65 (the water-to-octanol transfer energies for polar residues given by Jacobs and White (1989) average around 65% of the corresponding water-to-alkane values). We can also investigate bilayers with intermediate polarity tail regions, for which 0.65  f_q  1.00.

The sum totals of the environment-dependent residue energies, H(i), Q(i), and B(i), at different local helicities, are shown in Fig. 3, a and b, for alanine and aspartine, respectively. Alanine is a nonpolar residue, and aspartine is a strongly polar residue (an alkane-like tail region is assumed in each figure). In each figure, the top curve represents zero local helicity (V_H = 0), whereas the bottom curve represents the completely helical case (V_H = 1). For both residues, a small energy well close to the middle of the head region enables all residues to adsorb to the bilayer when the local helicity is zero. However, a large energy barrier makes insertion into the tail region unlikely for locally nonhelical residues. This barrier is primarily due to the loss of external hydrogen bonds, with an additional loss of polar interactions increasing the barrier height in the case of aspartine. In contrast, there is a dramatic change in the total energy curve when the residues have complete local helicity. In this case, the hydrogen bonding is entirely internal, and therefore indifferent to the water content. The alanine curve now has no barrier to insertion into the tail region, and the aspartine barrier is much smaller. The arrows in the figures illustrate the energetic driving force for helical self-assembly in the head region. After adsorption into the head region, a large energy barrier prevents the further progress of a nonhelical peptide into the tail region, but an increase in local helicity allows each residue to attain a lower energy state. Random conformational fluctuations in the peptide chain that result in greater helicity will therefore be favored. Given a high enough degree of local helicity, the nonpolar residues can spontaneously insert into the tail region, though the strongly polar residues may still be prevented from doing so by the residual energy barrier. So, although all residues favor adsorption and helical self-assembly within the head region of the model bilayer, the balance between polar and nonpolar residues (i.e., their energies and also their position within the peptide chain) will dictate which peptides insert into the tail region, and which do not.

View larger version (26K): [in this window] [in a new window]   FIGURE 3   Total energy functions (symmetrical about z = 0) for a range of local helicities, for (a) alanine and (b) aspartine. Local helicities shown are (from the top), 0%, 25%, 50%, 75%, and 100%. Vertical dashed lines denote the head region.

Magainin2 and M2 are 23-residue amphiphilic membrane peptides that interact quite differently with lipid bilayers (Bechinger et al., 1991). Magainin2 has the sequence
GIGKFLHSAKKFGKAFVGEIMNS, and M2 has the sequence
EKMSTAISVLLAQAVFLLLTSQR. Residues with no side chain (G) or a nonpolar side chain (q₀ = 0 kcal/mol) are shown in bold typeface, and those with a polar side chain (q₀ < 0 kcal/mol) are shown in regular typeface. The N-terminus is to the left of the sequence in each case. All the parameters used to model the different residues are given in Table 1.

	THE SIMULATION METHOD

TOP ABSTRACT INTRODUCTION THE SIMULATION MODELS THE SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

The canonical MC simulation method is used throughout this work (Allen and Tyldesley, 1989), the canonical ensemble being the appropriate choice for our system, in which the volume, temperature, and number of particles remain constant. In a closed system such as this, internal energy and Helmholtz free energy are considered instead of the more familiar enthalpy and Gibbs free energy terms of open systems. However, because the volume change accompanying real peptide adsorption and insertion into a membrane is negligible, the closed and open system terms are equal and interchangeable in this case. Periodic boundary conditions are applied in all directions of the three-dimensional simulation box, so that a particle exiting the box in any direction will reappear at the opposite side. Given a sufficiently large box (at least twice the length of the longest interaction potential between particles), periodic boundaries adequately represent infinite space for a model system. As usual, the minimum image convention for particle interactions is used in conjunction with the periodic boundaries. In this convention, interactions are calculated from the minimum separation between two particles, which may be measured across the simulation box, or across a periodic boundary. Appendix A describes how the chain of system configurations integral to the MC method is generated in our simulations.

MC simulations are most commonly used when the number of particles within the system fluctuates (e.g., in the grand canonical or Gibbs' ensembles). In such cases, the chain of system configurations exists "out of time"a configuration is not necessarily related in time to the previous or next configuration. This is largely because particles may be created within, or deleted from, the system (more accurately, these particles are transferred between the simulated portion of the real system, and somewhere outside the simulated region). In reality, it would take an unknown amount of time for such a configurational change to occur, if such a change was even possible (energetic or physical barriers may rule out altogether the large-scale movement associated with a creation or deletion). So, although the chain of configurations correctly samples phase space giving good ensemble averages, it yields no time-dependent information. In confined systems such as ours, where the particles move around but are never created or deleted from the system, the situation is slightly different. The exact time between successive configurations is still unknown for MC simulations in the canonical ensemble, preventing the calculation of time-dependent properties such as the diffusion coefficient, but successive configurations are much more closely related. If the configurational perturbations are sufficiently small, then it is unlikely that energetic or physical barriers will be bypassed between configurations, and the chain of states can be said to approximate a time-independent relaxation of the system. This means that, as well as enabling us to measure equilibrium properties of the system, we can also examine the relaxation pathway, i.e., the mechanism by which the system changes from its initial state to its equilibrated state.

Unless otherwise indicated, the system temperature was 305 K for all simulations, which is above the thermal phase transition of the lipid membrane. Each simulation run consisted of 50 million MC steps, and took ~4 h to run on a Microway Screamer-LX SuperCache-8 workstation with a 667 MHz DEC Alpha 21164 CPU (Microway, Kingston, MA). Because each MC step consists of three configurational changes, 150 million configurations were sampled during each simulation. The initial configuration was either an aqueous phase random coil, or a trans-bilayer helix. The random coil configuration was generated by a 50,000-step MC simulation of a peptide in the aqueous phase, using a lattice-type peptide structure as the starting configuration. The internal hydrogen bonding energy of the peptide was monitored to ensure that all memory of the lattice-type structure was lost by the end of the simulation. After ~10,000 steps, the internal hydrogen-bonding energy had decreased from an initial value of zero to its final equilibrium value (due to thermal fluctuations, this is actually an equilibrium range). A further 40,000 steps, in which the internal hydrogen bonding energy remained within the equilibrium range, ensured that the peptide was fully equilibrated to an aqueous random coil conformation. For its use in bilayer simulations, this random peptide conformation is shifted laterally, so that it is located entirely within the aqueous phase of the system, and has no initial interaction with the lipid bilayer. The other initial peptide conformation used here is a trans-bilayer helix, which was generated from the random conformation via a 50-million-step simulation using conditions favorable to insertion.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THE SIMULATION MODELS THE SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS APPENDIX A APPENDIX B REFERENCES

The qualitative insertion behavior of our model peptides is summarized in Table 2. We investigated a range of membrane tail polarities, indicated by the polarity factor, f_q, from octanol-like ( f_q = 0.65) to alkane-like ( f_q = 1.00). Aqueous phase and helical trans-bilayer initial configurations were used in each case. The results represent the percentage of runs for which a trans-bilayer helix was the final system configuration, with the number of runs indicated in parentheses.

Starting from a random, aqueous phase configuration, Magainin2 very rapidly adsorbs to the bilayer surface, and assembles into a loose helical structure within the head region. It can move laterally within this region and can translocate through the membrane at f_q = 0.65, but never irreversibly forms a trans-bilayer helix. Starting from a trans-bilayer helix, Magainin2 remains in that conformation in one out of five simulations at f_q = 0.65, and in no simulations at higher f_q. These results show that the equilibrium configuration for Magainin2 in our model bilayer (at all polarities) is a loose helix, adsorbed in the head region. The results from random initial configurations alone do not discount the possibility that Magainin2 could eventually insert, given a long enough simulation run. However, the instability of the trans-bilayer conformation suggests that the barrier preventing Magainin2 insertion is thermodynamic, rather than kinetic. Although the trans-bilayer conformation has a lower internal energy, U, than the adsorbed conformation, it has a higher Gibbs energy at 305 K (the simulation implicitly includes the system entropy and thus generates an equilibrium configuration in which the Gibbs energy, rather than the internal energy, U, is minimized). Therefore, even when it is initially a trans-bilayer helix, Magainin2 will not remain in this thermodynamically unstable state for long. The only exception was at f_q = 0.65, for which the Gibbs energy difference between the trans-bilayer and adsorbed conformations is smallest. Although still unlikely, it is possible for the metastable trans-bilayer conformation to survive for the duration of a 50-million-step simulation in a bilayer of this type. At very low temperatures, the entropy difference is negligible, and the internal energy, U, becomes the deciding thermodynamic factor. Under these conditions, Magainin2 should therefore be more stable in the trans-bilayer conformation than in the head region. This hypothesis was confirmed by a simulation at 50 K and f_q = 1.00, in which trans-bilayer Magainin2 remained in the trans-bilayer conformation throughout a 50-million-step run. It should be noted however, that our model bilayer bears no relation to a real bilayer at 50 K, and the only purpose of this low-temperature simulation was to confirm the entropic contribution to the equilibrium state of Magainin2 at 305 K.

The behavior of M2 is quite different. This peptide also rapidly adsorbs to the lipid bilayer from its initially random, aqueous phase structure, forming a loose helix in the head region. However, for all but the nonpolar, alkane-like bilayer ( f_q = 1.00), M2 inserts and forms a trans-bilayer helix in at least one simulation run. When the peptide inserts, the trans-bilayer form remains stable for the remainder of the simulation. In addition, when the initial configuration is a trans-bilayer helix, the peptide remains in this configuration throughout the simulation. For M2, the helical trans-bilayer form is the thermodynamically stable structure, having the lowest Gibbs energy and the lowest internal energy, U. Only kinetic factors prevent M2 from inserting from the aqueous phase in every simulation. For each bilayer with tail polarity, f_q  0.90, insertion occurs in at least one simulation, proving that the kinetic barrier to insertion is not insurmountable, given enough time. Only the bilayer with alkane-like tail polarity ( f_q = 1.00) resists peptide insertion in every simulation. The kinetic barrier to insertion may be insurmountable for this bilayer at 305 K, or insertion may simply require longer simulations. The reduction in insertion percentage as the tail polarity decreases ( f_q increases) suggests that the size of the kinetic barrier is inversely proportional to tail polarity. This is to be expected because the energy of each inserted polar residue is also inversely proportional to tail polarity. A decrease in tail polarity, therefore, results in higher energy intermediate species and a reduction in the percentage of simulations in which insertion is observed.

We have successfully reproduced the experimental behavior of Magainin2 and M2 in lipid bilayers for a range of tail polarities (the insertion of M2 is not directly observed for f_q = 1.00, but a trans-bilayer helix is shown to be the thermodynamically preferred form, and the possibility of insertion cannot be discounted for a very long simulation run). The success of our model, regardless of tail polarity (within the given range), prevents us from recommending a single f_q value for all bilayer simulations. However, f_q = 0.85 comes closest to the suggestions of Roseman (1988a) and Griffith et al. (1974), that the tail region polarity is somewhere between that of octanol and hexadecane. The following discussion of equilibrium properties will therefore concentrate on peptide behavior in this bilayer model.

An advantage of simulation over experiment is the relative ease with which the microscopic details of a system can be observed. A model that accurately reproduces the macroscopic behavior of a system can also provide a wealth of information at the molecular level. We are therefore able to examine the equilibrium structures of adsorbed Magainin2, adsorbed M2 (metastable state), and inserted M2 in great detail. The results in Table 2 consider the final peptide configuration of each simulation, which tells us whether the peptide inserted or merely adsorbed to the surface. In contrast, equilibrium properties require the sampling of configurational data throughout the simulation. Equilibrium properties are often presented as average values, with no information given about the data spread, but this can be extremely misleading for non-Gaussian data (i.e., data points are not symmetrically distributed about the mean). Here, we consider the probability distribution of chosen properties instead. For the three peptide configurations, Fig. 4 shows probability distribution curves for peptide length, and Fig. 5 shows probability distribution curves for peptide tilt. Length is defined as the separation of end residues, and tilt is defined as the angle , between a line connecting the end residues, and the z axis. In all cases, 0°    180°, and when the N-terminus inserts before the C-terminus,  > 90°. The difference between inserted and adsorbed species is very clear, even though the probability distribution curves are not identical for the two adsorbed species. The average length of inserted M2 is 33 Å, and the fluctuation in length is minimal. This corresponds to a perfect -helix, which is fairly rigid, and has average distance between residues of 1.5 Å (Baumgaertner, 1996). In contrast, the adsorbed conformations have a smaller average length (16 and 23 Å), and a much wider distribution, indicating a bent equilibrium state and greater flexibility. The tilt angles of adsorbed and inserted species are similarly distinct. The adsorbed peptides lie in the bilayer surface, perpendicular to the z axis, indicated by a peak in the probability distribution function at a tilt of 90°. The slight shoulders in these curves indicate very closely related equilibrium variants in which the N-terminus of Magainin2, and the C-terminus of M2 dip slightly into the bilayer. In contrast, the probability distribution function for inserted M2 peaks at 165°, which is almost parallel to the z axis and indicative of a trans-bilayer configuration. The peak width indicates fluctuations in the tilt angle for all species, illustrating the dynamic nature of the equilibrium configurations.

View larger version (26K): [in this window] [in a new window]   FIGURE 4   Peptide length distribution functions for Magainin2, adsorbed M2, and inserted M2.

View larger version (30K): [in this window] [in a new window]   FIGURE 5   Peptide tilt angle distribution functions for Magainin2, adsorbed M2, and inserted M2.

Although inserted M2 was initially expected to lie along the z axis, rather than 15° short of parallel, a series of simulations using different bilayer thicknesses provided an explanation. Figure 6 shows the peak values from the length and tilt probability distribution functions for bilayers of different thickness. The standard thickness used here and elsewhere is 36 Å (Milik and Skolnick, 1993, 1995; Baumgaertner, 1996). Although the peptide length remains constant, the tilt angle is a function of the bilayer thickness. Because the most stable peptide configuration is a perfect helix with its terminal residues at the center of each head region, it assumes this configuration whenever possible. In thinner bilayers, the 33-Å-long helix tends to tilt away from the z axis to accommodate the preferred configuration, rather than compressing or bending to remain parallel to the axis. The thinner the bilayer, the greater the tilt. For thicker bilayers, there is a maximum thickness for which inserted M2 is more stable than adsorbed M2. This occurs at 39 Å for our model. At 40 Å and beyond, the peptide can no longer span the bilayer while maintaining its preferred helical structure, and rapidly assumes the adsorbed configuration, eschewing any intermediate, partially inserted form. This configurational preference is indicated by the sudden change in tilt angle (to 90°) and peptide length (to 17 Å). It should be noted that, although the thickness of our bilayer model is rigidly fixed, the tails of a real bilayer are quite flexible, allowing it to stretch or compress in the z direction, to better accommodate a trans-bilayer peptide. Perturbations of this kind increase the energy of the bilayer, but may decrease the energy of the system by stabilizing the lower energy trans-bilayer configuration of the peptide (Dan and Safran, 1995). However, this type of local thinning or thickening of the bilayer around an inserted peptide is beyond the scope of our current model.

View larger version (27K): [in this window] [in a new window]   FIGURE 6   M2 length and tilt angle as a function of lipid bilayer thickness.

Figures 7 a, 8 a, and 9 a show the equilibrium z coordinates of all residues in adsorbed Magainin2, adsorbed M2, and inserted M2, respectively. Milik and Skolnick (1993), and Baumgaertner (1996) show similar plots, but their z coordinates are simple simulation averages. As noted above, averages omit useful information, and can sometimes be quite misleading. Instead, we generated probability distribution plots for the z coordinate of each residue, and used these data to produce the equilibrium z-coordinate ranges of these figures. As expected, no residue had a fixed z coordinate, but many had extremely broad, and usually asymmetrical probability distributions. Furthermore, some distributions contained more than one peak, indicating multiple locally stable residue locations. The circles and squares therefore represent the z coordinate of the highest peak in the probability distribution (the most likely z coordinate), and the error bars represent the z-coordinate range at half the main peak height. Multiple peaks were only plotted separately when the valley height between them dropped below half the main peak height. In these cases, the alternate, lower peak position is indicated with a dotted line. Residues with polar side chains (q₀ < 0) are shown as circles, and entirely hydrophobic residues are shown as squares.

View larger version (22K): [in this window] [in a new window]   FIGURE 7   Magainin2. (a) residue z-coordinate ranges, and (b) a representative equilibrium configuration. Squares denote nonpolar residues, circles denote polar residues. Dashed lines indicate the head region.

The equilibrium z coordinates of Magainin2 (Fig. 7 a) illustrate several phenomena; first, the four terminal residues at each end prefer to remain in the head region, even though one end is substantially hydrophobic in character. These residues are unable to satisfy all of their hydrogen-bonding requirements internally, and the presence of water molecules in the head region allows some external hydrogen bonding to occur. Second, wherever possible, the polar residues prefer to remain outside, or at the outer edge of the tail region. This is due to the energy increase experienced by a polar residue on entering the tail region, where it can no longer form strong polar interactions. Finally, although some residues are relatively fixed in the z direction, having a small z-coordinate range, others are quite mobile, with an equilibrium range of several Angstroms. In two cases, the residue has dual peaks, showing that it can occupy two essentially separate positions, "snapping" from one position to the other, but spending little time in between. For Magainin2, the distribution of polar residues throughout the peptide serves to anchor its middle section on the edge of the tail region, while the external hydrogen bonding interactions keep its ends in the midst of the head region. Figure 7 b shows a representative Magainin2 equilibrium configuration. Although the peptide spends most of its time in this type of conformation, large fluctuations can occur, and any part of the peptide can insert under the right conditions. Indeed, when the bilayer tail region is most accommodating ( f_q = 0.65), the peptide is able to pass entirely through the bilayer, reforming its equilibrium structure in the opposite head region. This process, known as translocation, occurred in 40% of our simulations involving octanol-like bilayer tails, and is described in more detail below.

The general orientation of Magainin2 and its position in the head region are in agreement with the experimental results of Bechinger et al. (1991), and the simulations of Milik and Skolnick (1993). However, the deeper insertion of the central section of our model peptide, relative to its ends, is the opposite of the Milik and Skolnick predicted conformation. This is due to the incomplete hydrogen bonding of inserted end residues in our model, and, although a minor feature here, such a difference could be fundamental to the behavior of this and other peptides in different circumstances (note: the results of Bechinger et al. do not contain enough detail to confirm either conformation).

Figure 8 a shows the equilibrium z coordinates of the adsorbed M2 conformation. This is a metastable equilibrium state, which shows long-term stability in many of our simulations. It is useful to examine this metastable state to understand the reasons for its stability, and the factors that prevent this conformation from inserting. Although M2 contains a comparable number of polar residues to Magainin2, they are distributed quite differently, being concentrated at the ends of the peptide, rather than evenly spread along the whole chain. The end residues are considerably more stable in the head region of the bilayer than the tail region, due to their polar side chains and reduced internal hydrogen-bonding capabilities. In contrast, the middle section of the peptide is extremely hydrophobic, containing only two residues with polar side chains, and has a lower energy within the tail region of the bilayer. These preferences are optimized in the metastable equilibrium structure, in which the ends remain in the head region, while the middle section extends up to 10 Å into the tail region. To accommodate the surrounding hydrophobic residues, the two centrally located polar residues are also dragged into the tail region at some small energy cost. As a result of this u-shaped conformation, the ends of the peptide are brought closer together, reducing the average length of the peptide, as defined above. Adsorbed M2 is therefore significantly shorter than adsorbed Magainin2, as seen in Fig. 4, though both have very broad length distributions, indicating a high degree of flexibility. Figure 8 b shows a representative equilibrium configuration of adsorbed M2. Like Magainin2, the adsorbed form of M2 shows a good deal of helicity. This metastable conformation for M2 has not been described elsewhere, although its conformational type, dubbed a "helical hairpin" has been characterized by Engelman and Steitz (1981), and modeled by Baumgaertner (1996).

View larger version (22K): [in this window] [in a new window]   FIGURE 8   Adsorbed M2. (a) residue z-coordinate ranges, and (b) a representative equilibrium configuration. Squares denote nonpolar residues, circles denote polar residues. Dashed lines indicate the head region.

Figure 9 a shows the equilibrium z coordinates of M2 after insertion. Inserted M2 forms a fairly rigid