Biophysical Journal Volume 85 December 2003 3431–3444 3431 Interactions of Hydrophobic Peptides with Lipid Bilayers: Monte Carlo Simulations with M2d Amit Kessel,* Dalit Shental-Bechor,* Turkan Haliloglu,y and Nir Ben-Tal* *Department of Biochemistry, George S. Wise Faculty of Life Sciences, Tel Aviv University, Ramat Aviv, Israel; and y Polymer Research Center and Chemical Engineering Department, Bogazici University, Bebek-Istanbul, Turkey ABSTRACT We introduce here a novel Monte Carlo simulation method for studying the interactions of hydrophobic peptides with lipid membranes. Each of the peptide’s amino acids is represented as two interaction sites: one corresponding to the backbone a-carbon and the other to the side chain, with the membrane represented as a hydrophobic profile. Peptide conformations and locations in the membrane and changes in the membrane width are sampled using the Metropolis criterion, taking into account the underlying energetics. Using this method we investigate the interactions between the hydrophobic peptide M2d and a model membrane. The simulations show that starting from an extended conformation in the aqueous phase, the peptide first adsorbs onto the membrane surface, while acquiring an ordered helical structure. This is followed by formation of a helical-hairpin and insertion into the membrane. The observed path is in agreement with contemporary understanding of peptide insertion into biological membranes. Two stable orientations of membrane-associated M2d were obtained: transmembrane (TM) and surface, and the value of the water-to-membrane transfer free energy of each of them is in agreement with calculations and measurements on similar cases. M2d is most stable in the TM orientation, where it assumes a helical conformation with a tilt of 148 between the helix principal axis and the membrane normal. The peptide conformation agrees well with the experimental data; average root-mean-square deviations of 2.1 Å compared to nuclear magnetic resonance structures obtained in detergent micelles and supported lipid bilayers. The average orientation of the peptide in the membrane in the most stable configurations reported here, and in particular the value of the tilt angle, are in excellent agreement with the ones calculated using the continuum-solvent model and the ones observed in the nuclear magnetic resonance studies. This suggests that the method may be used to predict the three-dimensional structure of TM peptides. INTRODUCTION Monte Carlo (MC) simulations have been demonstrated in recent studies to be an important tool in the investigation of peptide-membrane interactions (Milik and Skolnick, 1993, 1995; Baumgartner, 1996; Ducarme et al., 1998; Sintes and Baumgartner, 1998; Efremov et al., 1999a,b, 2002a,b; Lins et al., 2001; Maddox and Longo, 2002). Typically, the MC methods are based on reduced representation of the peptidemembrane system, in contrast to the all-atom representation used in molecular dynamics (MD) simulations. The reduced representation enables comprehensive sampling of peptide conformations and locations in the membrane in an accelerated manner. The MC methods differ from each other in various aspects. For example, the peptide is represented in atomic resolution in some of them and in residue resolution in others. Similarly, whereas most of these studies are based on approximating the membrane as a hydrophobic profile, in Baumgartner’s (1996) study the membrane was represented as a matrix of cylinders, corresponding to the lipid chains. A fundamental difference between the methods is that in some of them the peptide was taken as a rigid helix Submitted October 7, 2002, and accepted for publication June 12, 2003. Amit Kessel and Dalit Shental-Bechor contributed equally to this work. Address reprint requests to Nir Ben-Tal, Dept. of Biochemistry, The George S. Wise Faculty of Life Sciences, Tel Aviv University, Ramat Aviv 69978, Israel. Tel.: 972-3-640-6709; Fax: 972-3-640-6834; E-mail: bental@ ashtoret.tau.ac.il. Ó 2003 by the Biophysical Society 0006-3495/03/12/3431/14 $2.00 (e.g., Ducarme et al., 1998), whereas in others it was flexible. The methods also differ from each other in the potential used in the simulations. For example, in the method developed by Efremov et al. (1999a,b), an atomic solvation free-energy term was added to a standard force field. Another free-energy term was recently added to account for effects due to the transmembrane (TM) potential (Efremov et al., 2002a,b). The methods of Milik and Skolnick (1993, 1995) and Maddox and Longo (2002) were based on combining ad hoc potential energy terms that promote a-helix formation, with hydrophobicity scale. Encouragingly, all these methods qualitatively reproduced the relevant experimental data. We present here a new MC method that can also provide quantitative data, such as the free energy of transfer of the peptide from the aqueous phase into the membrane. A potential was constructed to this effect, and is founded on knowledge-based terms, derived from the proteins of known three-dimensional structure, to govern peptide-conformation changes, in combination with a computationally derived hydrophobicity scale. Each of these terms was first validated separately. The large electrostatic free energy associated with the transfer of unsatisfied hydrogen bonds from the aqueous phase into the hydrophobic environment of the membrane induces secondary structure formation in TM proteins (Kessel and Ben-Tal, 2002); TM proteins adopt helix bundle or b-barrel folds (White and Wimley, 1999). The hydrophobicity scale mentioned above was derived by calculating 3432 Kessel et al. the free energy of transfer of the amino acid residues from the aqueous phase into the membrane in the context of a TM polyalanine a-helix. Thus, the model was tailored for short peptides such as M2d, which are unlikely to form b-barrels, and assume helical structures. (It is noteworthy that recently Efremov and co-workers used their MC method to study the surface adsorption of the b-sheet proteins of the cardiotoxin family; see Efremov et al., 2002a.) In the work described in the accompanying article, we used continuum-solvent model calculations to characterize key thermodynamic aspects of the interactions between M2d, a model peptide of 23 residues, corresponding to the second TM segment of the acetylcholine d-subunit, and lipid bilayers. We determined the most probable peptidemembrane configurations, calculated the free energy of peptide-membrane association, and addressed aspects related to peptide structure and membrane physical properties. In the present study, we used off-lattice MC simulations to study the peptide-membrane association process. We focused on the pathway of M2d insertion into the lipid bilayer, and the effect of the bilayer and aqueous solution on peptide conformations. Moreover, by describing the lipid bilayer as a polarity profile, we were able to consider effects of the water-membrane interface on M2d-membrane interactions. FIGURE 1 Schematic representation of the virtual bond model. A segment between backbone units Cai
2 and Cai11 is shown. The side chain attached to the ith a-carbon is marked as Si. fi is the rotational angle of the ith virtual bond (connecting Cai and Cai
1 ). ui is the bond angle between virtual bonds i and i 1 1. usi (between li and lsi , where lsi is the side-chain virtual bond vector pointing from Cai to Si), and fsi (defined by four consecutive atoms Cai
2 , Cai
1 , Cai ; and Si) are the side-chain bonds and torsional angles, respectively. Internal energy The internal energy of any conformation F ¼ {ui, fi, and fsi ) can be decomposed into contributions from short- (ESR) and long-range (ELR) interactions, as EfFg ¼ ESR fFg 1 ELR fFg: The short-range internal energy was calculated as described in Bahar et al. (1997b) using the formulation of METHODS N
1 The free-energy difference between a peptide in the membrane and in the aqueous phase (DGtot) can be broken down into a sum of differences of the following terms: peptide conformation effects (DGcon); solvation free energy (DGsol); peptide immobilization effects (DGimm); lipid perturbation effects (DGlip); and membrane deformation effects (DGdef) (Kessel and Ben-Tal, 2002; White and Wimley, 1999), as DGtot ¼ DGcon 1 DGsol 1 DGimm 1 DGlip 1 DGdef : (1) The free-energy terms in Eq. 1 were calculated using the MC method described below, in which the peptide structure was approximated using a reduced representation and the water-membrane environment was represented as a structureless smooth hydrophobicity profile. Peptide representation Each residue i was represented by two interaction sites: its a-carbon atom ðCai Þ and its side chain interaction center Si (Fig. 1). The latter were selected on the basis of the specific structure and energy characteristics of the amino acids (Bahar and Jernigan, 1996). The peptide backbone was represented by the virtual bond model originally proposed by Flory and collaborators (Flory, 1969). A peptide of n residues has N–1 virtual bonds connecting successive a-carbons. Virtual bonds are highly stiff and were taken here as fixed at their equilibrium values, li, all of which are of length 3.81 6 0.03 Å. Thus, the peptide backbone conformation can be defined by the 2N–5 dimensional vector {u2, u3, . . . , un
1, f3, f4, . . . ,fN
1} corresponding to n–2 virtual bond angles (ui) at the ith carbon, and n
3 dihedral angles (fi) at the ith virtual bond. Similarly, assuming that the distance between Si and Cai ðlsi Þ is fixed at its equilibrium value, and that the angle between li and lsi ðusi Þ is also fixed at its equilibrium value, the conformation of side chain i can be expressed by the torsion angle ðfsi Þ, defined by li
1, li, and lsi : Biophysical Journal 85(6) 3431–3444 (2) N
1
1 ESR ðFÞ ¼ + Eðui Þ 1 + ½ð1=2ÞðEðfi Þ 1 Eðfi ÞÞ i¼2 i¼2
1 N
1
1 DEðfi ; fi Þ 1 + ½DEðui ; fi Þ i¼2 1 1 DEðui ; fi Þ: (3) The first summation in Eq. 3 refers to the distortion of bond angles, and the second is for the bond rotational angles, in which f
and f1 refer to the rotational angles of the virtual bonds preceding and succeeding the ith acarbon and the coupling between the latter angles. The last term accounts for the coupling between the rotational and bond angle distortions. ELR {F} in Eq. 2 was calculated based on Bahar and Jernigan’s (1996) knowledge-based potential using the expression N
5 N N
4 N ELR ðFÞ ¼ + + WBB ðrij Þ 1 + + WBS ðrij Þ i¼1 j¼i15 N
3 i¼1 j¼i14 N 1 + + WSS ðrij Þ; (4) i¼1 j¼i13 where rij is the distance between interaction sites i and j in conformation F. The backbone-backbone (BB), backbone-side chain (BS), and side chainside chain (SS) interactions were taken into account in WBB, WBS, and WSS, respectively. The functional form of WBB, WBS, and WSS has been derived from data in the Protein Data Bank (PDB); each of them has two minima—characteristics of interactions at the first and second shells in dense systems. Although the suitability of WBB, WBS, and WSS for studies of large and tightly packed proteins has been demonstrated (Bahar et al., 1997a; Haliloglu and Bahar, 1999; Kurt and Haliloglu, 1999, 2001), it is not obvious that this holds for peptides as well. Thus, we carried out MC simulations (of the type described Monte Carlo Simulations with M2d 3433 below) of the folding and stability of short polyalanine-like peptides in the aqueous phase (Shental-Bechor, Kirca, Ben-Tal, and Haliloglu, unpublished data). These preliminary studies showed that when all the long-range interactions are included, compact structures such as helical-hairpins are overly stabilized. Our analysis showed that this phenomenon, which is in conflict with the available experimental data, is attributed to a minimum, which corresponds to the second coordination shell in WBB, WBS, and WSS. Thus, we introduced a cutoff distance of 6.8 Å, corresponding to the first coordination shell, into Eq. 4, and included the contributions of the longrange interaction only for interacting sites below this distance. The exception to this rule is the WBB term that accounts for generic interactions between the backbone atoms of two closely spaced amino acids independent of their type. Interaction between residues i and i 1 5, and between residues i and i 1 6, were taken into account in WBB regardless of the distance cutoff. This means of accounting for the peptide chain’s conformational stiffness ensures the incorporation of proteinlike correlations over several consecutive residues (Levy and Karplus, 1979; Hao and Scheraga, 1995; Skolnick and Kolinski, 1999). 58 each. The entropy was estimated on the premise of independent bond rotations using the familiar ‘‘P ln P’’ relation of 20 72 S ¼ + + pi;j lnðpi;j Þ; (6a) j¼3 i¼1 where pij is the probability of virtual bond j to be at interval i. The value of pij is estimated as pij ¼ ni;j =N; (6b) Calculation of DGcon where ni,j is the fraction of conformations in which a virtual bond j is at interval i of the total number of conformations N. Since the termini are in general disordered, we limited the analysis to the helix core, comprised of 18 bond-rotation angles. We estimated the entropy of the peptide in three cases: in water, in TM orientation, and in surface orientation. In simulations around the surface orientation, the peptide is at equilibrium between the water and the membrane. We were interested only in membrane-associated conformations and therefore considered conformations with a negative DGSIL value (see below). This criterion ensures that at least a small portion of the peptide is in the membrane. The free-energy change due to membrane-induced conformational changes in the peptide (in kT units) can be calculated as DGsol, DGimm, and DGlip DGcon DHcon DScon ¼
; kT kT k (5) where DHcon and DScon are the enthalpy and the entropy changes associated with the transition of the peptide from the aqueous phase into the lipid bilayer, k is the Boltzmann constant, and T is the simulation temperature. A value of T ¼ 1.4 was determined by studying the folding and stability of polyalanine and polyalanine-like peptides of different lengths. Using this value, the experimentally determined Zimm-Bragg parameters and percent helicity were accurately reproduced (Shental-Bechor, Kirca, Ben-Tal, and Haliloglu, unpublished data). Typically, upon association with the membrane, the peptide assumes a well-defined secondary structure, e.g., an a-helix, to avoid the free-energy penalty due to the transfer of unsatisfied backbone hydrogen bonds from polar to apolar media (Kessel and Ben-Tal, 2002). This affects both DHcon and DScon. For example, the confinement of the peptide to the a-helical region in the Ramachandran space leads to negative contributions from the DHcon term due to the stable nature of the helical structure, but at the same time it also involves an entropy penalty. This penalty is the difference between Saq, the entropy in the aqueous phase, and Smem, the entropy in the membrane; i.e., DScon ¼ Saq–Smem. Both the DHcon and the DScon contributions were calculated using the Monte Carlo simulations and the reduced model representation of the peptide. In principle, the free energy itself could have been calculated directly using the internal energies of ample conformations. However, MC (and MD) sampling seeks out lower energy regions, and does not adequately sample the high-energy regions of phase space that make important contributions to the free energy. Thus, the direct calculation of the partition function from the generated conformations may lead to poorly converged and inaccurate results. On the other hand, the free-energy difference between the two equilibrium states is a state function and depends only on the end points. Therefore, we attempted to focus on the properties of the two equilibrium states; the peptide when it interacts with the membrane, and the peptide when it is in the aqueous phase. DHcon was approximated as the effective energy difference, DEcon, and was calculated by averaging over the internal energy of ample peptide conformations sampled in the latter two states. The DScon term was estimated from the difference between the ‘‘volumes’’ in conformation space that are accessible to the peptide in the respective states. To this end, the peptide’s conformation space was estimated from the local conformation-space volume accessible to its individual bonds. The rotation space of each virtual bond was divided into 72 discrete intervals of We developed a hydrophobicity scale based on Dgi, the free energies of transfer of the amino acids from the aqueous phase into lipid bilayers (Kessel and Ben-Tal, 2002). Unlike other hydrophobicity scales (e.g., Kyte and Doolittle, 1982), which assume that the free energies of transfer are proportional to some inherent property (of an individual atom, a chemical group, or an amino acid), this scale was derived using the continuum-solvent model (Honig and Nicholls, 1995), and accounts for the amino acids being located at the center of an a-helix. The scale, which includes free-energy contributions from peptide solvation and immobilization, and from lipid perturbation effects (Dg_SIL_i ¼ Dgsol_i 1 Dgimm_i 1 Dglip_i), is presented in Table 1. Due to the excessive free-energy penalty associated with charge transfer from the aqueous phase into the hydrocarbon region of the bilayer (Honig and Hubbell, 1984), the titratable residues were assumed to be neutral (taking into account the freeenergy penalty of neutralizing them in bulk water). Overall, the scale is similar to other hydrophobicity scales in that the hydrophobic and polar/ titratable amino acids are at the two extremes. However, it is less hydrophobic than other scales. This is mainly because, unlike other scales, it includes the free-energy penalty of inserting the helix backbone into the bilayer. Our previous studies have indicated that the inclusion of this penalty is crucial for reproducing the experimental values of the free energy of peptide transfer from the aqueous phase into lipid bilayers. Indeed, using this approach, we successfully reproduced the transfer free energy of polyalanine from the aqueous solution to the lipid bilayer (Ben-Tal et al., 1996). Moreover, preliminary tests have shown that the scale is significantly more potent in detecting TM helices in the sequence of membrane proteins compared to other hydrophobicity scales (Chen et al., 2002; Ben-Ami, Honig, and Ben-Tal, unpublished data). Incorporation of the hydrophobicity scale into the reduced model The free energy (Dg_SIL_i(z)) of transferring the ith residue (at a given conformation) from the aqueous phase to a distance z from the bilayer midplane can be decomposed into contributions from its backbone ðDgbi Þ and side chain ðDgsi Þ: We defined DGSIL¼ DGsol 1 DGimm 1 DGlip as the value of the free energy of transferring the whole peptide from the aqueous phase to a distance z from the bilayer midplane as b b s s DGSIL ¼ +ðpi ðzÞDgi 1 pi ðzÞDgi Þ: (7) i Biophysical Journal 85(6) 3431–3444 3434 Kessel et al. TABLE 1 Amino acid DGi (kcal/mol) I L F V A G C S T M W P Y Q H K N E D R N–H C¼O
2.6
2.6
1.5
1.2
0.2 0.0 10.4 10.8 11.1 11.3 11.3 12.8 14.3 15.4 16.8 17.4 17.7 19.5 111.5 119.8 11.8 12.5 FIGURE 2 The hydrophobicity profile. Hydrophobicity, denoted as p, as a function of distance from the bilayer midplane z. The distance between the bilayer midplane and the profile’s torque point is marked as z0. The sharpness of the curve is determined by h. A value of 0.5 (solid line) was used here. At h ¼ 5 (dashed line) the water-membrane interface virtually does not exist. A hydrophobicity scale representing free energies of transfer of each of the 20 amino acids from water into the center of the hydrocarbon region of a model lipid bilayer (Dgi). The scale was computationally derived, as described in Kessel and Ben-Tal (2002). The amino acid residues are presented using a single letter code. The values include the free-energy penalty due to the transfer of the backbone hydrogen bond from water into the membrane. The last two rows present an extra free-energy penalty associated with the transfer of unsatisfied backbone N–H and C¼O hydrogen bonds from water to the membrane. This penalty was added to the water-to-membrane transfer free energy of amino acids in conformations that are incompatible with hydrogen-bond formation. The prefactors in both terms, representing the hydrophobicity of the environment, are proportional to the distance (z) between the interaction site and the bilayer midplane. A sigmoidal function (La Rocca et al., 1999; Milik and Skolnick, 1993; Baumgartner, 1996; Ducarme et al., 1998; Maddox and Longo, 2002; Biggin et al., 1997) of the type (Fig. 2) b pi ðzÞ ¼ 1=ð1 1 expðhjz
zm jÞÞ (8) and a similar expression for psi were used. In Eq. 8, zm is the width of a lipid monolayer as defined in the subsection below. The value of h (h [ 0) determines the sharpness of the profile, i.e., the characteristic length of the membrane-water interface (Fig. 2); h ¼ 0.5 was used here, whereas h [ 5 corresponds to the steep profile of the slab model used in the continuum calculations of the accompanying article. The free energy of transfer of the peptide backbone from the aqueous phase into the membrane (at a given conformation) was calculated using the set of expressions in Eq. 9, b N
H Dgi ¼ gi b i Dg ¼ f 3 ðg b i Dg ¼ g N
H i C¼O i 1f 3g 1g C¼O i 1f 3g C¼0 ði ¼ 1; 2; 3Þ; Þði ¼ 4; 5; . . . ; n
3Þ; N
H i ði ¼ n
2; n
1; nÞ; (9a) (9b) (9c) where the prefactor f is proportional to backbone deviations from the optimal a-helical conformation as observed in Bahar and Jernigan (1996), Biophysical Journal 85(6) 3431–3444 2
2 f ¼ 1=2f2
exp½ð
1=s Þðfi
f0 Þ  2 1 2
exp½ð
1=s Þðfi
f0 Þ g: (10) That is, f is assigned a value of zero for residues in their ideal a-helical conformations obtained at f0 ¼
1208 (Bahar et al., 1997b); for these residues, the C¼O and N–H backbone groups partially neutralize each other. However, the value of f approaches 1 for residues that deviate significantly from the ideal a-helical conformation, e.g., residues that are in extended conformations. For these residues the free-energy penalty due to the transfer of both the C¼O and N–H backbone groups is taken into account. Obviously, the presence of Eq. 10 in the potential strongly enhances the formation of helical structures upon membrane association. s is the standard deviation in the distribution of angles around the optimal a-helical conformation. We estimated a value of s ¼ 308 based on Fig. 3 in Bahar et al. (1997b). It is noteworthy that the stretches of three residues at the N- and C-terminals are treated differently than the peptide core (Eq. 9). The freeenergy penalty associated with the transfer of the uncompensated hydrogen bonds of the N–H groups of the three residues at the N-terminal is taken into account regardless of the peptide conformation (Eq. 9 a). Likewise, the freeenergy penalty due to the transfer of the uncompensated hydrogen bonds of the C¼O groups of the three residues at the C-terminal is also taken into account, regardless of the peptide conformation (Eq. 9 c). Calculation of DGdef Insertion of a rigid hydrophobic inclusion into a lipid bilayer may result in a deformation of the lipid bilayer to match the width of the hydrocarbon region to the hydrophobic length of the inclusion, following the mattress model (Mouritsen and Bloom, 1984). The deformation involves a freeenergy penalty, DGdef, resulting from the compression or expansion of the lipid chains. DGdef has been calculated for lipid bilayers composed of lipids of various types using different methods and yielding similar values (e.g., Mouritsen and Bloom, 1984; Helfrich and Jakobsson, 1990; Fattal and BenShaul, 1993; Ben-Shaul et al., 1996; Dan and Safran, 1998; Nielsen and et al., 1998; May and Ben-Shaul, 1999). We rely on the calculations of Fattal and Ben-Shaul (1993) that are based on a statistical-thermodynamic molecular model of the lipid chains (14-carbon lipids were used), and fitted a harmonic potential of the form Monte Carlo Simulations with M2d 3435 FIGURE 3 Folding of M2d on the membrane surface, starting from an extended conformation. The figure presents three snapshots taken at 300 (A), 3300 (B), and 22,500 (C) MC cycles from a representative simulation, carried out using the regular MC sampling protocol. The membrane hydrophobicity profile is color-coded so that dark black represents the most highly hydrophobic region of the lipid chains, gray represents the membrane-water interface, and the aqueous phase is white. The peptide’s N-terminus is marked with an arrow. An internal-toexternal motion ratio of 10:1 was used. DGdef ¼ vðzm
z0 Þ 2 (11) to their calculations. zm and z0 are the actual and native widths of a monolayer. We used a value of z0 ¼ 15 Å based on the available experimental data (White and Wimley, 1999; Nagle and Tristram-Nagle, 2000) and the theoretical studies mentioned above. v is a harmonic force constant related to the membrane elasticity. A fit to the calculations of Fattal and Ben-Shaul (1993) gives v ¼ 0.22 kT/Å2. maximum variations (subjected to the adjustment) of the random perturbations of a, b, g, and d. In general, the maximum variations are adjusted for a reasonable acceptance ratio of the attempted moves. In the current implementation, the proportion of internal versus external moves is also important. These parameters were determined by trial and error, such that the system will have enough time for internal relaxation but will not be trapped too often in local energy minima. Monte Carlo protocols Generation of conformations New conformations were generated by simultaneously subjecting the generalized coordinates fi, fsi , and ui to random perturbations of the type Dðfi ; ui ; fsi Þ ¼ dkð2r
1Þ; (12) where r is a random number between 0 and 1, and dk is the maximum variation for the respective coordinates. The maximal variation in fi was taken as 38, and a value of 0.58 was used for ui and fsi : Our experience has been that the simulations are not effected by slight changes in the magnitude of dk (Shental-Bechor, Kirca, Ben-Tal, and Haliloglu, unpublished data). The peptide has a chance of changing its configuration within the membrane mainly by external motions as described below. However, it is noteworthy that a set of randomly chosen conformational changes may eventually lead to slight changes in the orientation of the peptide in the membrane. Generation of configurations The external rigid body rotational and translational motions were carried out to allow the peptide to change its configuration, i.e., location in, and orientation with respect to, the membrane. These motions were employed respectively as anew ¼ aold 1 2ðr
1Þdamax cos bnew ¼ cos bold 1 2ðr
1Þdðcos bmax Þ (13a) (13b) g new ¼ gold 1 2ðr
1Þdg max (13c) rnew ¼ rold 1 2ðr
1Þddmax ; (14) and where a, b, and g are three Euler angles describing the orientation, and r represents the Cartesian coordinates of the peptide’s geometric center. da, db, dg, and dd (¼58, 58, 58, and 0.02 Å, respectively) were chosen to be We followed a standard MC protocol in sampling the conformational space of the peptide in the aqueous phase; the acceptance of each move, which creates a new conformation from the present conformation, was based on the Metropolis criterion and the internal energy difference (DDE, Eq. 2) between the new and old states (Metropolis et al., 1953). M2d is a 23-residue peptide. Thus, each MC cycle was composed of N ¼ 23 random perturbations of randomly chosen generalized fi, fsi ; and ui coordinates. The equivalent protocol for sampling conformational/configurational space of the M2d-membrane system would involve using the same criterion and the free-energy difference DDGtot ¼ DDE 1 DDGsol 1 DDGimm 1 DDGlip 1 DDGdef ¼ DDE 1 DDGSIL 1 DDGdef : (15) However, although both experimental (Bechinger et al., 1991; Opella et al., 1999) and theoretical (Milik and Skolnick, 1993; Law et al., 2000; Maddox and Longo, 2002) studies indicate that M2d inserts into lipid bilayers and resides in a TM orientation, it is very well known that membrane insertion involves crossing a large free-energy barrier due to the electrostatic freeenergy penalty of transfer of at least one of the peptide’s polar termini from the aqueous phase into the hydrocarbon core of the lipid bilayer (Kessel and Ben-Tal, 2002; White and Wimley, 1999). Overcoming this barrier requires a large free-energy fluctuation, and our preliminary simulations indicated that it is very unlikely to occur during the simulation (data not shown). Thus, for the peptide-membrane system, we used a revised MC protocol composed of two stages, for efficient sampling of the free-energy surface. The first stage was designed for an efficient and rapid sampling of various peptide-membrane configurations, while considering the exact peptide conformation as of secondary importance. The goal at this stage was to detect peptide-membrane configurations of low free energy (e.g., when the peptide is in surface and TM orientations). To this end, we replaced DDGtot in Eq. 15 with DDGtot ¼ DDE 1 ln (DDGSIL 1 DDGdef) when applying the Metropolis criterion to the external motions. The use of a logarithmic function for the membrane-related free-energy contributions significantly reduces the barrier of peptide insertion. To further enhance the membraneinsertion probability, we used an internal-to-external-motion ratio of 1:1. Biophysical Journal 85(6) 3431–3444 3436 Obviously, this crude search is likely to provide only approximated peptidemembrane configurations. (It is noteworthy that utilization of the logarithm of the energy, rather than the energy itself in MC simulations, has been proven to be useful in protein structure predictions as well; see Zhang et al., 2002.) This crude search provided two main free-energy minima: one corresponded to surface-adsorbed and the other to TM-inserted M2d. In the following stage, a standard MC sampling with the Metropolis criterion and the free-energy difference of Eq. 15 (without lowering the barrier heights for the external motions) was used in a fine-tuned search in conformation/ configuration space in the vicinity of each of these two minima. Experimental structures used in the study The structure of M2d in the simulations was compared to the following structures of the peptide, which have been previously determined using nuclear magnetic resonance (NMR) techniques: 1. The structure of M2d in dodecylphosphocholine (DPC) micelles (Opella et al., 1999; PDB entry 1A11). 2. The structure of M2d in dimirystoylphosphatidylcholine (DMPC) bilayers (Opella et al., 1999; PDB entry 1CEK). 3. The molecular model of the AChr M2 channel structure (Opella et al., 1999; PDB entry 1EQ8). RESULTS Kessel et al. entropy are different for the two interval sizes, DScon is stable and does not depend on the interval size. Membrane association Folding of M2d on the surface Simulations of M2d, starting from a randomly generated conformation near the membrane surface, were carried out with the standard MC sampling procedure described in Methods. The simulations demonstrated the following path of peptide association with the lipid bilayer (Fig. 3). First, the extended chain adsorbed onto the bilayer surface (Fig. 3 A). As a result of its interaction with the lipid bilayer, the chain gradually assumed a highly ordered helical structure (Fig. 3, B and C). M2d has a hydrophobic core and two polar termini. As a result, it assumed a curved, bridgelike conformation, with both unstructured polar termini anchored in the watermembrane interface, and the helical hydrophobic core slightly immersed in the hydrocarbon region of the membrane (Fig. 3 C). The bananalike helical structure persisted throughout the simulation. A similar adsorption/folding pathway was observed in 13 different runs. M2d in water Insertion of M2d into the membrane Starting from an extended conformation, we simulated M2d in the aqueous phase using the standard Metropolis MC sampling scheme described in Methods. Eight independent runs of 5 million MC cycles each were carried out and the results were averaged over all of them. A very low average helicity of 0.2(60.05)% was obtained, indicating that M2d is a random coil in the aqueous phase. This is in agreement with the MD simulations of Law et al. (2000). The average effective conformational energy in the aqueous phase was hEconiaq 
49 (64.9) kT. The conformational entropy, Scon, of M2d in the aqueous phase was estimated from the conformation-space volume of its individual bonds, as explained in Methods. In short, the conformation space was divided into discrete intervals of 58 each. The analysis showed that each of these intervals was occasionally occupied by the unstructured peptide bonds during the simulations. We calculated the frequency that an interval of the torsion space was occupied by the peptide. The two bonds at each peptide terminus were flexible both in the aqueous phase and in the membrane so we limited the analysis to the 18 central bonds, and obtained a value of (Scon)aq ¼ 64.7 6 2.5 k. (It is noteworthy that the absolute values of (Scon)aq and (Scon)mem—i.e., the conformational entropy of the peptide in the aqueous phase and in the membrane, respectively—reflect the sampling procedure and are therefore meaningless; only the entropy change DScon ¼ (Scon)mem
(Scon)aq is meaningful.) To check the robustness of the entropy estimate to interval size, we divided the space into 12 intervals of 308. Although the absolute values of the M2d insertion into the lipid bilayer involves a high freeenergy barrier, resulting from the transfer of at least one of the highly polar segments at the peptide termini from the aqueous phase into the hydrophobic core of the bilayer. To facilitate the insertion process, we used the revised MC protocol as described in Methods above, which lowers this free-energy barrier, and allows the system to explore unfavorable configurations with high free energy. In addition, we reduced the desolvation free-energy penalty values of residues E1 and K2 by half. This modification is in accordance with the peptide sequence, which suggests that these two residues may be salt-bridged. The salt-bridge is likely to reduce the polarity of both these residues significantly (Honig and Hubbell, 1984). The path of M2d association with the lipid bilayer, obtained using the revised MC protocol, is shown in Fig. 4. The extended peptide (Fig. 4 A) first adsorbed onto the membrane surface, assuming a helical structure. As the central part of the chain diffused into the hydrocarbon region of the bilayer, the peptide acquired a bridgelike form and its conformation became more ordered (Fig. 4, B and C). As the simulation advanced, the peptide core became more immersed in the hydrocarbon region and its curved conformation turned into a helical-hairpin (Fig. 4, D and E). From this point and on, several attempts to cross the barrier between the surface and TM orientations were made. One of these attempts, involving a translocation of the (modified) N terminus across the membrane, was successful (Fig. 4, F and G) and resulted in TM orientation (Fig. 4 H ). Biophysical Journal 85(6) 3431–3444 Monte Carlo Simulations with M2d 3437 FIGURE 4 Insertion of M2d into the membrane. Snapshots from a representative simulation carried out using the revised MC protocol. (A) 0 MC cycles, (B) 240,000 MC cycles, (C) 600,000 MC cycles, (D) 1,100,000 MC cycles, (E) 1,200,000 MC cycles, (F) 1,250,000 MC cycles, (G) 1,280,000 MC cycles, and (H) 1,500,000 MC cycles. The membrane hydrophobicity profile is color-coded as in Fig. 3. The peptide inserts into the membrane with its N-terminus (marked by arrow) first. An internal-to-external motion ratio of 1:1 was used. Fig. 5 displays the change in the peptide’s internal energy with respect to the value in the aqueous phase, the corresponding solvation free-energy change, and the change in zc as a function of the number of MC cycles. Peptide folding on the bilayer surface involved a decrease both in the internal energy and in the solvation free energy. At a certain point, the system gained enough solvation free energy to overcome the barrier due to peptide insertion (marked by asterisks in Fig. 5, B and C), and M2d assumed a TM orientation, with zc fluctuating around the bilayer midplane. The simulations were carried out using the modified MC protocol described above and the large zc fluctuations observed for M2d in the TM orientation are due to the overly smoothed nature of the free-energy surface used. Local search around the low energy configurations The simulations described above show two stable peptidemembrane configurations: surface and TM. These simulations were carried out by the use of the revised MC sampling protocol, which produced a coarse free-energy surface. To derive more accurate free-energy values for the transfer of M2d from the aqueous phase into the TM and surface orientations in the bilayer, we carried out further simulations around these two orientations using the regular MC sampling. Nine different simulations were carried out around the TM orientation and 13 around the surface orientation. Table 2 shows the thermodynamic characteristics of the average surface and TM orientations, obtained from these simulations. As Table 2 demonstrates, both are stable, with the latter the most stable, in agreement with our continuumsolvent studies (see accompanying article), and the studies mentioned above. The thermodynamic properties of Table 2 are analyzed below in the Discussion. We carried out a clustering analysis of the structures both in the TM and surface orientations. The analysis focused on the peptide core (residues 3–21) alone, since the terminal segments were relatively flexible and it is difficult to tell if this is an inherent property of the peptide or an artifact of the model. The analysis showed that most surface conformations have a helical bananalike core, but they fluctuate significantly at the terminal segments. Approximately one-half of the conformations correspond to a single cluster (Fig. 6 A), whereas the rest are distributed over many other clusters, most of which are singletons. The TM conformations correspond to another cluster of regular helical structures (Fig. 6 B). The 10 NMR structures (1A11) used in the continuum-solvent study (see accompanying article), and the structure of the peptide in the bilayer (1CEK) would also belong to the latter cluster. Our analysis also showed that an ideal a-helix would reside in the same Biophysical Journal 85(6) 3431–3444 3438 Kessel et al. TABLE 2 TM orientation DGtot* (kT) DGcony (kT) DEz (kT)
TDScon§ (k) DGSIL{ (kT) DGSIL_bk (kT) DGSIL_s** (kT) DGdefyy (kT) hzizz (Å) abs(zc)§§ (Å) hui{{ (8)
10.2
5.4
35 29.6
7.6 8.35
15.91 2.8 11.53 1.6 13.95 6 6 6 6 6 6 6 6 6 6 6 5.6 5.6 5 2.6 0.1 0.12 0.01 0.06 0.04 0.2 0.3 Surface orientation
7.4
3.3
14.4 11.1
4.3 3.5
7.8 0.2 15.4 24.6 86 613.4 6 11.6 6 9.6 6 6.6 6 2.5 6 1.2 6 3.7 6 0.2 6 0.4 6 3.2 6 3.2 Thermodynamic parameters for the membrane association of M2d in TM and surface orientations. The values represent the average obtained after equilibration. *Total free energy (DGtot). y DGcon ¼ DE
TDScon. z Internal energy, calculated with reference to hEaqi ¼ 49 kT. § The conformational entropy component. { Total external free energy, including peptide solvation and immobilization and lipid perturbation effects: DGSIL ¼ DGsol 1 DGimm 1 DGlip. k Backbone external free energy. **Sidechain external free energy. yy Membrane deformation free energy. zz The width of a monolayer. §§ Absolute value of the z-coordinate of the centroid of the chain; midplane of the membrane is at z ¼ 0. {{ The projection angle of the peptide’s end-to-end distance vector r and the membrane normal z. The values in parentheses depict the standard deviation of different runs for given case. FIGURE 5 (A) The internal energy change (DE) with respect to aqueous phase (Eaq ¼ 49 kT) vs. the number of MC cycles. (B) The total external free-energy change (DGSIL) vs. the number of MC cycles. (C) The change in zc (the distance between the peptide’s centroid and the membrane midplane) as a function of the number of MC cycles. An internal-to-external motion ratio of 1:1 was used. cluster, suggesting that this cluster represents a fluctuating canonical a-helix. A detailed comparison between the average RMSD value of the predicted conformations in the TM and surface orientations and the 10 NMR structures is presented in Table 3. The comparison further emphasizes the similarity between the predicted and experimental structures. The agreement is particularly good for the TM structures, where the average RMSD is typically \2 Å. Solid-state NMR provides information about the orientation of the peptide in oriented lipid bilayer (Opella et al., 1999). With the use of this information we can infer the tilt of the peptide in the lipid bilayer, and the helix rotation about the principal axis. Our simulations show that the average tilt angle about the membrane normal is 148, which is in good agreement with the value of 128 reported by Opella et al. (1999). In the model of the AchR M2 pentameric channel (PDB 1EQ8), the wide mouth of the funnel is on the N-terminal part of the peptide and the pore lining residues are E1, S8, V15 L18, and Q22. To check the position of these residues in the simulations, we randomly sampled Biophysical Journal 85(6) 3431–3444 10 conformations, all of which matched the orientation suggested by Opella et al. (1999). An estimate of DGcon and DGtot The calculation of DScon, the entropic component of DGcon, was carried out using Eq 7. Out of the total of 72 intervals, only a few, corresponding to the rotation angles that are FIGURE 6 Clusters of M2d conformations. Two clusters, corresponding to the peptide in surface, A, and TM, B, orientations, obtained from one simulation. The conformations were clustered according to a similarity measure of RMSD \ 2.5 Å. (A) The largest of all the clusters that were obtained in the surface orientation corresponds to a bananalike helix structure. (B) The single cluster that was observed in the TM orientation corresponds to a canonical a-helix structure. Monte Carlo Simulations with M2d 3439 TABLE 3 NMR structure TM hRMSDi (Å) Surface hRMSDi (Å) 1A11_1 1A11_2 1A11_3 1A11_4 1A11_5 1A11_6 1A11_7 1A11_8 1A11_9 1A11_10 1CEK 1.93 2.30 1.77 2.20 1.80 1.93 2.09 1.84 2.60 1.76 1.3 3.54 3.91 3.85 3.77 3.91 3.54 3.84 3.81 3.25 3.63 Average root-mean-square deviation (hRMSDi) of peptide conformations, in the TM and surface orientations, compared to the NMR structures of Opella et al. (1999). The RMSD of each conformation was defined as 1/19 Si¼3
21 ðRai
Raix Þ2 , where Rai and Raix are the position vectors of Cai in the predicted and NMR structures, respectively, after superimposition. characteristic of helix conformations, were sampled by the membrane-associated peptide bonds in the TM orientation (Strans ¼ 35.1 (60.7) k). In the surface orientations, the peptide is less stable and a larger part of the conformational space is accessible (Ssurf ¼ 53.7 (66.1) k). In contrast, all the 72 intervals were accessible to the peptide bonds in the aqueous phase (see above) with high probability (Saq ¼ 64.7 (62.5) k). The entropy values of the TM and surface conformations are in accordance with the results of the clustering analysis. The TM conformations are more rigid, and we therefore obtained one cluster of conformations, and a low value of the entropy. The surface conformations are more flexible; hence they are distributed over many clusters, and the entropy value is high. The average values of the internal energy change (DE) due to the transfer of M2d into the TM and surface orientations are ;
35 kT and
14 kT, respectively (Table 2). Using Eq. 6 and the estimated DScon, DGcon¼ ;
5.4 kT and ;
3.3 kT is obtained for the TM and surface orientations (Table 2). Using these values and Eq. 1, we obtained DGtot values of
10.2 kT and
7.4 kT for the transfer of M2d from the aqueous phase into the membrane in the TM and surface orientations (Table 2). DISCUSSION A novel MC method was used here to study the interactions of M2d with lipid bilayers. In the following we discuss a few central features of this method and its main limitations in view of the complexity of peptide-membrane systems. We then compare it to other methods and conclude with a discussion of the significance of the results. Features and limitations of the model The MC protocol used here involves a rough exploration of the peptide-membrane conformational/configurational space using an approximated potential energy as described above. Peptide-membrane conformations/configurations, which appear to be stable, are then explored further by carrying out more MC sampling using the ‘‘exact’’ potential of mean force of Eq. 15. In the M2d example presented here, two favorable peptide-membrane configurations were detected: surface (Fig. 4 B) and TM (Fig. 4 H), as anticipated. Less obvious configurations may be discovered by using the same searching protocol to investigate the membrane association of other peptides. Our approach was similar to that of Efremov et al. (1999a,b, 2002a,b), inasmuch as we verified that, when each of the free-energy terms in Eq. 1 is taken separately, it reproduced the appropriate experimental data. The DGcon component, which was essentially derived from the available proteins of known three-dimensional structure, was shown to be potent in reproducing experimental data on the folding and stability of a series of polyalanine-like peptides (ShentalBechor, Kirca, Ben-Tal, and Haliloglu, unpublished data). Similarly, the membrane-related terms (DGSIL), which were derived from continuum-solvent model calculations of the transfer free energy of each of the amino acids from the aqueous phase into the membrane, were shown to be potent in detecting the TM helices of proteins of known threedimensional structure (Chen et al., 2002). The potential we employed favors the formation of the peptide backbone hydrogen bonds in a direct, and an indirect, manner. Hydrogen-bonding is beneficial due to favorable internal energy contributions, and to the fact that by assuming helical conformations, which involve satisfied backbone hydrogen bonds, the system can escape a heavy penalty resulting from the desolvation terms. Indeed, as in similar methods (Milik and Skolnick, 1993, 1995; Baumgartner, 1996; Ducarme et al., 1998; Efremov et al., 1999a,b, 2002b; Lins et al., 2001; Maddox and Longo, 2002) our potential promotes helix formation in the membrane. Nevertheless, differences in the helicity of M2d between the TM and surface orientations were observed in this study. The hydrophobic environment of the membrane strongly imposed helicity in the TM orientation and significantly limited the conformations accessible to the peptide (Fig. 6 B). In contrast, several conformations with various degrees of helicity were observed in the surface-adsorbed orientation. The fact that stable, nonhelical conformations have been observed in the simulations suggests that the potential of Eq. 1 does not overimpose helicity. One of the novelties of this study is that peptide-induced membrane deformation effects were included in the MC simulations. This was done by the inclusion of the DGdef harmonic component in the potential as described in Methods. The magnitude of this term depends on the values chosen for the unperturbed width of the lipid bilayer and the force constant, reflecting the membrane response to changes in its width. Both of these values may differ, depending on the lipid composition of the membrane. In this study we Biophysical Journal 85(6) 3431–3444 3440 chose values that are typical for biological membranes (White and Wimley, 1999; Nagle and Tristram-Nagle, 2000). The same approach was successfully used in continuumsolvent model studies of the interactions between different peptides and lipid bilayers (e.g., Kessel et al., 2000a,b; Bransburg-Zabary et al., 2002). Our simulations suggest the formation of a helical-hairpin structure during peptide insertion into the lipid bilayer (see further discussion below). The insertion of the helicalhairpin into the lipid bilayer is likely to disturb the organization of the latter. In our simulations, membrane perturbation effects were taken into account in the DGSIL term, which is based on the calculated hydrophobicity scale of Kessel and Ben-Tal (2002). However, the values of the scale correspond to the perturbation of a membrane by a single, narrow and rigid a-helix that interacts with the two leaflets, whereas the hairpin, which is wider and most likely less rigid, spans only one leaflet of the bilayer. Since the value of DGlip is small compared to other free-energy terms, even a significant change in its magnitude should not drastically affect the results. The systematic approach in the construction of the potential in Eq. 1 enabled us to provide an estimate of the total free energy of membrane association as well as an energy breakdown into components. However, it is important to take note that this systematic approach does not fully guarantee that these terms perform well in concert. For example, it may well be that the various free-energy components do not add up to produce a well-balanced, total transfer free energy; that is, unit conversion factors may be missing. In addition there are some cases in which parts of the energy components may be counted several times. For example, in the case of large TM proteins, where residues at the protein core are tightly packed, the hydrophobicity may be accounted for both in the internal side-chain-to-side-chain interaction terms and the external solvation term. Such ‘‘overcounting’’ should have a minor effect in this study, since short peptides such as M2d are completely surrounded by the solvent. M2d is insoluble in water, and tends to aggregate above a certain critical concentration. Our simulations include only one M2d molecule and hence correspond to infinitely dilute aqueous solutions. This complicates the comparison of our results to experimental data. Experimental studies of M2d’s association with lipid bilayers (e.g., Opella et al., 1999) require its solubilization in organic solvents, which are missing in the simulations. This suggests that the pathway of membrane association and insertion of M2d in typical experiments is most likely different than the one observed in our simulations. Comparison with other MC methods In general, our approach is closest to those of Milik and Skolnick (1993, 1995) and Maddox and Longo (2002). As in these studies, we used a residue-level resolution. However, Biophysical Journal 85(6) 3431–3444 Kessel et al. we used two interaction sites for each amino acid, rather than only one. This is a more realistic representation of the physicochemical nature of the amino acids than single interaction sites. Furthermore, it enables hydrophobic residues, such as Leu, to assume conformations in which the side chain interacts favorably with the hydrocarbon core of the membrane, whereas the polar backbone partitions into the membrane-water interface. The treatment of backbone hydrogen-bonding in our model is significantly different than that of Milik and Skolnick (1993, 1995) and Maddox and Longo (2002). In contrast to these two methods, the favorable internal energy contribution from hydrogen-bond formation in our model may change depending on the identity of the residue. Previous MC studies provided only qualitative, rather than quantitative data on the peptide-membrane system (Milik and Skolnick, 1993, 1995; Baumgartner, 1996; Ducarme et al., 1998; Sintes and Baumgartner, 1998; Lins et al., 2001; Maddox and Longo, 2002). The method of Efremov et al. (1999a,b, 2002a,b) is an exception, as the free energy of peptide-membrane association is often provided. However, these values are usually unrealistically large in magnitude. This is presumably an inherent property of the internal energy term used in their potential, since the value reported recently for the membrane adsorption of two cardiotoxins, where only small conformation changes were recorded upon membrane association, appears to be reasonably low in magnitude (Efremov et al., 2002a). The main advantage of our method over similar MC simulations is that it provides quantitative information. Average values of the free energy of peptide-membrane association were provided, and are comparable to those measured in similar systems. The corresponding standard deviations provide an indication of the stability and convergence of the simulations. The model also accurately reproduced the helical conformation of M2d and its tilt in the membrane in the TM orientation. A more detailed comparison of our results with the measurements is provided below. Comparison with continuum-solvent model calculations In our previous studies (e.g., the accompanying article), we described the lipid bilayer as a low-dielectric slab, embedded in a high dielectric medium (i.e., the aqueous solution). This simplistic model of the bilayer has been shown in our studies with different membrane-spanning peptides to have the capacity to describe the main thermodynamic and kinetic properties of the peptide-membrane system (Kessel et al., 2000a,b; Bransburg-Zabary et al., 2002; Kessel and Ben-Tal, 2002). However, it could not capture various aspects related to the interaction of the peptides with the water-membrane interface. This problem was addressed in the present simulations. Following the approach of Baumgartner (1996), Biggin et al. (1997), Ducarme et al. (1998), La Rocca et al. Monte Carlo Simulations with M2d (1999), Maddox and Longo (2002), and Milik and Skolnick (1993), we described the lipid bilayer using a polarity profile, which allowed a mean field description of the polar headgroup region of the lipid bilayer and its effect on the association of M2d with the membrane. The average TM orientation of M2d observed in the simulations is such that the polar termini of the peptide reside in the polar headgroup region of the membrane (Fig. 7), which is missing in the slab model used in the accompanying article. This accounts for the difference in the free energy of association with the membrane, as suggested by the two models. That is, the association free energy was less negative when the polarity of the lipid headgroups was considered, as compared to the value obtained from the slab model (Table 2 here versus Table 2 in accompanying article). The reason for the difference is that in the modified model but not in the slab model some of the polar groups at the peptide termini were exposed to regions of intermediate polarity. That is, these regions are more polar than the hydrocarbon region, but not as polar as the aqueous phase. The resulting electrostatic transfer free-energy penalty makes the solvation free energy less negative. FIGURE 7 The average position of M2d in the lipid bilayer in the TM orientation. In the simulations, each residue in the peptide was represented by two interaction centers (backbone and side chain). Here, each residue is presented as a single sphere, for clarity. The spheres are located at the Ca nuclei. The identities of the residues are noted by one letter code. The peptide is colored according to the polarity of the environment of each residue in the average position. The color code appears in the scale at the right side of the figure. 3441 The difference between the two models also accounts for the observed differences in the membrane curvature in the two studies: the continuum-solvent calculations suggested a reduction of ;10 Å in the width of the lipid bilayer in response to peptide insertion, whereas the MC simulations suggest a more likely value of 7 Å (Table 2). Again, this is most likely due to the fact that the interactions of M2d’s termini with the polar headgroups of membrane lipids are included in the MC, but not in the continuum-solvent model studies. Our MC simulations use an implicit description of the peptide, where each residue is represented by two interaction sites. The implicit description makes the simulations computationally feasible. However, it involves inaccuracies in the calculations, especially those related to the solvation free energy, which strongly depends on the location of each atom (particularly those that are polar). Conversely, the continuum-solvent model, in which the peptide is described explicitly, emphasizes the consideration of such effects. Insertion pathway The simulations show that M2d association with the lipid bilayer begins with adsorption on the membrane surface (Fig. 3). The following step involves rearrangement into a partially helical structure, which is then inserted into the membrane (Fig. 4). This insertion pathway is in accordance with the model proposed by Jacobs and White (1989). Our continuum-solvent model calculations (see accompanying article), and the studies of Milik and Skolnick (1993, 1995) and Maddox and Longo (2002) give further support to the Jacobs and White model. In our simulations, M2d insertion into the lipid bilayer involves the formation of a helical-hairpin intermediate. This intermediate is formed independently of the starting conformation and orientation of the peptide with respect to the bilayer. The helical-hairpin intermediate, first suggested by Engelman and Steitz (1981), is consistent with the aminoacid sequence of M2d and other membrane-spanning peptides, which are typically composed of a hydrophobic core with polar terminal segments. It was also observed in other MC studies (e.g., Milik and Skolnick, 1993; Maddox and Longo, 2002). In M2d, the central region is overall hydrophobic, but includes a few polar residues (e.g., S8, Q13). In the helicalhairpin conformation, these polar residues are transferred into the hydrocarbon region of the bilayer, which is energetically unfavorable. The inherent instability of the hairpin conformation inside the bilayer encourages major structural fluctuations that ultimately lead to membrane translocation of one of the polar termini, such that the peptide resides in a TM orientation. It should be noted that S8 and Q13 are also exposed to the hydrocarbon region of the bilayer in TM orientations. However, in the TM orientations, the destabilizing effect due to the lipid exposure of these two Biophysical Journal 85(6) 3431–3444 3442 residues is overcompensated by stabilizing internal energy and solvation effects due to a-helix formation. Recent MC simulations by Maddox and Longo (2002) suggested another insertion pathway. In this alternative pathway, the peptide first assumes a straight helical conformation along the membrane-water interface. One of its termini becomes inserted into the hydrocarbon region of the membrane, and finally sprouts from the other end. Such a pathway was never observed in our simulations and since the potential used in our simulations differs very much from that of Maddox and Longo, it is difficult to suggest a possible reason. We repeated the simulations with a modified C-terminus; residue R23 was assigned half its desolvation free-energy penalty. Although the peptide structure does not justify this modification, we used it to see whether the simulations were sensitive to the modification. Indeed, the modification led to the insertion of M2d into the membrane from its Cterminus, in contrast to insertion through the N-terminus, which was observed using the regular protocol (Fig. 4). This is a further demonstration of the dominant effect of the termini on peptide insertion. It also demonstrates that the revised MC scheme, which operates on a much less rough energy surface compared to the regular scheme, is still sensitive enough to distinguish between the energetic modifications at the termini, as in all cases the peptide prefers to insert into membrane from the modified terminus. Thermodynamic properties The average free-energy values obtained in our simulations demonstrate the TM orientation to be most favorable for M2d. This state is favored over the surface orientation by both the solvation (DGsol) and internal (DE) components of the total free energies. However, it involves the DGdef penalty, since peptide insertion into the membrane introduces local thinning of the membrane to match with the length of the peptide’s hydrophobic core, and a high entropy penalty. The conformation entropy penalty (DScon) of the TM orientation is larger than that of the surface orientation; in the TM orientation the virtual bonds are confined to a much smaller fraction of the conformation space as compared to the surface orientation. The two membrane-associated orientations appear different, both by estimation of the entropy and by cluster analysis of the conformations generated during the simulations (Fig. 6, A and B). DGcon values of
5.4 and
3.3 kT were obtained for membrane association of M2d in the TM and surface orientations, respectively. These values are very similar to the Zimm-Bragg value of ;
4 kT (Zimm and Bragg, 1959; Lifson and Roig, 1961; Chakrabartty and Baldwin, 1995; Scholtz and Baldwin, 1992) that was used in the companion article. Biophysical Journal 85(6) 3431–3444 Kessel et al. On average, the TM orientation is favored over the surface orientation by 3 kT, which is significant in itself, but is smaller in magnitude than the standard deviations in the calculated free-energy values of the two peptide-membrane configurations. In particular, the fluctuations in the total freeenergy values of the surface orientation are larger in magnitude than the average value itself:
7.4 (613.4) kT. The free energy of the TM orientation has a lower variation:
10.2 (65.6) kT. The large variation in the surface orientation is a result of the high flexibility of the peptide in this orientation, which is reflected in high variability of the internal energy, the conformational entropy, and the solvation components of the total free energy. Our previous studies of peptide-membrane interactions (e.g., Kessel et al., 2000a,b; Bransburg-Zabary et al., 2002) have demonstrated the dominance of solvation effects in these systems. The free-energy change due to conformation effects (DGcon) was assumed to be negligible compared to the solvation free energy, mainly because conformation changes, which accompanied the transfer of the peptide from the aqueous solution into the lipid bilayer, were confined to small regions on the peptide. M2d is likely to go through significant conformation changes during this transfer, which suggests comparable free-energy changes. Indeed, our simulations suggest that, in the M2d-membrane system, the value of DGcon is similar in magnitude to the value of DGSIL. DGdef, which is the only term disfavoring the transfer of the peptide into the membrane, has a significant value as well, resulting from the large deformation observed in this system. Our simulations demonstrate that, on average, M2d assumes a TM a-helix conformation with a slight tilt with respect to the membrane normal. The TM conformations produced during the simulations are very close to the NMR structures obtained by Opella et al. (1999), with root-meansquare deviations of \3 Å (usually \2 Å; see Table 3). The observed tilt of 148 (Table 2) compares well with the value of 158 that was calculated using the continuum-solvent model (see accompanying article) and the value of 128, which was measured in lipid bilayers using solid-state NMR (Opella et al., 1999). The partitioning of M2d between the aqueous phase and lipid bilayers was not measured, but the calculated water-to-membrane transfer free energy of the peptide (DGtot ¼
10.2 kT) compared well with values that were measured for similar peptides (White and Wimley, 1999). The overall good agreement between simulations and experiments suggests the potency of the methodology for studying the interactions between hydrophobic peptides and membranes. This work was supported by the Magnet ‘‘Pharmalogica’’ Consortium of the Israel Ministry of Industry and Trade and by North Atlantic Treaty Organization grant LAST-CLG 977836. T.H.’s research is supported by State Planning Organization of Turkish Republic grant 01K120280 and Bogazici University Research grant 02HA501, EA-TUBA-GEBIP-2001-1-1. Monte Carlo Simulations with M2d 3443 REFERENCES Honig, B. H., and W. L. Hubbell. 1984. Stability of ‘‘salt bridges’’ in membrane proteins. Proc. Natl. Acad. Sci. USA. 81:5412–5416. Bahar, I., B. Erman, T. Haliloglu, and R. L. Jernigan. 1997a. Efficient characterization of collective motions and inter-residue correlations in proteins by low-resolution simulations. Biochemistry. 36:13512–13523. Honig, B., and A. Nicholls. 1995. Classical electrostatics in biology and chemistry. Science. 268:1144–1149. Bahar, I., and R. L. Jernigan. 1996. Coordination geometry of nonbonded residues in globular proteins. Fold. Des. 1:357–370. Bahar, I., M. Kaplan, and R. L. Jernigan. 1997b. Short-range conformational energies, secondary structure propensities, and recognition of correct sequence-structure matches. Proteins. 29:292–308. Baumgartner, A. 1996. Insertion and hairpin formation of membrane proteins: a Monte Carlo study. Biophys. J. 71:1248–1255. Bechinger, B., Y. Kim, L. E. Chirlian, J. Gesell, J. M. Neumann, M. Montal, J. Tomich, M. Zasloff, and S. J. Opella. 1991. Orientations of amphipathic helical peptides in membrane bilayers determined by solidstate NMR spectroscopy. J. Biomol. NMR. 1:167–173. Ben-Shaul, A., N. Ben-Tal, and B. Honig. 1996. Statistical thermodynamic analysis of peptide and protein insertion into lipid membranes. Biophys. J. 71:130–137. Jacobs, R. E., and S. H. White. 1989. The nature of the hydrophobic binding of small peptides at the bilayer interface: implications for the insertion of transbilayer helices. Biochemistry. 28:3421–3437. Kessel, A., and N. Ben-Tal. 2002. Free energy determinants of peptide association with lipid bilayers. In Peptide-Lipid Interactions. Current Topics in Membranes. S. A. Simon and T. J. McIntosh, editors. Academic Press, San Diego. Kessel, A., D. S. Cafiso, and N. Ben-Tal. 2000a. Continuum solvent model calculations of alamethicin-membrane interactions: thermodynamic aspects. Biophys. J. 78:571–583. Kessel, A., K. Schulten, and N. Ben-Tal. 2000b. Calculations suggest a pathway for the transmembrane migration of a hydrophobic peptide. Biophys. J. 79:2322–2330. Kurt, N., and T. Haliloglu. 1999. Conformational dynamics of chymotrypsin inhibitor 2 by coarse-grained simulations. Proteins. 37:454–464. Ben-Tal, N., A. Ben-Shaul, A. Nicholls, and B. Honig. 1996. Free-energy determinants of a-helix insertion into lipid bilayers. Biophys. J. 70: 1803–1812. Kurt, N., and T. Haliloglu. 2001. Conformational dynamics of subtilisinchymotrypsin inhibitor 2 complex by coarse-grained simulations. J. Biomol. Struct. Dyn. 18:713–731. Biggin, P. C., J. Breed, H. S. Son, and M. S. Sansom. 1997. Simulation studies of alamethicin-bilayer interactions. Biophys. J. 72:627–636. Kyte, J., and R. F. Doolittle. 1982. A simple method for displaying the hydropathic character of a protein. J. Mol. Biol. 157:105–132. Bransburg-Zabary, S., A. Kessel, M. Gutman, and N. Ben-Tal. 2002. Stability of an ion channel in lipid bilayers: implicit solvent model calculations with gramicidin. Biochemistry. 41:6946–6954. La Rocca, P., Y. Shai, and M. S. Sansom. 1999. Peptide-bilayer interactions: simulations of dermaseptin B, an antimicrobial peptide. Biophys. Chem. 76:145–159. Chakrabartty, A., and R. L. Baldwin. 1995. Stability of a-helices. Adv. Protein Chem. 46:141–176. Chen, C. P., A. Kernytsky, and B. Rost. 2002. Transmembrane helix predictions revisited. Protein Sci. 11:2774–2791. Law, R. J., L. R. Forrest, K. M. Ranatunga, P. La Rocca, D. P. Tieleman, and M. S. Sansom. 2000. Structure and dynamics of the pore-lining helix of the nicotinic receptor: MD simulations in water, lipid bilayers, and transbilayer bundles. Proteins. 39:47–55. Dan, N., and S. A. Safran. 1998. Effect of lipid characteristics on the structure of transmembrane proteins. Biophys. J. 75:1410–1414. Levy, R. M., and M. Karplus. 1979. Vibrational approach to the dynamics of an a-helix. Biopolymers. 18:2465–2495. Ducarme, P., M. Rahman, and R. Brasseur. 1998. IMPALA: a simple restraint field to simulate the biological membrane in molecular structure studies. Proteins. 30:357–371. Lifson, S., and A. Roig. 1961. On the theory of helix-coil transition in polypeptides. J. Chem. Phys. 34:1963–1974. Efremov, R. G., D. E. Nolde, G. Vergoten, and A. S. Arseniev. 1999a. A solvent model for simulations of peptides in bilayers. I. Membranepromoting a-helix formation. Biophys. J. 76:2448–2459. Efremov, R. G., D. E. Nolde, G. Vergoten, and A. S. Arseniev. 1999b. A solvent model for simulations of peptides in bilayers. II. Membranespanning a-helices. Biophys. J. 76:2460–2471. Efremov, R. G., P. E. Volynsky, D. E. Nolde, P. V. Dubovskii, and A. S. Arseniev. 2002a. Interaction of cardiotoxins with membranes: a molecular modeling study. Biophys. J. 83:144–153. Efremov, R. G., P. E. Volynsky, D. E. Nolde, A. van Dalen, B. de Kruijff, and A. S. Arseniev. 2002b. Monte Carlo simulations of voltage-driven translocation of a signal sequence. FEBS Lett. 526:97–100. Engelman, D. M., and T. A. Steitz. 1981. The spontaneous insertion of proteins into and across membranes: the helical hairpin hypothesis. Cell. 23:411–422. Fattal, D. R., and A. Ben-Shaul. 1993. A molecular model for lipid-protein interaction in membranes: the role of hydrophobic mismatch. Biophys. J. 65:1795–1809. Flory, P. J. 1969. Statistical Mechanics of Chain Molecules. WileyInterscience, New York. Haliloglu, T., and I. Bahar. 1999. Coarse-grained simulations of conformational dynamics of proteins: application to apomyoglobin. Proteins. 31:271–281. Hao, M. H., and H. A. Scheraga. 1995. Statistical thermodynamics of protein folding: comparison of a mean-field theory with Monte Carlo simulations. J. Chem. Phys. 102:1334–1348. Helfrich, P., and E. Jakobsson. 1990. Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys. J. 57:1075–1084. Lins, L., B. Charloteaux, A. Thomas, and R. Brasseur. 2001. Computational study of lipid-destabilizing protein fragments: towards a comprehensive view of tilted peptides. Proteins. 44:435–447. Maddox, M. W., and M. L. Longo. 2002. A Monte Carlo study of peptide insertion into lipid bilayers: equilibrium conformations and insertion mechanisms. Biophys. J. 82:244–263. May, S., and A. Ben-Shaul. 1999. Molecular theory of lipid-protein interaction and the La-HII transition. Biophys. J. 76:751–767. Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. J. Teller. 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21:1087–1092. Milik, M., and J. Skolnick. 1993. Insertion of peptide chains into lipid membranes: an off-lattice Monte Carlo dynamics model. Proteins. 15: 10–25. Milik, M., and J. Skolnick. 1995. A Monte Carlo model of FD and PF1 coat proteins in lipid membranes. Biophys. J. 69:1382–1386. Mouritsen, O. G., and M. Bloom. 1984. Mattress model of lipid-protein interactions in membranes. Biophys. J. 46:141–153. Nagle, J. F., and S. Tristram-Nagle. 2000. Structure of lipid bilayers. Biochim. Biophys. Acta. 1469:159–195. Nielsen, C., M. Goulian, and O. S. Andersen. 1998. Energetics of inclusioninduced bilayer deformations. Biophys. J. 74:1966–1983. Opella, S. J., F. M. Marassi, J. J. Gesell, A. P. Valente, Y. Kim, M. OblattMontal, and M. Montal. 1999. Structures of the M2 channel-lining segments from nicotinic acetylcholine and NMDA receptors by NMR spectroscopy. Nat. Struct. Biol. 6:374–379. Scholtz, J. M., and R. L. Baldwin. 1992. The mechanism of a-helix formation by peptides. Annu. Rev. Biophys. Biomol. Struct. 21:95–118. Sintes, T., and A. Baumgartner. 1998. Membrane-mediated protein attraction. A Monte Carlo study. Physica A. 249:571–575. Biophysical Journal 85(6) 3431–3444 3444 Skolnick, J., and A. Kolinski. 1999. Monte Carlo approaches to the protein folding problem. Adv. Chem. Physics. 105:203–242. White, S. H., and W. C. Wimley. 1999. Membrane protein folding and stability: physical principles. Annu. Rev. Biophys. Biomol. Struct. 28: 319–365. Biophysical Journal 85(6) 3431–3444 Kessel et al. Zhang, Y., D. Kihara, and J. Skolnick. 2002. Local energy landscape flattening: parallel hyperbolic Monte Carlo sampling of protein folding. Proteins. 48:192–201. Zimm, B. H., and J. K. Bragg. 1959. Theory of the phase transition between helix and random coil in polypeptide chains. J. Phys. Chem. 31:526–535.