Monte Carlo techniques for simulating in lipid bilayers

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(X^N), <Q>, via generation of a set of configurations with Boltzmann probability P_B 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:

• Rotate/translate waters
• Rotate/translate whole lipids
• Change torsional states

II.1.a. Order of moves.

• The waters were moved using the preferential sampling of Owicki that selects waters near lipids with higher probability.
• Lipids for translation/rotation has been selected in the so-called shuffled cyclic order (Mezei (1981)): in full cycles, but the order in each cycle is obtained from a random permutation. It has been shown that this approach combines the benefits of cyclic (each molecule has equal chance to move) and random (avoids possible correlations) moves.
• For torsion angle change, a lipid was first selected, again in a shuffled cyclic order. However, the torsion angle itself was selected in a novel fashion using what we termed reverse cyclic order: the selection started form the end of a chain and progressed inward to the center of molecule. This order maximizes the chance that between two moves of a given torsion angle the part moved by that torsion has been changed.

In addition, torsion angle change attempt were subjected to an additional probability filter to reduce the frequency of the changes in torsion angles near the end of the chains as these torsions can change much more in a single move than the torsion angles further inward.

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 water translation/rotation stepsizes were tuned to give ca 30% acceptance rate as suggested by the results of Kincaid and Scheraga. Their conclusion was confirmed by our extensive earlier (unpublished) calculations showing that the 20-40% range is significantly better than the often quoted 50% rule. Since the convergence of the waters is much faster than that of the lipids, no further optimization was done on the water stepsizes.
• The displacement and rotation stepsizes were optimized in conjunction with the torsional stepsizes since they are closely coupled and significantly affect each other. First, short runs tuned the torsion angle stepsizes to give (nearly) uniform acceptance rate. Next, several runs were performed where the lipid rotation stepsize, relative frequency of translation/rotation and an overall factor scaling all torsion stepsizes were varied in the direction of maximum overall system change.

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.

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: Alpha_max=c/R_max^1/2.

II.3. Grand-canonical ensemble simulations

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

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

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 Å² (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 R_L was varied and the rest of the parameters were kept fixed. This gave for the extension biased sampling the following result:

In the table above, P^d_acc and P^t_acc 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:

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

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.

                        C2
                        |
C25-O23-P20-O24-C17-C14-N1-C6
|                       |
|                       C10
|
C28-O30-C31-C33-C45-C48-C51-C54-C57-C60-C63--C66--C69--C72--C75--C78
|
|
|
C36-O39-C40-C42-C82-C85-C88-C91-C94-C97-C100-C103-C106-C109-C112-C115

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.

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 (10⁷ 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.

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 10⁷ steps of canonical and grand-canonical ensemble simulations:

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:

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 10⁸ MC steps.

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.

REFERENCES

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).
Abstract

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

Back to the Mezei Lab home page