Monte Carlo techniques for simulating in lipid bilayers

Pál Jedlovszky and Mihaly Mezei

Department of Physiology and Biophysics
Mount Sinai School of Medicine, NYU
New York, NY 10019

Originally presented in April, 1999 at a SCRI workshop at Tallahassee, Florida.

A short version of this work is published in the Journal of Chemical Physics, 111, 10770 (1999).
     Full text

I. Introduction

The determination of canonical and microcanonical ensemble averages in homogeneous disordered systems via computer simulation is by now a well established technique, using either the probabilistic Monte Carlo method or the deterministic molecular dynamics method. Furthermore, for lipid bilayers molecular dynamics methodology has been established in the early 90's by several groups (Sok and Berendsen, Stouch, Heller et al., Huang et al., Venable et al.). Both approaches calculate the Boltzmann-weighted average of quantities Q(XN), <Q>, via generation of a set of configurations with Boltzmann probability PB since a simple average of Q over these configurations yields its Boltzmann average <Q> .

The explicit time dependence of molecular dynamics made it the preferred technique in most applications. By mimicking the naturally occurring process in a sense guarantees convergence and thus gives an edge over Monte Carlo where time-dependent properties can, at best, be obtained indirectly. It is argued here, both from first principles and from additional examples, that Monte Carlo can gain an edge for certain problems over molecular dynamics since being tied to the natural time can also be a drawback. This point is illustrated by the comment in the paper of Tieleman and Berendsen describing a nanosecond-long simulation: "The possibility of trapping water molecules in unfavorable internal positions during the initialization procedure may be a serious cause of artifacts in simulations of biomolecular macromolecules". It should be stressed, however, that the term Monte Carlo covers a large class of methods and judicious choice among them is required for realizing any serious gain over molecular dynamics.

It will be demonstrated that the classical Metropolis method (Metropolis et al.) can be efficiently adopted to the simulation of lipid bilayers, provided the various parameters are judiciously optimized including the consideration of molecular conformation in a novel manner. Furthermore, the applicability and usefulness of cavity-biased grand-canonical ensemble technique (Mezei (1980,1989)) will also be shown.

II. Description of the methodology

II.1. Optimization of the Metropolis moves.

The modeling of the lipid bilayer requires the following types of attempted moves:

The relative frequencies of these moves were subject to optimization using shorter (200,000 MC steps long) runs.

II.1.a. Order of moves.

II.1.b. Stepsize selection

The regular Metropolis method selects a range for each degree of freedom sampled and keeps it fixed during the simulation. There is a class of methods, however, where this range is allowed to vary depending on the conformation relative to the original range selected (Goldman, Mezei et al.). In Sec. II.2 we will describe such a method for the rotations and torsions of a lipid. In the discussion below, the stepsize refers either to an actual stepsize or to the reference stepsize, around which the actual stepsizes vary.

The sampling efficiency was characterized by the mean overall displacement of the center of the lipid molecules during the run (averaged over all lipids), DL; the mean overall rotation of the molecule around the laboratory x, y, and z axes (the bilayer is in the y-z plane), Dalpha(x) Dalpha(y), Dalpha(z), respectively; and the mean overall displacement due to torsion angle changes of the entire molecule, the headgroup chain and the two hydrocarbon chains, DT, DTHG, DTHC1, DTHC2, respectively. These torsion-related displacements were obtained for each lipid molecule by first overlaying the fixed core of the lipid (the glycerol carbon the headgroup is attached to and its neighbors) in its initial and final conformations and calculating the root mean square differences of the atoms in each chain separately.

Finally, it should be noted that as long as only one torsion angle is changed at a time or the change in each torsion affects a disjoint set of atoms, the distance of moving atoms from their rotation axis remains unchanged, thus detailed balance is satisfied thus the acceptance probability is the same as for sampling with fixed stepsize. If the change in one torsion does affect atoms moved by a different torsion during the same move, the acceptance probability has to be scaled with the products of the ratios of stepsizes (for the forward and reverse moves) and it even has to be checked if the reverse move is possible at all, as described in Mezei et al.. We found, however, that it is best to change one torsion at a time, so the simulations described here all used the original Metropolis acceptance probability.

II.2. Extension-biased sampling

The angle with which a molecule (or a part thereof) can be moved without bumping into a neighboring molecule clearly depends on the shape of the molecule. For a given axis of rotation, the larger the moving atoms fall from the axis of rotation, the more displacement do they undergo for a given angle change. This observation lead us to the development of the extension-biased sampling, that modulates the stepsize parameter based on the the maximum distance seen between the axis of rotation and the moving atoms. The figure shows a lipid molecule (DMPC) with a torsion selected around the A-B bond, showing the atom C of maximum extension, calculated as the distance R-C where R is on the line of the A-B bond and R-C and R-A form a right angle. This can be calculated with little computational effort and without calculating the energy, thus the problem of outright rejection of some moves associated with the range modulation (discussed above) will not deteriorate the performance significantly.

DMPC molecule

The actual functional form of the scaling has been established by experimentation. We considered scaling with the reciprocal of the maximum extension, as well as with its square root. We also made comparisons between considering the maximum extension and the average extension. The best results were obtained using scaling with the reciprocal of the square root of the maximum extension: Alphamax=c/Rmax1/2.

II.3. Grand-canonical ensemble simulations

Following Adams, grand-canonical ensemble simulations select a parameter B related to the excess chemical potential µ' as

µ' = kT (B - ln <N>)
with N being the number of molecules. Attempts to insert (at a random position) or delete (a randomly chosen) molecule are accepted by probabilities
Pins = min { 1, exp[-(E(N+1)-E(N))/kT+B] / (N+1)}
Pdel = min { 1, N exp[-(E(N)-E(N-1))/kT-B]},
respectively. Random insertions into dense fluids is very unlikely to result in acceptance since the new molecule almost always overlaps with one of the existing ones.

The cavity-biased technique (Mezei (1980,1989)) has extended the range of usefulness of this methodology by attempting insertions only into cavities of suitable size. The bias thus introduced has to be removed, of course. This can be achieved by the modified acceptance probabilities

Pins = min { 1,Pcav(N) exp[-(E(N+1)- E(N))/kT+B] / (N+1)}
Pdel = min { 1, N exp[-(E(N)-E(N-1))/kT-B] / Pcav(N)}.
As described by Mezei (1989), a list of cavities can be maintained during the simulation with little computational cost. The same calculation also provides Pcav(N), the probability of finding a cavity, used in the modified acceptance expressions above.

Insertion and deletion steps were performed in alternating order after every attempt to move a water or a lipid. This added less than 20% to the computation time.

Earlier work showed that the cavity-biased method enables the GCE simulation of aqueous systems: liquid water (Mezei (1989)), crystal hydrates (Resat and Mezei) and solvated proteins (Marrone et al.)

III. Results

All calculations were performed on a bilayer of 50 dimyristoylphosphatidylcholine (DMPC) molecules surrounded by water. This simulation cell was a 72 Å long hexagonal prism, under periodic boundary conditions. The edge of the hexagon was established as 24.8 Å to reproduce an experimental estimate of the area per headgroup, 63.9 Å2 (De Young and Dill). The initial configuration was obtained from a system equilibrated by molecular dynamics (Duong et al), and reequilibrated after setting all bond lengths and bond angles to their equilibrium values.

III.1. Optimization results

The setup optimization proceeded in three stages. In the first stage, the lipid rotation range RL was varied and the rest of the parameters were kept fixed. This gave for the extension biased sampling the following result:

RL Pdacc Ptacc Dalpha(x) Dalpha(y) Dalpha(z) DL DT DTHG DTHC1 DTHC2
0.4 0.265 0.201 1.65 2.22 2.08 0.18 7.27 3.13 4.46 4.25
0.5 0.213 0.201 1.63 2.28 2.38 0.18 7.21 3.20 4.38 4.37
0.6 0.183 0.201 1.60 1.97 1.90 0.15 7.21 3.01 4.62 4.37
0.7 0.157 0.193 1.42 2.11 2.20 0.15 7.00 3.15 4.40 3.90

In the table above, Pdacc and Ptacc are the acceptance probabilities of the translaton/rotation and torsional steps, respectively. Selecting 0.5 as the best value, the optimization proceeded to establish the best ratio of of the lipid translation rotation attempt to the torsion change attempts, wR/wT:

wR/wT Pdacc Ptacc Dalpha(x) Dalpha(y) Dalpha(z) DL DT DTHG DTHC1 DTHC2
0.12 0.210 0.199 1.81 2.25 2.40 0.20 7.40 2.93 4.90 4.25
0.16 0.213 0.201 1.89 2.41 2.44 0.20 7.27 3.10 4.37 4.12
0.20 0.216 0.193 1.91 2.60 2.61 0.22 7.46 3.14 4.61 4.55
0.22 0.225 0.199 1.56 2.57 2.65 0.23 7.30 3.00 4.85 3.99
0.24 0.211 0.198 2.30 2.60 2.78 0.23 7.21 3.19 4.57 4.13
0.26 0.221 0.205 2.16 2.89 2.87 0.21 7.31 2.91 4.79 4.22

In the final round, using the best wR/wT value, 0.24, the factor (Tf) scaling all torsion reference stepsizes were optimized:

Tf Pdacc Ptacc Dalpha(x) Dalpha(y) Dalpha(z) DL DT DTHG DTHC1 DTHC2
0.08 0.206 0.246 2.55 2.99 2.69 0.21 7.15 2.48 4.59 4.46
0.10 0.211 0.198 2.30 2.60 2.78 0.23 7.21 3.19 4.57 4.13
0.11 0.211 0.181 2.23 3.06 3.19 0.22 8.01 3.29 4.94 4.82
0.12 0.211 0.164 2.01 2.57 2.34 0.22 6.92 3.64 4.16 3.64
0.14 0.215 0.145 1.89 2.57 2.58 0.22 7.30 3.32 4.18 4.35

It is clear from these tables that the selection of the "best" parameter set requires some judgment calls since the various convergence indicators are affected differently by these choices.

The optimal stepsizes and selection weights are given below.

Atom numbering scheme:

|                       |
|                       C10
Bond Selection weight Alphamax (deg) c (deg A1/2)Acceptance rate
28-251.0 3.150.4950.28
25-231.0 3.850.6050.25
23-201.0 3.850.6050.28
20-241.0 6.300.9900.25
24-171.0 10.51.6500.19
17-141.0 10.51.6500.20
14-10.2 14.02.2000.16
1-100.1 28.04.4000.29
1-60.1 28.04.4000.29
1-20.1 28.04.4000.29
28-301.0 5.550.7150.21
30-311.0 4.550.7150.22
31-331.0 7.01.1000.16
33-451.0 5.250.8250.26
45-481.0 8.751.3750.16
48-511.0 6.300.9900.25
51-541.0 10.501.6500.17
54-571.0 9.801.5400.20
57-601.0 14.002.2000.16
60-631.0 12.251.9250.21
63-661.0 17.502.7500.17
66-690.4 24.503.8500.14
69-720.4 28.004.4000.15
72-750.2 31.504.9500.15
75-780.1 38.506.0500.35
28-361.0 3.500.5500.25
36-391.0 4.200.6600.26
39-401.0 3.850.6050.25
40-421.0 4.900.9350.25
42-821.0 5.950.9350.21
82-851.0 6.300.9900.23
85-881.0 8.751.3750.18
88-911.0 8.401.3200.20
91-941.0 12.251.9250.16
94-971.0 12.251.9250.18
97-1001.0 15.052.3650.16
100-1031.0 21.003.3000.14
103-1060.4 27.753.5750.15
106-1090.4 28.004.4000.14
109-1120.2 31.504.9500.15
112-1150.1 42.006.6000.33

III.2. Extension-biasing results

The table below compares the various optimization targets using the fixed stepsize (UNB) and the extension-biased (EXB) sampling. The stepsize and relative sampling frequencies have been independently optimized from both methods.

Method Pdacc Ptacc Dalpha(x) Dalpha(y) Dalpha(z) DL DT DTHG DTHC1 DTHC2
EXB 0.211 0.181 2.23 3.06 3.19 0.22 8.01 3.29 4.94 4.82
UNB 0.410 0.314 1.82 2.43 2.60 0.25 6.68 2.64 4.44 3.92

The results show that extension biasing indeed improves the efficiency of rotational sampling significantly. Also, the improvement obtained by extension biasing in the sampling of the torsional space of the headgroup region is even larger than for the hydrocarbon chains. This is gratifying since the headgroup region, unlike that of the hydrocarbon tails, is largely inhomogenous.

It is also of interest that the acceptance rates found to be optimal are rather different for the two kinds of sampling techniques. This points out the importance of careful optimization of the setup parametrization - reliance of "rules of thumb" can lead to significantly less efficient simulations.

III.3. Overall canonical ensemble results

A week long MD simulation on an R10000 processor (30 ps) was compared with a week-long MC simulation (107 MC steps) using the combination of techniques described in Sec. II.1 and Sec. II.2 and the nine most-changed lipids were extracted from both. The figure below shows their comparison.

MD changes MC changes

It is clear from this figure that the MC run has produced more conformational change in these lipids than the MD run.

III.4. Grand-canonical ensemble results

The improved capability of the GCE simulations for penetrating the bilayer can be seen from the comparison of water density profiles after 107 steps of canonical and grand-canonical ensemble simulations:

and grand-canonical ensemble density profiles

As mentioned above, the significance of GCE simulations will become especially evident when the bilayer contains additional, large molecules. There is however an important unique property of GCE simulations of pure bilayers: they can detect early inadequately equilibrated systems by allowing waters into the hydrocarbon region. In the example below, a lipid molecule "slipped" into the hydrocarbon region during the regularization process. The various distributions (density profile, order parameter distribution) showed nothing untoward. However, the hydrocarbon region became "wet" during the GCE run relatively early:

DMPC molecule

Click HERE for full size picture (200 Mb)

Currently work is in progress on this system in the canonical, grand-canonical and isothermal-isobaric ensemble, with independently varying prism length and edge. The figure below shows snapshots of each system after 108 MC steps.

