Homepage The model Output FAQ Citing & Credits

MCPep is based on the Monte Carlo simulation model presented in references (1, 2) and here:

Peptide representation

Membrane representation

Calculation of ΔGtotal

Calculation of ΔGcon

Calculation of ΔGdef

Calculation of ΔGCoul

Calculation of ΔGSIL, ΔGimm and ΔGlip

Generation of conformations

Generation of configurations



Peptide representation

Each residue i is represented by two interaction sites corresponding to its α-carbon atom (Ciα) and its side chain interaction center, Si (Fig. 1) (1). These sites are selected in advance on the basis of the structure and energy characteristics of the amino acid (3). The peptide backbone is represented by virtual bonds connecting consecutive α-carbon atoms, as proposed by Flory and colleagues (4); a peptide of n residues has N-1 virtual bonds. Virtual bonds are highly stiff, and their lengths are taken here as fixed at their equilibrium values of 3.81±0.03Å. The peptide backbone conformation is defined by the 2N-5 dimensional vector [θ2, θ3, ..., θn-1, φ3, φ4, ..., φN-1] including n-2 virtual bond angles (θi) and n-3 dihedral angles (φi). The distance between Si and Ciα, as well as the side chain virtual bond vector pointing from Ciα to Si named θiS, are fixed at their equilibrium values. Thus, the conformation of side chain i is expressed by the torsion angle (φiS).


Figure 1. Schematic representation of the virtual bond model. A segment between backbone units Ci-2α and Ci+1α is shown. The site attached to the ith α-carbon is marked as Si. φi is the rotational angle of the ith virtual bond, connecting Ciα and Ci-1α. θi is the bond angle between virtual bonds i and i+1. θiS is the side chain virtual bond vector pointing from Ciα to Si. φiS is defined by Ci-2α, Ci-1α, Ciα and Si.

(back to top)

Membrane representation

The membrane is described using two parameters: hydrophobicity and surface charge. The hydrophobicity of the membrane (p) is proportional to the distance (z) between the interaction site and the bilayer midplane according to the following sigmoidal function:


where η determines the sharpness of the hydrophobicity profile and is set in advance at η = 1Å-1
(1), and zm representes the width of the hydrophobic region of a membrane monolayer. Fig. 2 shows a representative hydrophobicity profile of a membrane with zm=15Å. The surface charge is located (zm +5)Å from the membrane midplane. (back to top)


Description: Graph1.tif

Figure 2. Membrane representation. p(z) was calculated using Eq.1, with zm=15Å and η=1Å-1 (solid line). The location of the surface charge is designated by the dashed lines; the hydrophobic region of the membrane is defined by the dotted lines.

(back to top)

Calculation of

The total free energy difference between a peptide in the aqueous phase and in the membrane () can be divided into several terms as follows (5, 6):


(back to top)

Calculation of

is the free energy change due to membrane-induced conformational changes in the peptide. At constant (absolute) temperature T, it can be calculated as follows:


where is the internal energy difference between the water- and membrane-bound states of the peptide. The internal energy is derived from a statistical potential based on available 3-dimensional (3D) protein structures (3, 7). The energy function assigns a score (energy) to each peptide conformation according to the conformation's abundance in the Protein Data Bank (PDB). Common conformations are assigned high scores (low energy), while rare conformations are assigned lower scores (higher energy).

ΔS refers to the entropy difference between the water and membrane-bound states, while the entropy (S) in each state is determined by the distribution of the virtual bond rotations in the reduced peptide representation. To this end, the rotation space of each virtual bond is divided into 72 discrete intervals of 5° each. The entropy is estimated using the familiar "P lnP" relation:


where pi,j is the probability of virtual bond j to be in the interval i. As the virtual bond rotations for the first and last two amino acids are not defined, these are omitted from the entropy calculation. The value of pi,j is estimated as


where ni,j is the fraction of conformations out of the total number of conformations N in which a virtual bond j is in the interval i. (back to top)

Calculation of

ΔGdef is the free energy penalty associated with fluctuations of the membrane thickness around its resting (native) value. Insertion of a rigid hydrophobic inclusion into a lipid bilayer may result in deformation of the lipid bilayer such that the width of the hydrocarbon region matches the hydrophobic length of the inclusion, following the mattress model (8). The deformation involves a free energy penalty, ΔGdef, resulting from the compression or expansion of the lipid chains. Different methods have been used to calculate ΔGdef for lipid bilayers composed of lipids of various types; these methods have yielded similar values (8-14). For the MC model, the estimation of Fattal and Ben-Shaul (9) was chosen to estimate ΔGdef. Their calculations were based on a statistical-thermodynamic molecular model of the lipid chains and the fit of harmonic potential, using the following form:


and z0 are the actual and native widths of a monolayer. ω is a harmonic force constant related to the membrane elasticity and is equal to ω= 0.22 kT/Å2 (9). (back to top)

Calculation of ΔGcoul

Δ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 φ (measured in units of kBT/e, where kB is the Boltzmann constant, T is the temperature, and e is the electron charge) depends on the peptide's distance from the membrane surface in an electrolyte solution (15). To this end, the solution is considered neutral, containing monovalent salt. The protonation state of the side chains of the titratable residues in the solution is set according to pH = 7.


where κ is the inverse of Debye length:


is the number of monovalent anions per unit volume in bulk, ε0 is the permittivity in vacuum, and εr is the dielectric constant in water (taken as 80). Φ is the potential on the plane of smeared charges; it depends on the membrane charge density σ and on the molarity of the solution [K+]:


A is the Avogadro number. σ is determined according to the fraction of anionic lipids in the membrane (denoted fa), the valence of the lipids' charge (Z) and the area occupied by one phospholipid in the membrane (A):


A titratable residue interacts Coulombically with the charged membrane only when the residue is in its charged form. However, due to the large desolvation free energy penalty associated with the transfer of a charge into hydrophobic environment, titratable residues typically are neutralized when approaching the nonpolar environment of the membrane. Therefore, the dependence of the charged state fraction of titratable residue i on its distance from the membrane midplane is arbitrarily described with a sigmoidal function χi(z), similar to the function for the membrane polarity profile p(z) (Eq. 1):


where η is the profile steepness, previously calibrated to 1
(2), and h is the distance between the membrane midplane and the torque point of the sigmoidal function, equal to (z0-2)Å (2). The electrostatic energy of interaction between titratable residue i and the charged membrane surface weighted by χi(z) can be calculated according to the following:


The energy of the electrostatic interaction of the whole peptide with the membrane is:


A full positive charge is assigned to the side-chain interaction site of Lys and Arg and to the α-carbon of the N-terminus. A full negative charge is assigned to the side-chain interaction site of Asp and Glu and to the α-carbon of an unamidated C-terminal. (back to top)

Calculation of ΔGsol, ΔGimm and ΔGlip

ΔGsol is the free energy of transfer of the peptide from water to the membrane. It accounts for electrostatic contributions resulting from changes in solvent polarity, as well as for nonpolar (hydrophobic) effects, which result both from differences in the 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 thickness.

The sum of ΔGsol, ΔGimm and ΔGlip is denoted as ΔGSIL, and in the MC simulations it is calculated by summing over the contributions of the individual amino acids. The free energy of transfer of a specific residue i, i.e., ΔgSIL_i(z) can be decomposed into contributions from its backbone ΔgbSIL and side chain ΔgsSIL. Thus, ΔGSIL can be determined according to the following:


The prefactors in the two terms ( and ) represent the hydrophobicity of the environments at the respective interaction sites, i.e., the α-carbon site and the side chain interaction center, Si. These values are calculated according to Eq. 1.

Δgis are calculated using the Kessel and Ben-Tal hydrophobicity scale (Table 1) (5). 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. In order to avoid the excessive penalty associated with the transfer of charged residues into the bilayer, in the model the titratable residues are neutralized gradually upon insertion into the membrane, so that a nearly neutral form is desolvated into the hydrophobic core. However, as described above, a gradual transition between the charged and neutral forms of titratable residues based on χi(z) (Eq.11) is introduced into the model. Therefore, for the neutral state of a titratable residue Δgis is derived from the hydrophobicity scale (Table 1); for the charged state of a titratable residue Δgis is taken as 64 kT (16):


is the polarity profile of the charged side chains solvation, a sigmoidal function similar to the hydrophobicity profile (Eq. 1):


The free energy of transfer of the peptide backbone from the aqueous phase into the membrane (Δgib) is calculated using the following set of equations:


where f is proportional to backbone deviations from the optimal α-helical conformation as observed by Bahar and Jernigan (3). It is assigned a value of zero for residues in their ideal α-helical conformations obtained at φ0=-120° (7); for these residues, the C=O and N-H backbone groups neutralize each other. The value of 1 is assigned to f for residues deviating significantly from the ideal α-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 are taken into account:


Obviously, the presence of Eq. 18 in the potential strongly enhances the formation of helical structures upon membrane association. γ is the standard deviation of the distribution of Φ around its optimal value characterizing α-helical conformations, estimated as λ=30°

It is noteworthy that the stretches of three residues at the N- and C-termini are treated differently than the peptide core (Eq. 17). The free energy penalty associated with the transfer of the uncompensated hydrogen bonds of the N-H groups of the three residues at the N-terminus is taken into account regardless of the peptide conformation. Likewise, the free energy penalty due to the transfer of the uncompensated hydrogen bonds of the C=O groups of the three residues at the C-terminus is also taken into account, regardless of the peptide conformation.


Amino acid

ΔGi (kcal/mol)














































Table 1. 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 (ΔgSIL_i). The scale was computationally derived, as described in Kessel and Ben-Tal (2002) (5). 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. (back to top)


Generation of conformations

New conformations are generated by simultaneous random perturbations of φi, θi, φis:


where r is a random number between 0 and 1; δk is the maximum variation of the respective coordinates: 3° for φi and 0.5° for θi and φis.

The peptide configuration is changed by external motions as described below. However, it is noteworthy that a set of randomly chosen conformational changes could also lead to slight changes in the peptide's orientation in the membrane. (back to top)


Generation of configurations

External rigid body rotational and translational motions arre carried out to change the peptide's configuration, i.e. its location in, and orientation with respect to, the membrane. These motions are represented respectively by:




where α, β, and γ are three Euler angles describing the orientation, and r represents the Cartesian coordinates of the peptide's geometric center. δαmax, δβmax, δγmax and δdmax (=5°, 5°, 5° and 0.02Å, respectively) were chosen to be maximum variations of the random perturbations of α, β, γ and d. These parameters were determined by trial and error, such that the system would have enough time for internal relaxation but would not be trapped too often in local energy minima.

Additionally, the width of the membrane's hydrophobic region is perturbed:


δz, the maximum variation of a perturbation, is set to 0.5Å. Overall, the membrane width zm is allowed to vary from its native value z0 by up to 20%
(17). (back to top)



To calculate the free energy difference between a peptide in the aqueous phase and the peptide in the membrane (Eq. 2), peptide simulations in water and in membrane environments are performed. A standard MC protocol is employed, and acceptance of each move is based on the Metropolis criterion and the free energy difference between the new and old states (18). In water, peptides are subjected solely to internal conformational modifications. In one MC cycle the number of internal modifications performed is equal to the number of residues in the peptide. Therefore, the acceptance criterion is based on ΔΔE (Eq. 3). In the membrane, each MC cycle includes additional external rigid body rotational and translational motions as described above. Thus, the acceptance criterion in this case is derived from the following free energy difference:


(back to top)



1. Kessel, A., D. Shental-Bechor, T. Haliloglu, and N. Ben-Tal. 2003. Interactions of hydrophobic peptides with lipid bilayers: Monte Carlo simulations with M2delta. Biophys J 85:3431-3444.

2. Shental-Bechor, D., T. Haliloglu, and N. Ben-Tal. 2007. Interactions of cationic-hydrophobic peptides with lipid bilayers: a Monte Carlo simulation method. Biophys J 93:1858-1871.

3. Bahar, I., and R. L. Jernigan. 1996. Coordination geometry of nonbonded residues in globular proteins. Fold Des 1:357-370.

4. Flory, P. J. 1969. Statistical Mechanics of Chain Molecules. Wiley-Interscience, New-York.

5. Kessel, A., and N. Ben-Tal. 2002. Free energy determinants of peptide association with lipid bilayers. In Current Topics in Membranes. T. J. M. Sidney A. Simon, editor. Academic Press. 205-253.

6. Ben-Tal, N., A. Ben-Shaul, A. Nicholls, and B. Honig. 1996. Free-energy determinants of alpha-helix insertion into lipid bilayers. Biophys J 70:1803-1812.

7. Bahar, I., M. Kaplan, and R. L. Jernigan. 1997. Short-range conformational energies, secondary structure propensities, and recognition of correct sequence-structure matches. Proteins 29:292-308.

8. Mouritsen, O. G., and M. Bloom. 1984. Mattress model of lipid-protein interactions in membranes. Biophys J 46:141-153.

9. 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.

10. Helfrich, P., and E. Jakobsson. 1990. Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys J 57:1075-1084.

11. 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.

12. Dan, N., and S. A. Safran. 1998. Effect of lipid characteristics on the structure of transmembrane proteins. Biophys J 75:1410-1414.

13. Nielsen, C., M. Goulian, and O. S. Andersen. 1998. Energetics of inclusion-induced bilayer deformations. Biophys J 74:1966-1983.

14. May, S., and A. Ben-Shaul. 1999. Molecular theory of lipid-protein interaction and the Lalpha-HII transition. Biophys J 76:751-767.

15. McLaughlin, S. 1989. The electrostatic properties of membranes. Annu Rev Biophys Biophys Chem 18:113-136.

16. Honig, B. H., and W. L. Hubbell. 1984. Stability of "salt bridges" in membrane proteins. Proc Natl Acad Sci U S A 81:5412-5416.

17. Marsh, D. 2008. Energetics of hydrophobic matching in lipid-protein interactions. Biophys J 94:3996-4013.

18. Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. 1953. Equation of State Calculations by Fast Computing Machines. J Chem Phys 21:1087-1092.