Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev Lett. Author manuscript; available in PMC 2010 August 18.
Published in final edited form as:
PMCID: PMC2923463

Simulations of Interacting Membranes in the Soft Confinement Regime


The liquid crystalline model biomembrane system consisting of a stack of interacting membranes is studied in the realistic soft confinement regime by the newly developed Fourier Monte Carlo simulation technique. In this regime experiment and simulations show that the functional form of the fluctuation pressure is more nearly exponential rather than the power law valid for the hard confinement regime. The simulations provide quantitative improvement over perturbation theory. It is shown that the harmonic theory that is routinely used to interpret x-ray scattering line shapes is valid.

Keywords: 87.10. + e, 02.70.Lq, 61.30.Cz

Stacks of lipid bilayers (see Fig. 1) are model systems for biomembranes that are much studied for two reasons. First, such stacks diffract fairly well and this facilitates determination of the structure of individual membranes, which is of primary interest in biophysics. However, these stacks are not crystals with the long-range order that is assumed in traditional biophysical analysis, but smectic liquid crystals with quasi-long-range order. Therefore, quantitative use of the scattering intensity for structure determination requires correction for the fluctuations endemic in such systems [1]. A harmonic theory [2,3] that predicts power law tails in the scattering line shapes fits membrane data very well [4,5], but the anharmonicities that are inherent in realistic potentials have remained a concern for quantitative interpretation [6,7], even though a renormalization group analysis suggested that such effects are small [8].

FIG. 1
Snapshot of a slice through a simulated stack of M = 8 two-dimensional L × L fluctuating membranes. Since internal membrane structure is irrelevant here, each membrane is depicted as a line. The average position of each membrane is shown by a ...

Stacks of bilayers are also much studied because they provide ideal environments to study fundamental interactions between bilayers, especially since the range of interbilayer distances a can be systematically varied by applying osmotic pressure P [9]. The corresponding theory [10] is an approximate first-order perturbation theory that again relies on harmonic assumptions, such as the normality of the probability distribution function of the interbilayer spacing. While this theory has been a valuable guide, it is too inaccurate to extract fundamental interbilayer interactions from P(a) data.

Both these issues are addressed using Monte Carlo simulations with realistic intermembrane potentials for the biologically relevant regime where the interbilayer water spacing a is of order 5–30 Å and each membrane is flexible with bending modulus Kc [similar, equals] 10–25kBT. In this regime, called the soft confinement regime [10], it is usually supposed that the primary interbilayer interactions for dipolar lipids are the attractive van der Waals (vdW) potential and the mysterious, but easily measured, repulsive hydration potential,


where z is the local distance between two membranes [11]. These interactions are significantly anharmonic, to the extent that the potential of a membrane midway between two neighboring membranes (at the dashed positions in Fig. 1) may have a maximum instead of a minimum. The contrasting regime, sometimes called the hard confinement regime, consists of only excluded volume or steric interactions between neighboring membranes. That regime is appropriate when a is of order 100 Å because the hydration force is short range (λ ≈ 2 Å). For hard confinement the effective interbilayer force is the entropic fluctuation pressure which decays as a−3 [12]. Fluctuation forces also exist in the soft confinement regime, but they have a different functional form which is eludicated by these simulations.

In addition to the interbilayer interactions in Eq. (1), the energy of each membrane includes a bending term, proportional to the square of the local curvature of the membrane. Let um(x, y) be the local displacement of the mth membrane from its average position as shown on Fig. 1. Periodic boundary conditions are imposed in the plane of each membrane and also along the stack, so that uM [equivalent] u0. The membranes can collide, but cannot overlap, so that um+1 + aum, where a is the average distance between membranes. The Hamiltonian of the stack is then


The simulation method, called the Fourier Monte Carlo method, was developed for single membranes between hard walls [13] and is easily extended to stacks of membranes. Each membrane in the stack is represented by a complex array of dimensions N × N of Fourier displacement amplitudes. Instead of moving one lattice site at a time, moves are made in Fourier space and a whole membrane is displaced in each move. This allows larger moves and faster equilibration, without incurring large increases in the bending energy. One difference with our previous simulations [13] is that a fixed osmotic pressure P ensemble is employed instead of the previous fixed a ensemble, so that a is obtained as a function of P rather than vice versa. Of course, use of the P ensemble is fundamentally no different, but it does have better convergence properties that we now discuss.

Simulations performed systematically as a function of lattice size, density of lattice points, and number of membranes in the stack show that accurate results for infinite, continuous membranes in infinite stacks can be obtained at one P in real time of the order of one day on a Pentium Pro PC. The most sensitive finite-size parameter is the “density” of each membrane N/L, since when N is varied from 6 to 32, the root mean square fluctuation in nearest neighbor spacing, defined as Δ, can easily change by 40%, and the changes in a are also significant. This is shown in Fig. 2, which also shows that accurate values can be obtained by extrapolation. By comparison, variations with lateral system size L at fixed N/L are negligible (≈0.2% for L ≥ 700 Å), as are variations with M (≈1% for M ≥ 8) for a and Δ.

FIG. 2
Effect of finite density (N = 6,…,32) on a (solid squares) and Δ (open circles and right hand ordinate) for realistic interaction parameters given in [14] and for P = 104 erg/cm3.

It may be noted that stacks of several membranes (M ≈ 4) have been previously considered [1517], but mostly for the critical phenomenon of unbinding; this occurs in the limit of large average membrane spacing, where the van der Waals interaction is the main one, in addition to the spatial constraints. We have also performed simulations in the hard confinement regime and obtained results for the Helfrich fluctuation free energy cflT2/Kca2 with cfl ≈ 0.1, in agreement with [15,17].

Figure 3 shows results for P and Δ as functions of a for realistic values of the potentials [5,14]. Also shown are results for a simpler case when the attractive van der Waals force is absent so that the potential experienced by each membrane has a minimum in the middle between its neighbors and is therefore more like a harmonic potential. Figure 3 also shows results based on the perturbation approximation [10] which was developed for a single membrane between hard walls. After comparing to harmonic theory for multiple membranes [5], we adjusted [10] for the case of multiple membranes by putting a factor of 4/π into the relations for the fluctuational free energy and for Δ2 and then followed essentially the same procedure as in [10]. With this factor, agreement between the perturbation theory and the simulations is quite good for small a where the total interaction is harmoniclike because the repulsive hydration potential is dominant and the fluctuations are relatively small. The simulations and the theory also agree for P(a) when there is no attractive van der Waals interaction. However, at larger a, the simulated Δ increases with a faster than the perturbation approximation for either intermembrane potential. A large difference between the simulation results and the perturbation theory occurs when the potential has an attractive van der Waals part. The theory predicts the bound value a0 = 17.9 Å for P = 0, more than 2 Å smaller than the true value a0 = 20.2 ± 0.1 Å obtained by simulation. It is also of interest to compare the values of μ [equivalent](Δ/a)2 to the hard confinement values. In Fig. 3 the range of μ is 0.06–0.12, in good agreement with experiment [5] but considerably smaller than the hard confinement estimates 0.16–0.21 [10,12,18]. Therefore, neither hard confinement theory nor harmonic perturbation theory is satisfactory for the most important case of bound bilayers at low P.

FIG. 3
Simulation results (symbols) and perturbation theory (lines) for Δ and log10 P versus a for membranes bound at P = 0 for the parameter set in [14] with attractive interaction (solid squares and lines) and for membranes unbound at P = 0 with no ...

Existing approximations [10,19] assume that the distribution of distance between two nearest neighbor membranes is essentially Gaussian. Simulation results in Fig. 4 demonstrate that the actual PDF differs substantially from a normal distribution. Although the distribution is more non-normal for small N, it is nevertheless clear that even the limiting distribution is not normal. A distinguishing feature of the true PDF is a rapid decay at small distances. This validates the use of a cutoff [11] near 2 Å to avoid the formal divergence of the vdW potential. Another feature evident in Fig. 4 is the asymmetry of the actual PDF which shows that large fluctuations to larger intermembrane spacings are more probable than large fluctuations to smaller spacings. Overall, the shape of PDF is consistent with that of the interaction potential.

FIG. 4
Probability distribution function (PDF) of the nearest neighbor distance for parameters in [14] and P = 104 erg/cm3 for different membrane densities.

We turn next to the issue of whether harmonic fluctuation theory [2,3] should be expected to be reliable when interpreting detailed line shapes [4,5] from stacks of bilayers and smectic liquid crystals that have strongly anharmonic interactions. The most important quantity for determining the line shapes for powder samples [3] is the correlation function in the z direction, which is essentially the k dependence of the mean square relative displacement of two membranes of the stack Δ2(k) [equivalent] [left angle bracket][u(m + k) — u(m)]2[right angle bracket], where m and m + k are the membrane indices, and averages over m are performed for simulation efficiency. Figure 5 shows profiles of Δ(k), obtained for stacks with various numbers M of membranes. Convergence with M suggests that values of Δ(k) are sufficiently accurate for k < M/4. However, to minimize the finite size effect in a comparison with harmonic theory, Fig. 5 compares the results of the simulation with M = 32 with the exact harmonic result, also for M = 32. In the harmonic theory the bare interbilayer interactions are approximated with a compression modulus B. In Fig. 5 a value of B = 1.9 × 1013 erg/cm4 was chosen to match the large k end of the M = 32 curve. The resulting Δ(k) profile allows one to see that Δ [equivalent] Δ(1) is in fact a good proxy for describing the long-range correlations, since the difference between Δ, implied by the “harmonic curve,” and the actual Δ for the stack is about only 0.2 Å, i.e., relatively small compared to Δ = 4.6 Å.

FIG. 5
Root mean square fluctuations Δ(k) between kth neighbor membranes for a stack with different numbers M of membranes and for the parameter set in [14] (except that L = 1400 Å) at P = 3.16 × 105 erg/cm3. Also shown is Δ( ...

How is it that the harmonic theory works quite well for the correlation functions in the preceding paragraph and not so well for P and Δ in Fig. 3? The answer is that the perturbation theory does not yield the best value of B; for the example in Fig. 5 the theory yields a larger B = 5.4 × 1013 erg/cm4 which accounts for the smaller value of Δ in Fig. 3.

We turn finally to the entropic fluctuation pressure in a stack of membranes, which is defined to be the difference between the applied pressure and the pressure due to the bare van der Waals and hydration interactions. Perturbation theory [10], experiment [5], and simulations on a single membrane between hard walls [13] all agree that the decay of the fluctuation pressure is closer to exponential, with a decay length λfl, although the value of λfl found in both experiment and the previous simulations is larger than the perturbation theory prediction λfl = 2λ. Figure 6 shows that simulations of stacks also give essentially an exponential functional form with a decay length λfl that is greater than the value 2λ predicted by perturbation theory.

FIG. 6
Simulation results for Pfl vs a for the parameter set in [14] and also for H = 0. The slope is λfl = 4.34 Å.

A long-range goal is to obtain values of the interbilayer interaction parameters; the traditional analysis [9] uses osmotic pressure P(a) data, which has recently been supplemented by fluctuation Δ(a) data [5]. One of the main results of this paper indicates that the Δ data are indeed valid, even though the analysis of the basic x-ray scattering data is based on a harmonic theory. However, the intrinsic anharmonic nature of realistic interactions between bilayers in stacks makes it difficult to devise quantitatively accurate analytic or perturbation theories. We show here that the Fourier Monte Carlo method is sufficiently fast that it provides a viable alternative. Indeed, it is now possible to consider using it as part of a comprehensive data analysis program to determine the best values of the fundamental interaction parameters in the biologically relevant soft confinement regime.


We thank Horia Petrache for many useful discussions and especially for his insight that at large distances, smaller deviations from harmonicity may be expected. This research was supported by the U.S. National Institutes of Health Grant No. GM44976.


[1] Nagle JF, et al. Biophys. J. 1996;70:1419. [PubMed]
[2] Caillé A. C. R. Acad. Sci. Paris Ser. B. 1972;891:274.
[3] Zhang R, Suter RM, Nagle JF. Phys. Rev. E. 1994;50:5047–5060. [PubMed]
[4] Zhang R, et al. Biophys. J. 1996;70:349. [PubMed]
[5] Petrache HI, Gouliaev N, Tristram-Nagle S, Zhang R, Suter RM, Nagle JF. Phys. Rev. E. 1998;57:7014.
[6] Lemmich J, et al. Phys. Rev. E. 1996;53:5169. [PubMed]
[7] Zhang R, et al. Phys. Rev. Lett. 1995;74:2832. [PubMed]
[8] Grinstein G, Pelcovits RA. Phys. Rev. Lett. 1981;47:856.
[9] Rand RP, Parsegian VA. Biochim. Biophys. Acta. 1989;988:351.
[10] Podgornik R, Parsegian VA. Langmuir. 1992;8:557.
[11] The formal divergence of the van der Waals potential at z = 0 is nonphysical, being masked by stronger steric interactions. In the simulations, a constant potential is used when |z| is smaller than where V (z) has a maximum. For these values of z, which are typically less than 2 Å, the hydration potential is so large that few membrane pairs approach this closely. We also note that the 1/z2 dependence is simplified but is quantitatively accurate for the soft confinement region where the water spacing is smaller than the membrane thickness.
[12] Helfrich W. Z. Naturforsch. 1978;33A:305.
[13] Gouliaev N, Nagle JF. Phys. Rev. E. 1998;58:881.
[14] In most parts of this Letter we use the following parameters: A = 109 erg/cm3; H = 5 × 10−14 erg; Kc = 0.5 × 10−12 erg; λ = 1.8 Å; T = 323 K; L = 700 Å; M = 8.
[15] Gompper G, Kroll DM. Europhys. Lett. 1989;9:59.
[16] Netz RR, Lipowsky R. Phys. Rev. Lett. 1993;71:3596. [PubMed]
[17] Netz RR, Lipowsky R. Europhys. Lett. 1995;29:345.
[18] Our simulations for hard confinement yield μ = 0.16.
[19] Sornette D, Ostrowsky N. J. Chem. Phys. 1986;84:4062.