Snapshots in 3 ensembles

Click HERE for full size picture (210 Mb)

There is still some sign of excessive penetration of waters into the hydrocarbon region. This is in agreement with the change in simulation cell's shape that was produced by the istothermal-isobaric ensemble run: there is still "too much" space in the bilayer. It is expected that starting from the new box, the water penetration with the GCE ensemble will be normal.

IV. Conclusions

The calculations performed so far have shown that judiciously optimized Monte Carlo simulations, combined with some novel sampling techniques can successfully compete with molecular dynamics for the modeling of lipid bilayers. Furthermore, it has been shown that the cavity biased technique enables the simulation of lipid bilayers in the grand-canonical ensemble. Such simulations have the advantage of ensuring full penetration of water as well as eliminate the risk of waters remaining trapped in pockets. They can also be used as diagnostics for inadequate equilibration.


D.J. Adams, Grand canonical ensemble Monte Carlo for a Lennard-Jones fluid, Molec. Phys., 29, 311 (1975).

L.R. De Young, and K.A. Dill, Biochemistry, 27, (1988).

Duong, T-Ha, Mehler, E., and Weinstein, H. (1999) Molecular dynamics simulation of membrane bilayers and the stabilization of transmembrane helix-7 of the 5-HT2A receptor. J. Comput. Phys., 151:358-387.

S. Goldman, A simple new way to help speed up Monte Carlo convergence rates: energy-scaled displacement Monte Carlo, J. Chem. Phys., 79, 3938 (1983); Erratum: 84, 1952 (1986).

P. Huang, J.J. Perez, and G.H. Loew, Molecular dynamics simulation of phospholipid bilayers, J. Biomol. Struct. Dyn., 11, 927 (1994).

H. Heller, M. Schaefer, and K. Schulten, Molecular dynamics simulation of a bilayer of 200 lipids in the gel and in the liquid crystal phases, J. Phys. Chem., 97, 8343 (1993).

R.H. Kincaid, and H.A. Scheraga, Acceleration of Convergence in Monte Carlo Simulations of Aqueous Solutions Using the Metropolis Algorithm. Hydrophobic Hydration of Methane J. Comp. Chem., 3, 525-547 (1982)

T.J. Marrone, H. Resat, C.N. Hodge, C.-H. Chang, and J.A. McCammon, Solvation studies of DMP323 and A76928 bound to HIV protease: Analysis of water sites using grand canonical Monte Carlo simulations, Protein Science, 7, 573-579 (1998).

N.A. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, Equation of state calculation by fast computing machines, J. Chem. Phys., 21, 1087 (1953).

M. Mezei, A Cavity-Biased (T,V,mu) Monte Carlo Method for the Computer Simulation of Fluids, Molec. Phys., 40, 901-906 (1980);

M. Mezei, Grand-canonical ensemble Monte Carlo study of dense liquids: Lennard-Jones, soft spheres and water. Molec. Phys., 61, 565 (1987); erratum: 57, 1207 (1989).

M. Mezei, On the Selection of the Particle to be Perturbed in the Monte Carlo Method, J. Comp. Phys., 39, 128-136 (1981).

M. Mezei, K.A. Bencsath, S. Goldman, and S. Singh, The detailed balance energy-scaled displacement Monte Carlo algorithm, Molecular Simulation 1, 87 (1987).

J.C. Owicki, "Optimization of sampling algorithms in Monte Carlo calculations of fluids," in Computer Modeling of Matter, P.G. Lykos, ed., (American Chemical Society, Washington, D.C., 1978).

R.M. Sok, and H.J.C. Berendsen, Molecular dynamics simulation of the transport of small molecules across a polymer membrane, J. Chem. Phys. 96, 4699 (1992).

T.R. Stouch, Lipid membrane structure and dynamics studied by all-atom molecular dynamics simulations of hydrated phospholipid bilayers, Molecular Simulation 10, 335 (1993).

D.P. Tieleman, and H.J.C. Berendsen, A molecular dynamics study of the pores formed by Escherichia coli OmpF porin in a fully hydrated palmitoyloleylphospatidycholine bilayer. Biophys. J., 74, 2786-2801 (1998).

H. Resat, and M. Mezei, Grand Canonical Monte Carlo Simulation of Water Positions in Crystal Hydrates. J. Am. Chem. Soc., 116, 7451-7452 (1994).

H. Resat, and M. Mezei, Grand Canonical Ensemble Monte Carlo Simulation of the dCpG/Proflavine Crystal Hydrate. Biophys. J. 71, 1179-1190 (1996).

R.M. Venable, Y.Zhang, B.J. Hardy, R.W. Pastor, Science, 262, 223 (1993)

Back to the Mezei Lab home page