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

Each residue i is represented by two interaction sites corresponding to its α-carbon atom (C_i^α) and its side chain interaction center, S_i (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 [θ₂, θ₃, ..., θ_n_-1, φ₃, φ₄, ..., φ_N_-1] including n-2 virtual bond angles (θ_i) and n-3 dihedral angles (φ_i). The distance between S_i and C_i^α, as well as the side chain virtual bond vector pointing from C_i^α to S_i named θ_i^S, are fixed at their equilibrium values. Thus, the conformation of side chain i is expressed by the torsion angle (φ_i^S).

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

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:

(1)

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

Figure 2. Membrane representation. p(z) was calculated using Eq.1, with z_m=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.

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):

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 p_i,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 p_i,j is estimated as

where n_i,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)

ΔG_def 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, ΔG_def, resulting from the compression or expansion of the lipid chains. Different methods have been used to calculate ΔG_def 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 ΔG_def. Their calculations were based on a statistical-thermodynamic molecular model of the lipid chains and the fit of harmonic potential, using the following form:

z_m and z₀ 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/Å² (9). (back to top)

ΔG_coul 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 k_BT/e, where k_B 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.

c^b is the number of monovalent anions per unit volume in bulk, ε₀ 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⁺]:

N_A is the Avogadro number. σ is determined according to the fraction of anionic lipids in the membrane (denoted f_a), 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 (z₀-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)

ΔG_sol 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. ΔG_imm is the free energy penalty resulting from the confinement of the external translational and rotational motion of the peptide inside the membrane. ΔG_lip 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 ΔG_sol, ΔG_imm and ΔG_lip is denoted as ΔG_SIL, 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., Δg_{SIL_i}(z) can be decomposed into contributions from its backbone Δg^b_SIL and side chain Δg^s_SIL. Thus, ΔG_SIL 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, S_i. These values are calculated according to Eq. 1.

Δg_i^s 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 Δg_i^s is derived from the hydrophobicity scale (Table 1); for the charged state of a titratable residue Δg_i^s 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 (Δg_i^b) is calculated using the following set of equations:

(17)

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 φ₀=-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° (7).

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	ΔG_i (kcal/mol)
I	-2.6
L	-2.6
F	-1.5
V	-1.2
A	-0.2
G	0.0
C	+0.4
S	+0.8
T	+1.1
M	+1.3
W	+1.3
P	+2.8
Y	+4.3
Q	+5.4
H	+6.8
K	+7.4
N	+7.7
E	+9.5
D	+11.5
R	+19.8
N-H	+1.8
C=O	+2.5

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 (Δg_{SIL_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)

New conformations are generated by simultaneous random perturbations of φ_i, θ_i, φ_i^s:

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 φ_i^s.

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)

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 δd_max (=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.

δz, the maximum variation of a perturbation, is set to 0.5Å. Overall, the membrane width z_m is allowed to vary from its native value z₀ 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: