Monte Carlo simulations of peptide-membrane interactions with the MCPep web server.

The MCPep server (http://bental.tau.ac.il/MCPep/) is designed for non-experts wishing to perform Monte Carlo (MC) simulations of helical peptides in association with lipid membranes. MCPep is a web implementation of a previously developed MC simulation model. The model has been tested on a variety of peptides and protein fragments. The simulations successfully reproduced available empirical data and provided new molecular insights, such as the preferred locations of peptides in the membrane and the contribution of individual amino acids to membrane association. MCPep simulates the peptide in the aqueous phase and membrane environments, both described implicitly. In the former, the peptide is subjected solely to internal conformational changes, and in the latter, each MC cycle includes additional external rigid body rotational and translational motions to allow the peptide to change its location in the membrane. The server can explore the interaction of helical peptides of any amino-acid composition with membranes of various lipid compositions. Given the peptide's sequence or structure and the natural width and surface charge of the membrane, MCPep reports the main determinants of peptide-membrane interactions, e.g. average location and orientation in the membrane, free energy of membrane association and the peptide's helical content. Snapshots of example simulations are also provided.

INTRODUCTION

Biological membranes separate between the interior of cells and organelles and the exterior environment. The membranes are composed mostly of lipid and protein molecules of various types (1). Roughly speaking, the lipid molecules are packed heterogeneously and organized in a fluid mosaic, mutually affecting each other. Biological membranes and protein– and peptide–lipid interactions have been investigated extensively with a wide range of biophysical and computational techniques (2). The goal of these studies has been to decipher the mechanism of action of membrane-active peptides, such as antimicrobial and viral-fusion peptides, and to deduce general principles of folding, energetics, stability, structure and function of membrane proteins.

Several computational approaches are commonly applied for investigation of protein– and peptide–lipid interaction (3–5). Our group developed a model that performs Monte Carlo (MC) simulations of the interaction of helical peptides with lipid bilayers. The simulation model has been tested on a variety of antimicrobial peptides, including magainin2 (6), penetratin (6), melittin (7), NKCS (8) and novicidin (9). These studies have furthered our understanding of the action mechanisms and selectivity of the antimicrobial peptides. The simulations correlated very well with the empirical data, and in some cases, they were used to guide further experimental effort (7–9). Additionally, the model was used to investigate the channel-forming peptide M2δ from the acetylcholine receptor δ-subunit (10). The energetically preferable configuration and the 3-dimensional (3D) structure of the M2δ peptide in the membrane were correctly predicted. We also used the model to predict stable conformations of the outer membrane 45 (OM45) protein on interaction with the mitochondrial outer membrane (11).

Here, we present MCPep (http://bental.tau.ac.il/MCPep/): a web implementation of the MC simulation model. MCPep can be used to study short membrane-active peptides, such as antimicrobial and fusion peptides, short membrane proteins (see the Case Study described later in the text) or fragments of large membrane proteins. The server takes as input the peptide’s sequence or structure, fractions of acidic and zwitterionic lipids in the membrane, the native width of the hydrocarbon region of the membrane and the ionic strength of the aqueous solution. The output includes the peptide helicity and average orientation in the membrane, the free energy of its membrane-association and the decomposition of the free energy. MCPep also presents snapshots of example simulations. In addition, the conformations obtained from the simulations are clustered, and the centroid conformations of the largest clusters are provided in Protein Data Bank (PDB) format.

THE MC MODEL

The peptide is described in a reduced way; each amino acid is represented by two interaction sites, corresponding to its α-carbon and side chain (10). The interaction sites of sequential α-carbons are connected by virtual bonds. The membrane hydrophobicity is represented as a smooth profile, corresponding to the native thickness of the hydrocarbon region (10).

The total free energy difference between a peptide in the aqueous phase and in the membrane (ΔGtotal) is calculated as (12,13) follows: where ΔGcon is the free energy change owing to membrane-induced conformational changes in the peptide. At constant temperature (T), it can be calculated as where ΔE is the internal energy difference between the peptide in the aqueous phase and in the bilayer. The internal energy is a statistical potential derived from available 3-dimensional (3D) protein structures (14). The energy function assigns a score (energy) to each peptide conformation according to its abundance in the PDB. Common conformations are assigned high scores (low energy), whereas rare conformations are assigned lower scores (higher energy). ΔS refers to the entropy difference between the aqueous solution and membrane-bound states, whereas the entropy (S) in each state is determined by the distribution of the virtual bond rotations in the reduced peptide representation.

ΔGdef is the free energy penalty associated with fluctuations of the membrane thickness around its native (resting) value, calculated following the estimation of Fattal and Ben-Shaul (15). Their calculations were based on a statistical-thermodynamic molecular model of the lipid chains. Their model fits a harmonic potential of the form: where ΔL is the difference between the native and actual width of the membrane; ω is a harmonic-force constant related to the membrane elasticity and is equal to (15).

ΔGCoul stands for the Coulombic interactions between titratable residues of the peptide and the (negative) surface charge of the membrane. It is calculated using the Gouy-Chapman theory that describes how the electrostatic potential depends on the distance from the membrane surface in an electrolyte solution (6).

ΔGsol is the free energy of transfer of the peptide from the aqueous phase into the membrane. It takes into account electrostatic contributions resulting from changes in solvent polarity, as well as nonpolar effects, both resulting from differences in van der Waals interactions of the peptide with the membrane and aqueous phases and from solvent structure effects. ΔGimm is the free energy penalty resulting from the confinement of the external translational and rotational motion of the peptide inside the membrane. ΔGlip is the free energy penalty resulting from the interference of the peptide with the conformational freedom of the aliphatic chains of the lipids in the bilayer while the membrane retains its native width. The latter three terms, i.e. ΔGsol, ΔGimm, ΔGlip, collectively referred to as ΔGSIL are calculated using the Kessel and Ben-Tal hydrophobicity scale (12). The scale accounts for the free energy of transfer of the amino acids, located in the center of a polyalanine α-helix, from the aqueous phase into the membrane midplane. To avoid the excessive penalty associated with the transfer of charged residues into the bilayer, in the model, the titratable residues are neutralized gradually on insertion into the membrane, so that a nearly neutral form is desolvated into the hydrocarbon region.

To calculate the free energy change in Equation 1, the peptide is simulated in the aqueous phase and membrane environment. In the aqueous phase, the peptide is subjected only to internal conformational modifications. In one MC cycle, the number of internal modifications performed is equal to the number of residues in the peptide. In the membrane, each MC cycle includes additional external rigid body rotational and translational motions to allow the peptide to change its location in the membrane and its orientation with respect to the membrane normal. A standard MC protocol is used, and acceptance of each move is based on the Metropolis criterion and the free energy difference between the new and old states (16).

THE WEB SERVER

Input

The user is asked to provide the peptide sequence or upload an initial 3D structure in PDB format. The user should also specify the fraction of acidic lipids (mol percent) in the membrane and ionic strength of the aqueous solution. Additionally, the user may set the native width of the hydrocarbon core of the membrane according to the lipid types. To assist the user in doing so, the web server provides a link to recent X-ray studies of phospholipids [reviewed in reference (17)]. The user can also specify an e-mail address where a link to the results page will be sent once the simulations are completed (optional).

Processing

If a peptide sequence is provided as initial input, an initial canonical α-helix (phi = −58, psi = −47) model structure of the peptide is constructed, with side chains modeled using Scwrl 4.0 (18); otherwise, the structure provided by the user is taken as the initial conformation. For simulations in the aqueous phase (i.e. without the membrane), the initial structure is used as it is. The simulations are carried out in three to five independent runs of 500 000–900 000 MC cycles each; the recommended number and length of the runs are set as defaults but the user can use alternatives. In a membrane, a helical peptide typically adheres to one of the two following configurations: transmembrane (TM) orientation, with the principal axis of the helix positioned approximately along the membrane normal; or surface orientation, with the axis approximately parallel to the membrane surface. The transition between the two configurations is associated with a high desolvation free energy barrier (2). Thus, for simulations in the membrane environment, each of the two configurations is used as the initial orientation. The number of independent simulations in each configuration and the number of cycles in each individual simulation are the same as in simulations in the aqueous phase. However, it should be noted that a simulation that starts in TM orientation may end in surface orientation if the peptide is amphipathic, too short to span the bilayer or too polar. The transition from a surface to TM orientation is highly unlikely owing to the high free energy penalty associated with the insertion of the polar backbone of the helix termini into the membrane.

Output

The web server’s output includes the helical content of the peptide, i.e. the total number of residues that are in helix conformation, in the aqueous phase and in the membrane, calculated as described in reference (19). A residue is defined as being in a helical conformation if the rotational angle of each of its two adjacent bonds lies within the interval typical for a helix (−120° ± 30°) (14). As the rotational angles of the peptide’s termini cannot be defined, the two residues at the N-terminus and the two residues at the C-terminus are neglected. Energetically favorable orientations of the peptide in the membrane are reported, including the average locations of the peptide’s α-carbons. The free energy of membrane association and the decomposition of the free energy (Equation 1) are also presented. The membrane thickness is allowed to vary from its native value by up to 20%. The (perturbed) thickness along the simulation is also provided. In addition, snapshots of example simulations, both in the aqueous phase and in the membrane, are presented. The conformations obtained from the simulations are clustered with a preset cut-off for the root mean square deviation between the α-carbons. The default cut-off value is 3 Å, but the user is free to alter the value. The centroid conformation of the largest cluster from each simulation is reported (in PDB format), along with its average orientation in the membrane.

CASE STUDY

The synaptic vesicle protein synaptobrevin2 (syb2, also known as vesicle-associated membrane protein 2 (VAMP2)) is part of the so-called Soluble N-ethylmaleimide sensitive fusion protein Attachment Protein Receptor (SNARE) complex, responsible for fusion of neuronal vesicles with the presynaptic membrane in the neuron. Although the process is essential for synaptic neurotransmitter release, the exact mechanism of vesicle fusion is not yet clear. Here, we studied the membrane interaction of syb2 to illustrate the functionalities of the MCPep server. We used a peptide corresponding to the sequence of the putative TM and juxtamembrane regions of the protein, namely, residues 75–116 (UniProt id P63027). The fraction of the charged lipids in the membrane was set to 30%, similar to that of neuronal vesicles (20). The results showed that both TM and surface orientations of syb2 in the membrane are stable, with the former being somewhat more favorable energetically (Table 1). The helical content of the peptide increased on interaction with the membrane, as anticipated (Figure 1). The average TM and surface orientations of syb2, as predicted by MCPep, are presented in Figure 2. In both surface and TM configurations, the hydrophobic residues are inserted into the membrane core, whereas the hydrophilic residues partition into the aqueous phase (Figures 2 and 3).

Table 1.

Thermodynamic characteristics for syb2(75-116) in TM and surface configurations

Quality	TM orientation	Surface orientation
ΔGtotal (kT)	−53.2 ± 0.3	−44.4 ± 0.3
ΔGCon (kT)	−1.5 ± 0.3	5.3 ± 0.3
TΔS (kT)	−32.0 ± 0.0	−33.1 ± 0.3
ΔEint (kT)	−33.5 ± 0.3	−27.8 ± 0.3
ΔGSIL (kT)	−40.8 ± 0.1	−37.9 ± 0.1
ΔGdef (kT)	0.4 ± 0.00	0.3 ± 0.0
ΔGCoul (kT)	−11.3 ± 0.0	−12.1 ± 0.0
Mwidth (Å)	27.6 ± 0.0	30.1 ± 0.0
Zcenter (Å)	10.6 ± 0.0	18.7 ± 0.0
Tilt (°)	29.4 ± 0.2	73.2 ± 0.2

All values are reported as means ± standard error. The free energy terms are defined in Equation 1. Mwidth: the width of the membrane hydrophobic core. Zcenter: the average distance of the peptide’s center of mass from the membrane midplane. The Z-axis is the membrane normal, and the origin coincides with the membrane midplane. Tilt: the angle between the N'-to-C' vector of the peptide’s helical core and membrane normal.

Figure 1.

The average helical content of syb2 versus residue number as predicted by the MCPep server. (a) syb2 in the aqueous phase. (b) syb2 in TM configuration within a membrane containing 30% anionic lipids. The helical content increases on interaction with the membrane, as anticipated.

Figure 2.

The average location of the amino acids of syb2 in the membrane in (a) surface and (b) TM orientations. The membrane includes 30%-charged lipids. The horizontal dashed lines designate the location of the phosphate groups of the lipid polar heads. The hydrophobic residues (G, A, V, L, I, F) are in orange, polar residues (M, C, T, S, W, Y, H, Q, N) in green and charged amino acids (R, K, E) in blue.

Figure 3.

The centroid conformation of the largest cluster of syb2 in TM orientation. The peptide is colored according to the hydrophobicity scale. The average location of the phosphate heads of the two leaflets of the lipid bilayer is represented by the red rectangles.

We suggest that the exceptionally high stability of the TM configuration (−53kT; Table 1) is required for syb2’s function. Syb2 is located in vesicles, whereas its SNARE-complex partners, i.e. syntaxin and Soluble N-ethylmaleimide sensitive fusion protein Attachment Protein (SNAP)-25, are located in the presynaptic plasma membrane (21,22). Syntaxin and two molecules of SNAP-25 associate with syb2 via a process known as zippering, which starts at the mostly unstructured N-termini of the proteins. On association, the proteins form a remarkably stable tight bundle of four parallel coiled-coil α-helices, where each helix is from a different SNARE partner. As zippering continues towards the C-termini, syb2 pulls the vesicle towards the plasma membrane (21,22). This process eventually leads to the fusion of vesicular and plasma lipid bilayers. Clearly, vesicle pulling can only occur if both syb2-vesicle association and syb2/syntaxin/SNAP-25 interaction are strong. Indeed, the stability free energy of the SNARE complex is estimated to be between 18 (23) and 35 (24) kT, which is indicative of the strong forces required for membrane fusion to occur (21). In this respect, the strong anchoring of syb2 to the membrane (−53kT; Table 1) is compatible with these estimates. For comparison, we also simulated the interaction of the putative TM region of syntaxin (UniProt id Q16623, residues 247–288) with a membrane containing 30% charged lipids and obtained a similar value of −47kT. The free energy of membrane association of both proteins is comparable with the free energy required for hemi-fusion, that is the fusion of the outer leaflets of the bilayers, estimated at about 40–50kT (24).

DISCUSSION

MCPep allows rapid simulations of helical peptides in association with biological membranes. To minimize the computational burden, the server relies on reduced representation of the peptide. For similar reasons, it uses an implicit model of the membrane, which is described as a hydrophobic profile of preset thickness, and a surface charge (6,10). Thus, the computational model captures only certain characteristics of the complex peptide–lipid system. Clearly, other properties that might affect peptide–lipid interactions are missing in the model. One of these is the membrane curvature. The flat representation of the membrane surface can roughly describe a cell membrane or a large vesicle. However, the curvature of small vesicles or micelles, which are comparable in size to the peptide length, cannot be overlooked and might play a significant role in the interaction with peptides. Another example is the phase of the lipids, which affects the fluidity and rigidity of the bilayer. The MC model describes lipids in their crystalline-liquid phase, which are more susceptible to stretching and bending than gel-phase lipids (25). The liquid phase bilayer adjusts more easily to the presence of membrane proteins. Lipid de-mixing is another characteristic that is not incorporated into the model. In this process, which occurs in interactions with cationic peptides, negatively charged lipids migrate to the peptide–membrane interaction zone (26–28). As a consequence, the local charge density in the interaction zone increases, which may induce an increased local concentration of the peptide on the membrane’s surface.

The reduced representation also precludes investigation of specific interactions. For instance, it does not explicitly take into account hydrogen bonds and salt bridges between the side chains of polar amino acids or between the amino acids and the lipids’ polar heads. Thus, the model is not suitable to describe peptide–peptide interactions, and it does not explicitly account for packing of the lipid chains against the peptide. Yet, the use of MC (rather than molecular dynamics) simulations facilitates efficient sampling in conformational/configurational space, making this approach suitable for studying thermodynamic quantities. Indeed, the previous studies using our MC model have demonstrated that computations using the approximations described earlier in the text are feasible and that many important characteristics of peptide–membrane interactions are captured (6–11).

Many of the observations derived from the simulations depend on the water-to-membrane partition free energy of the amino acids, i.e. the hydrophobicity scale used in MCPep. We used the Kessel and Ben-Tal scale, derived from continuum (implicit)-solvent model calculations of the water-to-membrane transfer free energy of polyalanine-based peptides (12). The calculations were based on PARameter SEt (PARSE), a parameter set that was derived specifically to reproduce the experimental partition free energy of small molecules between water and liquid alkanes (29). The continuum-solvent model with PARSE has been used successfully to study thermodynamic and kinetic aspects of the transfer of several peptides and protein fragments into the lipid bilayer (10,13,30–33). The Kessel and Ben-Tal hydrophobicity scale was subsequently exploited within the MC simulation model (6–10). The scale’s hydrophobicity values are significantly lower in magnitude than those of other scales because the Kessel and Ben-Tal scale includes the free energy penalty owing to the transfer of the polar helix backbone from the aqueous phase into the membrane core (12). The scale is also unique in its extreme asymmetry: the free energy penalty associated with the water-to-membrane transfer of a single highly polar residue, such as Lysine in its neutral form, is approximately three times larger in magnitude than the gain owing to the transfer of a highly hydrophobic residue, such as Leucine (12).

PLANS

The methodology described herein can potentially be applied towards predicting TM protein structures. The most accurate models of protein structures are achieved through homology modeling, in which a model is inferred from a known high-resolution structure of a homologous protein (34). A crucial step in the modeling process is alignment of the query and template sequences (35), required to identify the TM helices of the query protein and match them to the template. Alignment is especially challenging if the query and target proteins share low sequence identity (36). The MCPep server could be used for identification of TM helices. In addition, it could predict the membrane boundaries for existing structures of TM proteins. Further investigation is required to explore these possibilities.

CONCLUSIONS

The MCPep server is capable of exploring the interactions of helical peptides of up to 50 residues with membranes of various types. Its results provide valuable mechanistic insights into membrane–peptide interactions, as demonstrated by the case study. It is preferable, when possible, to correlate the results with empirical data, and, thus, the conditions simulated in calculations should be similar to those used in the experiments. Model parameters that can be matched to experimental conditions include the fraction of acidic lipids in the membrane, the ionic strength of the surrounding aqueous solution and the native width of the membrane hydrocarbon core. The server can be easily handled even by non-expert users.

FUNDING

The Israel Science Foundation [1331/11 to N.B-T.] and NATO Science Programme Collaborative Linkage Grant (T.H. and N.B.-T.); The European Commission under the sixth Framework Program through the Marie-Curie Action: BIOCONTROL, contract number MCRTN—33439 to Y.G.; DPT (State Planning Organization, Turkey) Project No: 2009K120520 to T.H. Funding for open access charge: Prof. Nir Ben-tal, the university of Tel Aviv.

Conflict of interest statement. None declared.