Testing procedures for extracting fluctuation spectra from lipid bilayer simulations Joseph C. Albert, Lucas T. Ray, and John F. Nagle Citation: The Journal of Chemical Physics 141, 064114 (2014); doi: 10.1063/1.4892422 View online: http://dx.doi.org/10.1063/1.4892422 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Predicting solute partitioning in lipid bilayers: Free energies and partition coefficients from molecular dynamics simulations and COSMOmic J. Chem. Phys. 141, 045102 (2014); 10.1063/1.4890877 Coarse-grain simulations of active molecular machines in lipid bilayers J. Chem. Phys. 138, 195101 (2013); 10.1063/1.4803507 Thermal fluctuations in shape, thickness, and molecular orientation in lipid bilayers J. Chem. Phys. 135, 244701 (2011); 10.1063/1.3660673 Simulations of edge behavior in a mixed-lipid bilayer: Fluctuation analysis J. Chem. Phys. 126, 045105 (2007); 10.1063/1.2430714 Coarse-grained simulations of lipid bilayers J. Chem. Phys. 121, 11942 (2004); 10.1063/1.1814058
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
THE JOURNAL OF CHEMICAL PHYSICS 141, 064114 (2014)
Testing procedures for extracting fluctuation spectra from lipid bilayer simulations Joseph C. Albert,a) Lucas T. Ray, and John F. Nagleb) Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
(Received 19 May 2014; accepted 25 July 2014; published online 13 August 2014) To address concerns about how to obtain the height-height spectrum from simulations of biomembranes, we emulated the fluctuations in real space using exact input spectra. Two different methods that have given different results in the literature were then used to extract spectra from the emulated fluctuations that were then compared to the exact input spectra. A real space method shows systematic, but small deviations attributed to splines introducing an artifactual filter. A direct Fourier method obtains accurate results when the in-plane placement of the emulated particles is uncorrelated with the out-of-plane undulations, but systematic underestimates occur when the particle placement is more realistically correlated with the undulations. Although quantitative corrections cannot be estimated from our one-dimensional model, the results are qualitatively consistent with the direct Fourier method underestimating the 1/q2 spectral dependence that is characteristic of a tilt degree of freedom in simulations. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4892422] I. INTRODUCTION
Biomembranes are two-dimensional structures that undergo thermal fluctuations. Defining the unperturbed membrane to lie in the xy plane, the out-of-plane, transverse, small amplitude, undulational fluctuations are described by u(x, y), the displacements in the z direction. These fluctuations have biological consequences1, 2 and they are used to obtain quantitative structure from experimental x-ray scattering.3 Simulations of biomembranes in equilibrium are an important tool for validating the continuum theories which provide predictions of experimental Fourier spectra.4–6 Brandt et al.7 performed large scale molecular dynamics simulations from which the average undulational Fourier spectrum Su (q) ∼ |u(q)|2 was obtained employing a method8–11 called the Direct Fourier (DF) method that differs from the alternative Real Space (RS) method that has been somewhat more commonly used for lipid bilayer simulations.4, 5, 12, 13 The DF method obtained a Fourier spectrum SuDF (q) that conformed to the classical tensionless undulation spectrum14, 15 for small q, Su (q) ∼ 1/q 4 .
(1)
SuDF (q)
included an apparently additive term For larger q, equal to a suitably normalized spectrum Sρ (q) obtained from in-plane, longitudinal, fluctuations ρ(x, y) of the lipid molecules. Therefore, it was suggested that subtraction of Sρ from SuDF (q) gave the true undulation spectrum Suexact (q). The major conclusion from this procedure was that there was no crossover from the tensionless 1/q4 undulation spectrum to the 1/q2 regime that is found using the RS method. A 1/q2 regime had been predicted by earlier theories of molecular a) [email protected] b) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/141(6)/064114/8/$30.00
protrusion4, 5, 16 and is now being predicted by theories of molecular tilting.12, 13, 17 However, there is no proof of simple additivity, SuDF (q) = Suexact (q) + Sρ (q),
(2)
of the transverse and longitudinal fluctuations in the DF spectrum, so a 1/q2 term could have been suppressed by the subtraction. The DF method has been criticized for not agreeing with RS results and with tilt theories.13, 17 This is an issue that is difficult to resolve with existing molecular dynamics simulations because limitations on the size and time scales constrain the small end of the q range, thereby making it difficult to disentangle terms with different power law behavior. This compromises the ultimate goal of determining values of the bending modulus Kc and the tilt modulus Kθ from the undulation spectra Su (q). This paper addresses this issue by devising emulated data that input an exactly known spectrum Suexact (q) and then analyzing these data with the DF method to test whether subtraction of Sρ (q) from SuDF (q) returns the known Suexact (q) input. In this paper, the system is simplified to a one-dimensional analogue that we construct to have the same q dependence as the two-dimensional case; Appendix A provides theoretical background for this simplification. Our numerical results show that, when lateral (x-y) and undulation (z) fluctuations are completely uncorrelated, the DF method effectively reproduces Suexact . However, upon inclusion of a simple and necessary geometric correlation, the DF method, while continuing to give accurate results for wide ranges of the parameters, generally underestimates the undulation spectrum in the q range where a 1/q2 term would be observable. Interestingly, the RS method also suppresses the spectrum, but only for larger q compared to the DF method, so it is generally freer from artifacts that affect the interpretation of the Su (q) spectrum.
141, 064114-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
Transverse displacement u(x) (nm)
064114-2
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
using
10
N 1 u(xn )e−iqm xn , u(qm ) = N n=1
5 0
(6)
where the 1/N normalization factor is required for consistency with Eq. (5). We now discuss and quantify the b parameter in Eq. (4). Appendix A shows that the thermally averaged height-height spectrum is given by
-5 -10 0
100
200
300
400
500
600
700
800
|u(qm )|2 =
Lateral distance x (nm) FIG. 1. Typical u(x) for one set of random φ m , Kc = 20 kT, u = 0, u2 1/2 = 6.7 nm.
II. TRANSVERSE FLUCTUATIONS
We begin with a 1-dimensional continuum model of length L with periodic boundary conditions. The mth harmonic is thus related to its wavenumber qm and wavelength λm by qm = 2π/λm =
2π m . L
(3)
From Eq. (1) the amplitudes of these normal modes are postulated to be b eiφm , u(qm ) = qm2
(4)
where b is a numerical constant to be determined shortly and φ m is a random phase. Real space undulations u(x) are then obtained using u(x) =
u(qm )eiqm x =
m
b ei(qm x+φm ) , 2 q m m
(5)
where the φ m are generated as random numbers over a uniform distribution in the range [0, 2π ]. For real u(x), the sums include negative m with φ −m = −φ m . Figure 1 shows u(x) for one set of random φ m . To emulate thermal averages, many random sets are generated, each corresponding to an MD snapshot. Note that, while the phases are random, b is an exact number, so |u(q)|2 is precisely known and is the same for each emulated microstate. This eliminates an uncertainty inherent in MD simulations and accounts for our use of the word emulations instead of simulations.
III. DISCRETIZATION
As lipid membranes are comprised of discrete molecules, an MD simulation only gives u(xn ) at discrete values of xn , n = 1, . . . , N. Crucially, the lateral coordinates xn may not be evenly spaced, so the surface u(x) cannot be described by a simple discrete Fourier transform. The RS method addresses this problem by interpolating u(xn ) onto a uniform x grid which is then Fourier analyzed to obtain u(qm ). The DF method calculates the u(qm ) directly from the simulated u(xn )
kT , Kc N 2 a 2 qm4
(7)
where kT is the thermal energy, Kc is the bending modulus, and a = L/N is the mean distance between points. Following Brandt et al.,7 it is convenient to compare different values of N on the same plot by defining the normalized undulation spectrum Su (qm ) = N 2 |u(qm )|2 .
(8)
This normalization will be used both for the DF spectrum SuDF (qm ) and for the exact input spectrum Suexact (qm ). With this normalization, Suexact (qm ) = b2 /qm4 ,
(9)
where b2 =
kT , Kc a 2
(10)
and the b in Eqs. (4) and (5) equals b/N. Appropriate biophysical choices for the parameters are Kc /kT = 20 and a = 0.8 nm, giving b = 0.28 nm−1 . This choice places Su (q) on a similar scale as Fig. 3 in Brandt et al.7 Although the xn are not located on a uniform grid in a simulation of a fluid phase bilayer, it is illuminating to consider a kind of “crystalline” case for which the xn are chosen to be precisely at xn = (n − 1)a, for n from 1 to N. As is well known in solid state physics, Fourier spectra of lattices consist of Brillouin zones that repeat in q space where the boundary of the first Brillouin zone is π /a = qB > qm > −qB . Unlike the continuum case where the value of q is unlimited, for the discrete system, any q outside the first Brillouin zone is equivalent to a q in the first Brillouin zone that differs by integer multiples of 2π /a. Figure 2 shows the result of applying Eq. (6) to Eq. (5) where the independent q modes were restricted to the first Brillouin zone with amplitudes given by Eq. (4). The inset emphasizes the repeating Brillouin zones that are replicated by Eq. (6). The main figure emphasizes that Eq. (6) reproduces the input |u(qm )|2 exactly. One also sees that our real space method essentially applies a filter that underestimates the true spectrum near qB . We used a monotone piecewise cubic interpolation18 to obtain a smooth u(x) on a grid, where the result did not depend significantly on the grid size when it was smaller than about one third of the lattice spacing. The lateral density consists of delta function particles located at positions xn . The Fourier components of the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-3
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
We end this section by noting that there is no advantage to extending the input values of |u(qm )|2 for qm outside the first Brillouin zone; when many sets of φ m phases are averaged, their amplitudes are simply added to the values in the first Brillouin zone; these additions only subject SuDF (q) to negligible random noise near the zone boundaries.
6
10
4
10
-15 -10
-5
0
5
q
10
15
20
25
30
2
10
IV. LATERAL FLUCTUATIONS exact
Su 0
DF
10
Su Sρ
RS
Su
-2
10
-4
10
0.01
0.1
-1
q (nm )
1
10
FIG. 2. A “crystalline” case with particles placed laterally on a uniform grid −4 and with the input parameters in Figure 1. The exact input qm undulation modes are shown as black squares and the direct Fourier method produces the red ×’s that are periodic in successive Brillouin zones as emphasized in the inset. The in-plane structure factor Sρ (q) is a sum of delta functions located at 2 π h/a for integer h, indicated by the large blue squares. Also shown with green diamonds is a real space result SuRS obtained using a cubic spline interpolation.
longitudinal lateral density fluctuations are then ρ(qm ) =
N 1 −iq x e m n, N n
(11)
and the thermal average of |ρ(qm )|2 is the usual longitudinal structure factor. Following Brandt et al.,7 the corresponding normalized in-plane structure factor is defined to be Sρ (qm ) = N 2 u2 |ρ(qm )|2 ,
(12)
where u2 is the thermal average of the average over lattice sites of the mean square transverse displacements in real space, N 1 2 2 |u(xn )| , (13) u ≡ N n=1 which corresponds to the factor introduced by Brandt et al. and which gives Sρ the same physical dimensions and high-q limit as Su. 7 For the crystalline case, Sρ (qm ) is a delta function comb, namely, zero except when q = 2π h/a for integer h. Therefore, subtraction does not affect the exact agreement of the DF result with the exact spectrum except at these peak values. These values of q correspond to q = 0 in the first Brillouin zone for which u(q = 0) is undefined in Eq. (4). Physically, the q = 0 mode describes uniform translation in z, which is excluded in simulations by the placement of z = 0 to correspond to the center of mass in the z direction for each snapshot. The question now becomes whether the DF method is a good approximation for fluid phase simulations and over what range of q it may apply. For that, we turn in Sec. IV to adding lateral displacements to the uniformly spaced particle positions.
The main consideration for adding lateral fluctuations is how to choose the xn values. We have tried three different models. (i) Independent random displacement δ n of each xn from its gridded value na neglects any correlation caused by interactions between nearest neighbors. It can be shown that this model gives a q2 dependence to the Sρ (qm ) spectra. As this is unlike the MD simulations and as the model does not respect interactions between neighboring particles, we do not show any results for this artificial model. (ii) Assuming harmonic interactions between neighbors spaced a apart, the usual phonon spectrum for discrete particles has a spectrum that depends on q as [sin(qa/2)]−2 . The appropriate emulation input is N/2 kT cos(qm na + θm ) , (14) δn ≡ xn − na = Klat N m=1 sin(qm a/2) where the θ m are random phases and Klat is the harmonic spring constant. Appendix B raises the interesting, but not crucial, issue of discrete versus continuum input spectra. (iii) We have mainly generated the xn values using a stochastic method, reminiscent of the so-called paracrystalline model,19, 20 which, unlike model (i) but like model (ii), correlates the displacements between neighboring particles. Particle 1 is placed at x1 = 0 and then subsequent positions, n = 2, . . . , N + 1, are iteratively computed using x˜n = x˜n−1 + a + δn ,
(15)
where δ n is a random number generated over a gaussian distribution centered at 0 that has root mean square value σ lat . This creates a dependence between nearest neighbor separations, not unlike real neighboring atoms. However, this does not guarantee periodic boundary conditions because x˜N+1 is only accidentally located at (N + 1)a, so each x˜n is rescaled by the factor f =
Na . x˜N+1 − x˜1
(16)
For any of the three methods of generating the lateral fluctuations, one set of δ n generated one “snapshot” of the lateral fluctuations of discrete particles. Many random sets of δ n were similarly generated while also generating sets of random φ m for the transverse fluctuations and averages were then taken over a number W of these “snapshots.” Once the xn were chosen, we considered two options for choosing the value of u(xn ). Option 1 uses Eq. (5); this keeps the particles on the postulated, continuum u(x) surface. Option 2 sets u(xn ) equal to its value at x = na calculated on the grid; this takes the particle off the original u(x) surface, thereby involving protrusions, although these protrusions are highly correlated
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-4
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
7
1
10
u(x) ρ(x)
q
DF
-4
FIG. 3. Exaggerated schematic emphasizing that non-zero derivatives in u(x) lead to larger lateral density ρ(x).
RS
q
-1
0
V. RESULTS
Figure 4 shows the results of adding lateral fluctuations to a system that has a spectrum consisting of a sum of 1/q4 and
0.01
10
1
exact
[Su - Su
-2
3
This has the effect that fractional position of the particle xˆn along the length of the curve u(x) with total length L is equal to the fractional position of xn along the flat line of length L. Note that the periodic boundary condition xˆN+1 = L follows by summing Eq. (17) over n because xN+1 = L from Eq. (16), so this construction enforces periodic boundary conditions and it ensures that nearest-neighbor spacings are calculated tangent to the curve u(x). The input descriptors for each emulation are the number of particles N (we only show N = 1000), the average spacing between particles a (we use a = 0.8 nm, corresponding to an area/lipid of 0.64 nm2 (Ref. 21)), the bending modulus Kc , the average root mean square nearest neighbor lateral fluctuation σ lat , and the number W (always 104 ) of our uncorrelated snapshots. (If we suppose that 10 ps snapshots are reasonably uncorrelated, W = 104 would be obtained from a 100 ns simulation.) Some output descriptors are the root mean square transverse fluctuation σ trans = u2 1/2 (from Eq. (13), noting that xˆn are used when correlations are included), and the root mean square longitudinal deviation, (xn − na)2 1/2 , which is typically about 15 times greater than σ lat for N = 1000.
exact
]/Su
0
5
10
with the lateral fluctuations. We found that the differences in SuDF (q) using these two options are negligible for physically realistic values of Kc because there is little difference in the u(xn ) and the u(na) values; subsequent results are given only for option (1). There is, however, an important oversight in the above procedures for generating the xn . A reasonable expectation is that the average distance between neighbors will be the same along the u(x) curve. As shown for the exaggerated example in Fig. 3, this leads to a higher density ρ(x) for larger slopes in u(x). This introduces correlations between u(x) and ρ(x) not represented in previous methods of generating sets of {xn } and {u(xn )}, but which must exist in a physical system. This effect was achieved in our emulations by first generating a surface u(x) as described in (5) and a set of lateral positions {xn } of the methods detailed above. New lateral positions by one xˆn were then calculated iteratively using xˆ n+a L (17) α(x)dx = (xn+1 − xn ) , L xˆn where α(x) = 1 + (du/dx)2 and the total length is L α(x)dx. (18) L =
exact
[Su - Sρ- Su
Su
exact
Su
DF
0.1
exact
]/Su
q
1
Sρ
10
DF
Su - Sρ Su 0.01
RS
0.1
1
10
q FIG. 4. Results for paracrystalline lateral fluctuations and particle placement uncorrelated with undulations, averaged over W = 104 snapshots, with σ lat = 0.012 nm, to a 1-d system of N = 1000 particles with average separation a = 0.8 nm and with Kc /kT = 20, Kθ = 0.11 kT /nm2 ∼ 0.46 mN/m, and σ trans = 6.7 nm. The inset shows the relative errors of the DF and the RS methods normalized to Suexact . The exact input adds the q−2 and the q−4 terms shown by the solid lines.
1/q2 terms. With increasing q, the exact input values of Su (q) decrease and become as small as Sρ (q) at a crossover value qc ≈ 1 nm−1 near the center of a crossover region (0.7 < q < 2 nm−1 ) for SuDF ; after tracking Suexact for smaller q, SuDF tracks Sρ for larger q. The choice of N = 1000 was made to extend the small q region by a decade compared to currently feasible atomic level or coarse-grained simulations. The choice of σ lat in Figure 4 was made so that the crossover region occurs at similar values of qc as the simulations of Brandt et al.7 However, the factor u2 in Sρ scales as N2 as shown in Appendix C. Therefore, in order to obtain the same qc , it is necessary that σ lat in Eq. (12) be small. This means that the nearest neighbor peak in Sρ at q = 2π /a is much sharper than seen by Brandt et al. We also obtained a broader peak with the same qc by reducing N and increasing σ lat . In the main panel of Fig. 4 the differences SuDF − Sρ become√noisy as q increases; as expected, this noise level scales as 1/ W . Furthermore, SuDF − Sρ appears to deviate strongly from Suexact , but the logarithmic plot is very misleading in this regard because it does not show negative deviations and it spreads out the distribution for small deviations. The inset shows that the deviations (SuDF − Sρ − Suexact ) generally average to Suexact . However, in the inset we have chosen to divide these differences by Suexact because that is the quantity one wishes to extract, and that emphasizes that the statistical deviations using the DF method become rather too large for suitable averaging as q exceeds 2 nm−1 when Suexact becomes very small. The RS result, by comparison, remains smooth but exhibits systematic deviations for large q as shown in the main panel of the figure. Emulations with particle placement uncorrelated with undulations have been obtained for many different values of Kc , different size N, with and without a
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-5
Albert, Ray, and Nagle
8
J. Chem. Phys. 141, 064114 (2014) 8
3
10
2
7
10
10 DF
exact
[Su - Sρ- Su
]/Su
7
10
1
6
10
0
exact
RS
exact
[Su -Su
exact
]/Su
6
10
0
5
10
4
10
5
0.01
0.1
q
1
2
S(q)
10
1
10
10
-1
Su
DF
10
-2
Sρ
10 10
1
q
1
Su
exact
Su
DF
Sρ
0
10 DF
DF
Su - Sρ
-3
0.1
2
exact
10
0.01
3
10
Su
0
exact
]/Su
4
10
-3
3
exact
[Su - Sρ- Su
-10
10
-2
10
S(q)
DF
-1
Su - Sρ
-1
10
10
-4
Su
-2
10
RS
10 0.01
0.1
1
10
0.01
0.1
1
-1
q (nm )
q−2 tilt term and with different amounts of lateral fluctuations σ lat . The DF extracted spectrum does not show systematic deviations from the exact input. Fig. 5 shows that correlating the particle placement with the undulations has two effects. The first is in Sρ for q < 1 nm−1 . As shown in Fig. 3, projecting equally spaced points on the undulating curve leads to additional in-plane inhomogeneity and this accounts for a larger Sρ , increasingly so for the small q modes that have larger Su amplitudes. The second effect is that the DF subtraction now gives slightly smaller values than the exact input for q > 1 nm−1 , as best seen in the inset to Fig. 5. As this is in the direction of underestimating a q−2 term, we explored this effect further for other values of the system parameters. This underestimation became quite dramatic for smaller Kc /kT as shown in Fig. 6. Instead of following the input spectrum that has a q−2 tilt contribution, the DF spectrum follows the q−4 line up to q = 1 nm−1 and then plunges to negative values as emphasized by the inset to Fig. 6. In contrast, the RS method tracks the exact spectrum to larger q, roughly the same as in Fig. 4. This underestimation is qualitatively similar to what is shown in Fig. 4(a) in Brandt et al.7 There it was suggested that the negative deviations that appeared for larger q might have been due to calculating Sρ in the projected (x, y) plane. The suggested unperformed alternative was to calculate Sρ on the curve, which we will call Sρ ∗ . Calculating Sρ ∗ is essentially the same as calculating Sρ with no correlations of particle placement with undulations; although there should be a small shift in the q values because the average distance between points exceeds a, the near constancy of the uncorrelated Sρ should impart little difference compared to Sρ ∗ . Unfortunately, and somewhat surprisingly, Fig. 7 shows that this
q (nm ) FIG. 6. Results for with correlated particle placement with Kc /kT = 5 and Kθ = 3 kT /nm2 ∼ 13 mN/m that gives σ trans = 13.3 nm. Other input parameters are the same as in Fig. 5.
refinement of subtracting Sρ ∗ from the correlated SuDF gives unacceptably large systematic positive deviations from the exact input values. VI. DISCUSSION
Lipid bilayers have a non-zero thickness and there are distributions of different atoms along the bilayer normal. The phosphates are the favorite for describing the upper and lower membrane surfaces in simulations. The terminal methyls on the hydrocarbon chains are best for describing the center of the membrane. Large differences in SuDF (q) result from this choice as seen in Fig. 2(a) in Braun et al.22 Both choices gave the usual q−4 behavior in SuDF (q) for small q. But for q in the
8
10
7
10
60 50
6
10
40 30
5
10
20 4
DF
[Su - Sρ*- Su
exact
]/Su
exact
10
10
S(q)
FIG. 5. Results for particle placement correlated with undulations and paracrystalline lateral fluctuations, averaged over W = 104 snapshots, with σ lat = 0.002 nm, to a 1-d system of N = 1000 particles with average projected separation a = 0.8 nm and with Kc /kT = 20, Kθ = 12 kT /nm2 ∼ 50 mN/m and σ trans = 6.7 nm. The inset shows the differences between SuDF − Sρ and the input values Suexact normalized to Suexact . The thin solid line shows the q−4 part of the exact input.
10
-1
3
10
Su
exact
Su
DF
0 0.01
0.1
q
1
Sρ*
2
10
DF
Su - Sρ*
1
10
0
10
-1
10
-2
10
0.01
0.1
1
10
-1
q (nm ) FIG. 7. Same parameters as in Fig. 6 except that the in-plane spectrum is Sρ ∗ calculated along the undulating curves instead of projected onto the x axis.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-6
Albert, Ray, and Nagle
range 1 − 3 nm−1 , using the phosphate groups gave a strongly increasing SuDF (q), whereas using the terminal methyl groups gave a nearly constant SuDF (q). The latter behavior is similar to our emulations here, which is not surprising as our model has a single surface, similar to the terminal methyl surface in MD simulations. While it could have been of interest to study a model that has two surfaces emulating the phosphate groups in MD bilayer simulations, this paper has focussed on the reliability of methods for extracting the true spectra for which a single surface suffices. We obtained many results, such as Fig. 4, for which the DF subtraction method reproduced the exact input to high accuracy when the particles were distributed with average equal spacing on the projected x axis. However, that procedure does not correlate the particle spacing with the undulations. In the uncorrelated case, Eq. (2) indeed follows as we show in Appendix D. However, as it is more realistic to distribute the particle positions on the undulating curve with equal average spacing as indicated in Fig. 3, we performed emulations with such correlations. For the typical value Kc /kT = 20 the DF subtraction also was rather accurate, but small systematic deviations occurred as shown in Fig. 5. Fig. 6 shows that these deviations became quite large for Kc /kT = 5, indicating that the DF subtraction method is incapable of dealing quantitatively with correlations between undulations and density fluctuations, even when the alternative subtraction in Fig. 7 was tried. The underestimation of the DF method occurs in a crucial q range of the spectrum that is required to extract the effect of tilt. The underestimated DF spectrum in Fig. 6 is similar to a simulated DF spectrum, from which it was concluded that tilt has no effect.7 Although that simulation had Kc /kT = 19 and our emulations at Kc /kT = 20 do not show a similarly large effect, we do not suggest that estimates of accuracy from one-dimensional emulations should be simply transferred to two-dimensional simulations. The more significant finding in this paper is that the DF method can lead to significant underestimation of the effect of tilt on the spectrum. Figures 4 and 6 also show that a RS spline method begins to deviate from the exact input as the boundary of the first Brillouin zone is approached, consistent with RS methods introducing an effective filter as noted by Brandt et al.7 However, our results using the RS method indicate that it provides a good approximation to the spectrum up to values of q large enough to quantify either a pure q−4 spectrum or one involving a substantial q−2 additive term. This suggests that the RS method is less likely to lead to artifacts than the DF method.
ACKNOWLEDGMENTS
We acknowledge Frank Brown for persistent criticism of the DF method which motivated this study. Summer support for J.C.A. by the Undergraduate Enrichment Fund of the CMU Physics and continuing support for J.F.N. by the National Institute of General Medical Sciences of the National Institutes of Health under Award No. R01GM44976 are gratefully acknowledged.
J. Chem. Phys. 141, 064114 (2014)
APPENDIX A: COMPARISON TO 2 DIMENSIONS
Let us consider differences and similarities between the theory in one dimension vs two, beginning with a review of the two-dimensional case whose bending energy is L Kc E= ( 2 u(x, y))2 dxdy. (A1) 2 0 Straightforward Fourier transformation gives K E= c |u(qm qn )|2 (qm2 + qn2 )2 L2 2 m n =
Kc L 2 |u(q)|2 q 4 , 2 q
(A2)
and using the equipartition theorem then gives the wellknown result, |u(q)|2 =
kT 1 1 , Kc (N a)2 q 4
(A3)
where we have also used Na = L, where a is the linear spacing between points and N is the number of points in a linear direction with N2 being the total number of points. The normalized spectrum is then given by Su (qm ) in Eq. (9) with b2 given by Eq. (10). Now let us compare the one-dimensional analogue 2 L 2 d u(x) K E = c1 dx, (A4) 2 0 d 2x which immediately shows that Kc1 has the units of energy times distance, unlike Kc in Eq. (A1) which has units of energy. Again Fourier transform of u(x)
Kc1 L 2 iqm x 2 −iqn x E= −u(qm )qm e −u(qn )qn e dx 2 0 m n (A5) =
L Kc1 u(qm )u(qn )qm2 qn2 ei(qm −qn )x dx 2 m n 0
=
Kc1 L |u(qm )|2 qm4 . 2 m
(A6)
To conform to the two-dimensional case, we now set Kc1 = Kc L. Then, E=
Kc L 2 |u(qm )|2 qm4 . 2 m
(A7)
Applying the equipartition theorem to each mode gives the same result as Eq. (A3) and the subsequent renormalization gives the same Su (qm ) in Eq. (9) with b2 given by Eq. (10). An alternative way to derive the one-dimensional analogue is to assign an infinitely stiff bending modulus Kcy to curvature in the y direction of a two-dimensional membrane while retaining a finite Kcx in the x direction. Then, 2 L
2 2 1 1/2 d u(x, y) 1/2 d u(x, y) E= Kcx + Kcy dxdy, 2 0 d 2x d 2y (A8)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
Any mode with non-zero qn has infinite energy and therefore zero amplitude, so the equipartition theorem directly gives |u(qm )|2 =
1 1 kT , 2 Kcx (N a) qm4
(A10)
equivalent to Eq. (A3). These similarities support using the d = 1 model for our study.
APPENDIX B: COMPARISON OF DISCRETE AND CONTINUUM MODELS FOR FLUCTUATIONS
3
10
1
10
-1
10
0.1
-1
(B1)
−2
As Figure 9 shows, this leads to a q contribution when analyzed for power law behavior. However, the coefficient of the q−2 contribution is only 1/6 as large as the coefficient of the q−4 term, so it only exhibits noticeable devia-
Sρ Continuum Phonons DF Su
4
10
a ~ 0.8 nm
tions from q−4 behavior when q exceeds 1 nm−1 as shown in Figure 9. APPENDIX C: MEAN SQUARE DEVIATIONS IN d = 1 and d = 2
The mean square deviation u2 plays a role in comparing Sρ (q) with Su (q). It is useful to know how this quantity depends upon the size of the system, in both d = 1 and d = 2 dimensions. In one dimension, using Eq. (5) and b = b/N 1 L u2 ≡ |u(x)|2 dx L 0 N/2 b2 1 2b2 N 2 a 4 1 = ∼ N 2. 2 q4 4 4 N (2π ) m m m m=1
(C1)
The final sum depends only weakly upon N, going from 1 at N = 2 to ζ (4) ≈ 1.0823, so u2 scales as N2 . In two dimensions L 1 2 u ≡ |u(x, y)|2 dxdy (L)2 0
Paracrystalline DF Su
10
qB=π/a
FIG. 9. Comparison of undulation spectra when a continuum model is used (solid black line) versus a discrete model (blue squares). The difference (open green circles) is well represented over the range available to simulations by a q−2 term (dashed red line).
=
6
1
q (nm )
Figure 8 addresses the generation of lateral fluctuations with the paracrystalline method versus the phonon method. When Eq. (14) was used for the phonon method, the results were too close to those for the paracrystalline method to be visually distinguishable. Instead, the comparison in Figure 8 is to the continuum phonon case which replaces sin(qa/2) in Eq. (14) with the long wavelength limiting qa/2 behavior. Using the continuum phonon spectrum gives the artifact of a discontinuous derivative at the Brillouin zone boundary compared to the discrete phonon model and the paracrystalline model. The comparison of discrete and continuum phonon generation raises a similar issue regarding transverse fluctuations. If one uses the discrete version for the bending energy, then Suexact (q) = b2 (a/2sin(qa/2))4 .
-4
continuum q discrete discrete - continuum -2 q
5
10 2
and the Fourier transformation gives 1 1/2 1/2 |u(qm qn )|2 (Kcx qm2 + Kcy qn2 )2 (L)2 . (A9) E= 2 m n
undulation spectrum |u(q)|
064114-7
Sρ
=
2
10
N/2 1 4b2 N 2 a 4 ∼ N 2, (2π )4 n,m=1 (n2 + m2 )2
(C2)
so the mean square amplitude diverges as the square of the linear size N and L for both the d = 1 and d = 2 systems.
0
10
-2
10
0.01
0.1
-1
q (nm )
1
10
FIG. 8. Comparison of results obtained with paracrystalline generation and discrete phonon generation of lateral fluctuations are indistinguishable in this figure. The phonon result shown is for phonon generation using a continuum model. Input parameters are the same as in Figure 4 except that σ lat = 0.0018 nm for the phonon case.
APPENDIX D: DERIVATION OF CONDITIONS FOR VALIDITY OF THE DIRECT FOURIER METHOD
In this appendix, we first show that the direct Fourier spectrum SuDF (q) is equal to the true fluctuation spectrum uexact (q) convolved with the lateral density structure factor ρ(q) and multiplied by a constant. We further show that in relevant limits this expression reduces to the additive relation
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-8
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
(2) posited by Brandt et al. when ρ(q) is uncorrelated with uexact (q). SuDF is defined by 2 N DF 2 1 −iqxn u(xn )e Su (q) ≡ N (D1) . N n=1 This can be written as an integral of the summand times N delta functions
2 N ∞ u(x)e−iqx δ(x − xn ) dx . SuDF (q) = N 2 x=0 n=1 (D2) By the convolution theorem, the Fourier transform of this product is the convolution of the individual transforms of u(x) and the sum of delta functions, denoted w(x). Consequently, 2 ∞ DF 2 u(q )w(q − q )dq , (D3) Su (q) = N −∞ where w(q) is w(q) =
L
x=0
1 N
N
Separating the terms with ρ(0) factors, noting that ρ(0) = 1, and recalling that |u(q)| = |u( − q)|, ⎛ 2|u(qp )|2 SuDF (qm ) = N 2 ⎝|u(qm )|2 + qp >0,qp =qm
×
(|ρ(qm + qp )|2 + |ρ(qm − qp )|2 ) 2
. (D7)
In q regions where ρ does not have significant curvature and can be approximated as the mean of ρ(q ± ) for small , one obtains ⎞ ⎛ 2|u(qp )|2 ⎠ SuDF (qm ) ≈ N 2 ⎝|u(qm )|2 + |ρ(qm )|2 qp >0 =q
≈ Suexact (qm )(1 − |ρ(qm )|2 ) + Sρ (qm ).
(D8)
In our emulations and likely for actual simulations, |ρ(q)|2 is small compared to unity, so this reduces to Eq. (2).
δ(x − xn ) e−iqx dx
1 Y.
n=1
N 1 −iqx n = ρ(q). = e N n=1
(D4)
The ρ(q) in Eq. (D4) is exactly the lateral density fluctuation components in Eq. (11). Thus, the direct Fourier spectrum SuDF in Eq. (D1) is equal to the square of the convolution of the true undulation spectrum u(q) with the lateral density spectrum ρ(q), multiplied by a factor of N2 . The analysis so far has been identical to that done for non-equispaced time series by Ref. 23. Unfortunately, as noted in Ref. 23, deconvolving Eq. (D3) to get the u is not always feasible when working with actual, noisy data. We next demonstrate that Eq. (D3) reduces to Eq. (2) proposed in Ref. 7, if lateral fluctuations are sufficiently small and statistically uncorrelated with undulations. Crucially, if either of these conditions is violated, the additive relation Eq. (2) will not necessarily hold. Substituting Eq. (D4) into (D3) and using the assumption of periodic boundary conditions on u(q) gives 2 DF 2 u(qp )ρ(qm − qp ) . (D5) Su (qm ) = N qp Assuming the u and ρ modes are uncorrelated, then SuDF (qm ) = N 2 |u(qp )|2 |ρ(qm − qp )|2 .
(D6)
Shibata, J. J. Hu, M. M. Kozlov, and T. A. Rapoport, Annu. Rev. Cell Dev. Biol. 25, 329 (2009). 2 L. V. Chernomordik and M. M. Kozlov, Annu. Rev. Biochem. 72, 175 (2003). 3 Y. Liu and J. F. Nagle, Phys. Rev. E 69, 040901 (2004). 4 E. Lindahl and O. Edholm, Biophys. J. 79, 426 (2000). 5 R. Goetz, G. Gompper, and R. Lipowsky, Phys. Rev. Lett. 82, 221 (1999). 6 J. F. Nagle, M. S. Jablin, S. Tristram-Nagle, and K. Akabori, “What are the true values of the bending modulus of simple lipid bilayers?,” Chem. Phys. Lipids (in press). 7 E. G. Brandt, A. R. Braun, J. N. Sachs, J. F. Nagle, and O. Edholm, Biophys. J. 100, 2104 (2011). 8 W. K. den Otter and W. J. Briels, J. Chem. Phys. 118, 4712 (2003). 9 J. Stecki, J. Chem. Phys. 120, 3508 (2004). 10 J. Stecki, Adv. Chem. Phys. 144, 157 (2010). 11 J. Stecki, J. Chem. Phys. 137, 116102 (2012). 12 E. R. May, A. Narang, and D. I. Kopelevich, Phys. Rev. E 76, 021913 (2007). 13 M. C. Watson, E. G. Brandt, A. J. Welch, and F. L. H. Brown, Phys. Rev. Lett. 109, 028102 (2012). 14 W. Helfrich, Z. Naturforsch., C 28, 693 (1973). 15 W. Helfrich and R. M. Servuss, Il Nuovo Cimento D 3, 137–151 (1984). 16 G. Brannigan and F. Brown, Biophys. J. 90, 1501 (2006). 17 M. C. Watson, E. S. Penev, P. M. Welch, and F. L. H. Brown, J. Chem. Phys. 135, 244701 (2011). 18 F. N. Fritsch and R. E. Carlson, SIAM J. Numer. Anal. 17, 238 (1980). 19 R. Hosemann and S. N. Bagchi, Direct Analysis of Diffraction by Matter (North-Holland, Amsterdam, 1962). 20 A. Guinier, X-Ray Diffraction (W. H. Freeman, San Francisco, 1963). 21 J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. Acta, Rev. Biomembr. 1469, 159 (2000). 22 A. R. Braun, E. G. Brandt, O. Edholm, J. F. Nagle, and J. N. Sachs, Biophys. J. 100, 2112 (2011). 23 T. J. Deeming, Astrophys. Space Sci. 36, 137 (1975).
qp
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
THE JOURNAL OF CHEMICAL PHYSICS 141, 064114 (2014)
Testing procedures for extracting fluctuation spectra from lipid bilayer simulations Joseph C. Albert,a) Lucas T. Ray, and John F. Nagleb) Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
(Received 19 May 2014; accepted 25 July 2014; published online 13 August 2014) To address concerns about how to obtain the height-height spectrum from simulations of biomembranes, we emulated the fluctuations in real space using exact input spectra. Two different methods that have given different results in the literature were then used to extract spectra from the emulated fluctuations that were then compared to the exact input spectra. A real space method shows systematic, but small deviations attributed to splines introducing an artifactual filter. A direct Fourier method obtains accurate results when the in-plane placement of the emulated particles is uncorrelated with the out-of-plane undulations, but systematic underestimates occur when the particle placement is more realistically correlated with the undulations. Although quantitative corrections cannot be estimated from our one-dimensional model, the results are qualitatively consistent with the direct Fourier method underestimating the 1/q2 spectral dependence that is characteristic of a tilt degree of freedom in simulations. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4892422] I. INTRODUCTION
Biomembranes are two-dimensional structures that undergo thermal fluctuations. Defining the unperturbed membrane to lie in the xy plane, the out-of-plane, transverse, small amplitude, undulational fluctuations are described by u(x, y), the displacements in the z direction. These fluctuations have biological consequences1, 2 and they are used to obtain quantitative structure from experimental x-ray scattering.3 Simulations of biomembranes in equilibrium are an important tool for validating the continuum theories which provide predictions of experimental Fourier spectra.4–6 Brandt et al.7 performed large scale molecular dynamics simulations from which the average undulational Fourier spectrum Su (q) ∼ |u(q)|2 was obtained employing a method8–11 called the Direct Fourier (DF) method that differs from the alternative Real Space (RS) method that has been somewhat more commonly used for lipid bilayer simulations.4, 5, 12, 13 The DF method obtained a Fourier spectrum SuDF (q) that conformed to the classical tensionless undulation spectrum14, 15 for small q, Su (q) ∼ 1/q 4 .
(1)
SuDF (q)
included an apparently additive term For larger q, equal to a suitably normalized spectrum Sρ (q) obtained from in-plane, longitudinal, fluctuations ρ(x, y) of the lipid molecules. Therefore, it was suggested that subtraction of Sρ from SuDF (q) gave the true undulation spectrum Suexact (q). The major conclusion from this procedure was that there was no crossover from the tensionless 1/q4 undulation spectrum to the 1/q2 regime that is found using the RS method. A 1/q2 regime had been predicted by earlier theories of molecular a) [email protected] b) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/141(6)/064114/8/$30.00
protrusion4, 5, 16 and is now being predicted by theories of molecular tilting.12, 13, 17 However, there is no proof of simple additivity, SuDF (q) = Suexact (q) + Sρ (q),
(2)
of the transverse and longitudinal fluctuations in the DF spectrum, so a 1/q2 term could have been suppressed by the subtraction. The DF method has been criticized for not agreeing with RS results and with tilt theories.13, 17 This is an issue that is difficult to resolve with existing molecular dynamics simulations because limitations on the size and time scales constrain the small end of the q range, thereby making it difficult to disentangle terms with different power law behavior. This compromises the ultimate goal of determining values of the bending modulus Kc and the tilt modulus Kθ from the undulation spectra Su (q). This paper addresses this issue by devising emulated data that input an exactly known spectrum Suexact (q) and then analyzing these data with the DF method to test whether subtraction of Sρ (q) from SuDF (q) returns the known Suexact (q) input. In this paper, the system is simplified to a one-dimensional analogue that we construct to have the same q dependence as the two-dimensional case; Appendix A provides theoretical background for this simplification. Our numerical results show that, when lateral (x-y) and undulation (z) fluctuations are completely uncorrelated, the DF method effectively reproduces Suexact . However, upon inclusion of a simple and necessary geometric correlation, the DF method, while continuing to give accurate results for wide ranges of the parameters, generally underestimates the undulation spectrum in the q range where a 1/q2 term would be observable. Interestingly, the RS method also suppresses the spectrum, but only for larger q compared to the DF method, so it is generally freer from artifacts that affect the interpretation of the Su (q) spectrum.
141, 064114-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
Transverse displacement u(x) (nm)
064114-2
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
using
10
N 1 u(xn )e−iqm xn , u(qm ) = N n=1
5 0
(6)
where the 1/N normalization factor is required for consistency with Eq. (5). We now discuss and quantify the b parameter in Eq. (4). Appendix A shows that the thermally averaged height-height spectrum is given by
-5 -10 0
100
200
300
400
500
600
700
800
|u(qm )|2 =
Lateral distance x (nm) FIG. 1. Typical u(x) for one set of random φ m , Kc = 20 kT, u = 0, u2 1/2 = 6.7 nm.
II. TRANSVERSE FLUCTUATIONS
We begin with a 1-dimensional continuum model of length L with periodic boundary conditions. The mth harmonic is thus related to its wavenumber qm and wavelength λm by qm = 2π/λm =
2π m . L
(3)
From Eq. (1) the amplitudes of these normal modes are postulated to be b eiφm , u(qm ) = qm2
(4)
where b is a numerical constant to be determined shortly and φ m is a random phase. Real space undulations u(x) are then obtained using u(x) =
u(qm )eiqm x =
m
b ei(qm x+φm ) , 2 q m m
(5)
where the φ m are generated as random numbers over a uniform distribution in the range [0, 2π ]. For real u(x), the sums include negative m with φ −m = −φ m . Figure 1 shows u(x) for one set of random φ m . To emulate thermal averages, many random sets are generated, each corresponding to an MD snapshot. Note that, while the phases are random, b is an exact number, so |u(q)|2 is precisely known and is the same for each emulated microstate. This eliminates an uncertainty inherent in MD simulations and accounts for our use of the word emulations instead of simulations.
III. DISCRETIZATION
As lipid membranes are comprised of discrete molecules, an MD simulation only gives u(xn ) at discrete values of xn , n = 1, . . . , N. Crucially, the lateral coordinates xn may not be evenly spaced, so the surface u(x) cannot be described by a simple discrete Fourier transform. The RS method addresses this problem by interpolating u(xn ) onto a uniform x grid which is then Fourier analyzed to obtain u(qm ). The DF method calculates the u(qm ) directly from the simulated u(xn )
kT , Kc N 2 a 2 qm4
(7)
where kT is the thermal energy, Kc is the bending modulus, and a = L/N is the mean distance between points. Following Brandt et al.,7 it is convenient to compare different values of N on the same plot by defining the normalized undulation spectrum Su (qm ) = N 2 |u(qm )|2 .
(8)
This normalization will be used both for the DF spectrum SuDF (qm ) and for the exact input spectrum Suexact (qm ). With this normalization, Suexact (qm ) = b2 /qm4 ,
(9)
where b2 =
kT , Kc a 2
(10)
and the b in Eqs. (4) and (5) equals b/N. Appropriate biophysical choices for the parameters are Kc /kT = 20 and a = 0.8 nm, giving b = 0.28 nm−1 . This choice places Su (q) on a similar scale as Fig. 3 in Brandt et al.7 Although the xn are not located on a uniform grid in a simulation of a fluid phase bilayer, it is illuminating to consider a kind of “crystalline” case for which the xn are chosen to be precisely at xn = (n − 1)a, for n from 1 to N. As is well known in solid state physics, Fourier spectra of lattices consist of Brillouin zones that repeat in q space where the boundary of the first Brillouin zone is π /a = qB > qm > −qB . Unlike the continuum case where the value of q is unlimited, for the discrete system, any q outside the first Brillouin zone is equivalent to a q in the first Brillouin zone that differs by integer multiples of 2π /a. Figure 2 shows the result of applying Eq. (6) to Eq. (5) where the independent q modes were restricted to the first Brillouin zone with amplitudes given by Eq. (4). The inset emphasizes the repeating Brillouin zones that are replicated by Eq. (6). The main figure emphasizes that Eq. (6) reproduces the input |u(qm )|2 exactly. One also sees that our real space method essentially applies a filter that underestimates the true spectrum near qB . We used a monotone piecewise cubic interpolation18 to obtain a smooth u(x) on a grid, where the result did not depend significantly on the grid size when it was smaller than about one third of the lattice spacing. The lateral density consists of delta function particles located at positions xn . The Fourier components of the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-3
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
We end this section by noting that there is no advantage to extending the input values of |u(qm )|2 for qm outside the first Brillouin zone; when many sets of φ m phases are averaged, their amplitudes are simply added to the values in the first Brillouin zone; these additions only subject SuDF (q) to negligible random noise near the zone boundaries.
6
10
4
10
-15 -10
-5
0
5
q
10
15
20
25
30
2
10
IV. LATERAL FLUCTUATIONS exact
Su 0
DF
10
Su Sρ
RS
Su
-2
10
-4
10
0.01
0.1
-1
q (nm )
1
10
FIG. 2. A “crystalline” case with particles placed laterally on a uniform grid −4 and with the input parameters in Figure 1. The exact input qm undulation modes are shown as black squares and the direct Fourier method produces the red ×’s that are periodic in successive Brillouin zones as emphasized in the inset. The in-plane structure factor Sρ (q) is a sum of delta functions located at 2 π h/a for integer h, indicated by the large blue squares. Also shown with green diamonds is a real space result SuRS obtained using a cubic spline interpolation.
longitudinal lateral density fluctuations are then ρ(qm ) =
N 1 −iq x e m n, N n
(11)
and the thermal average of |ρ(qm )|2 is the usual longitudinal structure factor. Following Brandt et al.,7 the corresponding normalized in-plane structure factor is defined to be Sρ (qm ) = N 2 u2 |ρ(qm )|2 ,
(12)
where u2 is the thermal average of the average over lattice sites of the mean square transverse displacements in real space, N 1 2 2 |u(xn )| , (13) u ≡ N n=1 which corresponds to the factor introduced by Brandt et al. and which gives Sρ the same physical dimensions and high-q limit as Su. 7 For the crystalline case, Sρ (qm ) is a delta function comb, namely, zero except when q = 2π h/a for integer h. Therefore, subtraction does not affect the exact agreement of the DF result with the exact spectrum except at these peak values. These values of q correspond to q = 0 in the first Brillouin zone for which u(q = 0) is undefined in Eq. (4). Physically, the q = 0 mode describes uniform translation in z, which is excluded in simulations by the placement of z = 0 to correspond to the center of mass in the z direction for each snapshot. The question now becomes whether the DF method is a good approximation for fluid phase simulations and over what range of q it may apply. For that, we turn in Sec. IV to adding lateral displacements to the uniformly spaced particle positions.
The main consideration for adding lateral fluctuations is how to choose the xn values. We have tried three different models. (i) Independent random displacement δ n of each xn from its gridded value na neglects any correlation caused by interactions between nearest neighbors. It can be shown that this model gives a q2 dependence to the Sρ (qm ) spectra. As this is unlike the MD simulations and as the model does not respect interactions between neighboring particles, we do not show any results for this artificial model. (ii) Assuming harmonic interactions between neighbors spaced a apart, the usual phonon spectrum for discrete particles has a spectrum that depends on q as [sin(qa/2)]−2 . The appropriate emulation input is N/2 kT cos(qm na + θm ) , (14) δn ≡ xn − na = Klat N m=1 sin(qm a/2) where the θ m are random phases and Klat is the harmonic spring constant. Appendix B raises the interesting, but not crucial, issue of discrete versus continuum input spectra. (iii) We have mainly generated the xn values using a stochastic method, reminiscent of the so-called paracrystalline model,19, 20 which, unlike model (i) but like model (ii), correlates the displacements between neighboring particles. Particle 1 is placed at x1 = 0 and then subsequent positions, n = 2, . . . , N + 1, are iteratively computed using x˜n = x˜n−1 + a + δn ,
(15)
where δ n is a random number generated over a gaussian distribution centered at 0 that has root mean square value σ lat . This creates a dependence between nearest neighbor separations, not unlike real neighboring atoms. However, this does not guarantee periodic boundary conditions because x˜N+1 is only accidentally located at (N + 1)a, so each x˜n is rescaled by the factor f =
Na . x˜N+1 − x˜1
(16)
For any of the three methods of generating the lateral fluctuations, one set of δ n generated one “snapshot” of the lateral fluctuations of discrete particles. Many random sets of δ n were similarly generated while also generating sets of random φ m for the transverse fluctuations and averages were then taken over a number W of these “snapshots.” Once the xn were chosen, we considered two options for choosing the value of u(xn ). Option 1 uses Eq. (5); this keeps the particles on the postulated, continuum u(x) surface. Option 2 sets u(xn ) equal to its value at x = na calculated on the grid; this takes the particle off the original u(x) surface, thereby involving protrusions, although these protrusions are highly correlated
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-4
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
7
1
10
u(x) ρ(x)
q
DF
-4
FIG. 3. Exaggerated schematic emphasizing that non-zero derivatives in u(x) lead to larger lateral density ρ(x).
RS
q
-1
0
V. RESULTS
Figure 4 shows the results of adding lateral fluctuations to a system that has a spectrum consisting of a sum of 1/q4 and
0.01
10
1
exact
[Su - Su
-2
3
This has the effect that fractional position of the particle xˆn along the length of the curve u(x) with total length L is equal to the fractional position of xn along the flat line of length L. Note that the periodic boundary condition xˆN+1 = L follows by summing Eq. (17) over n because xN+1 = L from Eq. (16), so this construction enforces periodic boundary conditions and it ensures that nearest-neighbor spacings are calculated tangent to the curve u(x). The input descriptors for each emulation are the number of particles N (we only show N = 1000), the average spacing between particles a (we use a = 0.8 nm, corresponding to an area/lipid of 0.64 nm2 (Ref. 21)), the bending modulus Kc , the average root mean square nearest neighbor lateral fluctuation σ lat , and the number W (always 104 ) of our uncorrelated snapshots. (If we suppose that 10 ps snapshots are reasonably uncorrelated, W = 104 would be obtained from a 100 ns simulation.) Some output descriptors are the root mean square transverse fluctuation σ trans = u2 1/2 (from Eq. (13), noting that xˆn are used when correlations are included), and the root mean square longitudinal deviation, (xn − na)2 1/2 , which is typically about 15 times greater than σ lat for N = 1000.
exact
]/Su
0
5
10
with the lateral fluctuations. We found that the differences in SuDF (q) using these two options are negligible for physically realistic values of Kc because there is little difference in the u(xn ) and the u(na) values; subsequent results are given only for option (1). There is, however, an important oversight in the above procedures for generating the xn . A reasonable expectation is that the average distance between neighbors will be the same along the u(x) curve. As shown for the exaggerated example in Fig. 3, this leads to a higher density ρ(x) for larger slopes in u(x). This introduces correlations between u(x) and ρ(x) not represented in previous methods of generating sets of {xn } and {u(xn )}, but which must exist in a physical system. This effect was achieved in our emulations by first generating a surface u(x) as described in (5) and a set of lateral positions {xn } of the methods detailed above. New lateral positions by one xˆn were then calculated iteratively using xˆ n+a L (17) α(x)dx = (xn+1 − xn ) , L xˆn where α(x) = 1 + (du/dx)2 and the total length is L α(x)dx. (18) L =
exact
[Su - Sρ- Su
Su
exact
Su
DF
0.1
exact
]/Su
q
1
Sρ
10
DF
Su - Sρ Su 0.01
RS
0.1
1
10
q FIG. 4. Results for paracrystalline lateral fluctuations and particle placement uncorrelated with undulations, averaged over W = 104 snapshots, with σ lat = 0.012 nm, to a 1-d system of N = 1000 particles with average separation a = 0.8 nm and with Kc /kT = 20, Kθ = 0.11 kT /nm2 ∼ 0.46 mN/m, and σ trans = 6.7 nm. The inset shows the relative errors of the DF and the RS methods normalized to Suexact . The exact input adds the q−2 and the q−4 terms shown by the solid lines.
1/q2 terms. With increasing q, the exact input values of Su (q) decrease and become as small as Sρ (q) at a crossover value qc ≈ 1 nm−1 near the center of a crossover region (0.7 < q < 2 nm−1 ) for SuDF ; after tracking Suexact for smaller q, SuDF tracks Sρ for larger q. The choice of N = 1000 was made to extend the small q region by a decade compared to currently feasible atomic level or coarse-grained simulations. The choice of σ lat in Figure 4 was made so that the crossover region occurs at similar values of qc as the simulations of Brandt et al.7 However, the factor u2 in Sρ scales as N2 as shown in Appendix C. Therefore, in order to obtain the same qc , it is necessary that σ lat in Eq. (12) be small. This means that the nearest neighbor peak in Sρ at q = 2π /a is much sharper than seen by Brandt et al. We also obtained a broader peak with the same qc by reducing N and increasing σ lat . In the main panel of Fig. 4 the differences SuDF − Sρ become√noisy as q increases; as expected, this noise level scales as 1/ W . Furthermore, SuDF − Sρ appears to deviate strongly from Suexact , but the logarithmic plot is very misleading in this regard because it does not show negative deviations and it spreads out the distribution for small deviations. The inset shows that the deviations (SuDF − Sρ − Suexact ) generally average to Suexact . However, in the inset we have chosen to divide these differences by Suexact because that is the quantity one wishes to extract, and that emphasizes that the statistical deviations using the DF method become rather too large for suitable averaging as q exceeds 2 nm−1 when Suexact becomes very small. The RS result, by comparison, remains smooth but exhibits systematic deviations for large q as shown in the main panel of the figure. Emulations with particle placement uncorrelated with undulations have been obtained for many different values of Kc , different size N, with and without a
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-5
Albert, Ray, and Nagle
8
J. Chem. Phys. 141, 064114 (2014) 8
3
10
2
7
10
10 DF
exact
[Su - Sρ- Su
]/Su
7
10
1
6
10
0
exact
RS
exact
[Su -Su
exact
]/Su
6
10
0
5
10
4
10
5
0.01
0.1
q
1
2
S(q)
10
1
10
10
-1
Su
DF
10
-2
Sρ
10 10
1
q
1
Su
exact
Su
DF
Sρ
0
10 DF
DF
Su - Sρ
-3
0.1
2
exact
10
0.01
3
10
Su
0
exact
]/Su
4
10
-3
3
exact
[Su - Sρ- Su
-10
10
-2
10
S(q)
DF
-1
Su - Sρ
-1
10
10
-4
Su
-2
10
RS
10 0.01
0.1
1
10
0.01
0.1
1
-1
q (nm )
q−2 tilt term and with different amounts of lateral fluctuations σ lat . The DF extracted spectrum does not show systematic deviations from the exact input. Fig. 5 shows that correlating the particle placement with the undulations has two effects. The first is in Sρ for q < 1 nm−1 . As shown in Fig. 3, projecting equally spaced points on the undulating curve leads to additional in-plane inhomogeneity and this accounts for a larger Sρ , increasingly so for the small q modes that have larger Su amplitudes. The second effect is that the DF subtraction now gives slightly smaller values than the exact input for q > 1 nm−1 , as best seen in the inset to Fig. 5. As this is in the direction of underestimating a q−2 term, we explored this effect further for other values of the system parameters. This underestimation became quite dramatic for smaller Kc /kT as shown in Fig. 6. Instead of following the input spectrum that has a q−2 tilt contribution, the DF spectrum follows the q−4 line up to q = 1 nm−1 and then plunges to negative values as emphasized by the inset to Fig. 6. In contrast, the RS method tracks the exact spectrum to larger q, roughly the same as in Fig. 4. This underestimation is qualitatively similar to what is shown in Fig. 4(a) in Brandt et al.7 There it was suggested that the negative deviations that appeared for larger q might have been due to calculating Sρ in the projected (x, y) plane. The suggested unperformed alternative was to calculate Sρ on the curve, which we will call Sρ ∗ . Calculating Sρ ∗ is essentially the same as calculating Sρ with no correlations of particle placement with undulations; although there should be a small shift in the q values because the average distance between points exceeds a, the near constancy of the uncorrelated Sρ should impart little difference compared to Sρ ∗ . Unfortunately, and somewhat surprisingly, Fig. 7 shows that this
q (nm ) FIG. 6. Results for with correlated particle placement with Kc /kT = 5 and Kθ = 3 kT /nm2 ∼ 13 mN/m that gives σ trans = 13.3 nm. Other input parameters are the same as in Fig. 5.
refinement of subtracting Sρ ∗ from the correlated SuDF gives unacceptably large systematic positive deviations from the exact input values. VI. DISCUSSION
Lipid bilayers have a non-zero thickness and there are distributions of different atoms along the bilayer normal. The phosphates are the favorite for describing the upper and lower membrane surfaces in simulations. The terminal methyls on the hydrocarbon chains are best for describing the center of the membrane. Large differences in SuDF (q) result from this choice as seen in Fig. 2(a) in Braun et al.22 Both choices gave the usual q−4 behavior in SuDF (q) for small q. But for q in the
8
10
7
10
60 50
6
10
40 30
5
10
20 4
DF
[Su - Sρ*- Su
exact
]/Su
exact
10
10
S(q)
FIG. 5. Results for particle placement correlated with undulations and paracrystalline lateral fluctuations, averaged over W = 104 snapshots, with σ lat = 0.002 nm, to a 1-d system of N = 1000 particles with average projected separation a = 0.8 nm and with Kc /kT = 20, Kθ = 12 kT /nm2 ∼ 50 mN/m and σ trans = 6.7 nm. The inset shows the differences between SuDF − Sρ and the input values Suexact normalized to Suexact . The thin solid line shows the q−4 part of the exact input.
10
-1
3
10
Su
exact
Su
DF
0 0.01
0.1
q
1
Sρ*
2
10
DF
Su - Sρ*
1
10
0
10
-1
10
-2
10
0.01
0.1
1
10
-1
q (nm ) FIG. 7. Same parameters as in Fig. 6 except that the in-plane spectrum is Sρ ∗ calculated along the undulating curves instead of projected onto the x axis.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-6
Albert, Ray, and Nagle
range 1 − 3 nm−1 , using the phosphate groups gave a strongly increasing SuDF (q), whereas using the terminal methyl groups gave a nearly constant SuDF (q). The latter behavior is similar to our emulations here, which is not surprising as our model has a single surface, similar to the terminal methyl surface in MD simulations. While it could have been of interest to study a model that has two surfaces emulating the phosphate groups in MD bilayer simulations, this paper has focussed on the reliability of methods for extracting the true spectra for which a single surface suffices. We obtained many results, such as Fig. 4, for which the DF subtraction method reproduced the exact input to high accuracy when the particles were distributed with average equal spacing on the projected x axis. However, that procedure does not correlate the particle spacing with the undulations. In the uncorrelated case, Eq. (2) indeed follows as we show in Appendix D. However, as it is more realistic to distribute the particle positions on the undulating curve with equal average spacing as indicated in Fig. 3, we performed emulations with such correlations. For the typical value Kc /kT = 20 the DF subtraction also was rather accurate, but small systematic deviations occurred as shown in Fig. 5. Fig. 6 shows that these deviations became quite large for Kc /kT = 5, indicating that the DF subtraction method is incapable of dealing quantitatively with correlations between undulations and density fluctuations, even when the alternative subtraction in Fig. 7 was tried. The underestimation of the DF method occurs in a crucial q range of the spectrum that is required to extract the effect of tilt. The underestimated DF spectrum in Fig. 6 is similar to a simulated DF spectrum, from which it was concluded that tilt has no effect.7 Although that simulation had Kc /kT = 19 and our emulations at Kc /kT = 20 do not show a similarly large effect, we do not suggest that estimates of accuracy from one-dimensional emulations should be simply transferred to two-dimensional simulations. The more significant finding in this paper is that the DF method can lead to significant underestimation of the effect of tilt on the spectrum. Figures 4 and 6 also show that a RS spline method begins to deviate from the exact input as the boundary of the first Brillouin zone is approached, consistent with RS methods introducing an effective filter as noted by Brandt et al.7 However, our results using the RS method indicate that it provides a good approximation to the spectrum up to values of q large enough to quantify either a pure q−4 spectrum or one involving a substantial q−2 additive term. This suggests that the RS method is less likely to lead to artifacts than the DF method.
ACKNOWLEDGMENTS
We acknowledge Frank Brown for persistent criticism of the DF method which motivated this study. Summer support for J.C.A. by the Undergraduate Enrichment Fund of the CMU Physics and continuing support for J.F.N. by the National Institute of General Medical Sciences of the National Institutes of Health under Award No. R01GM44976 are gratefully acknowledged.
J. Chem. Phys. 141, 064114 (2014)
APPENDIX A: COMPARISON TO 2 DIMENSIONS
Let us consider differences and similarities between the theory in one dimension vs two, beginning with a review of the two-dimensional case whose bending energy is L Kc E= ( 2 u(x, y))2 dxdy. (A1) 2 0 Straightforward Fourier transformation gives K E= c |u(qm qn )|2 (qm2 + qn2 )2 L2 2 m n =
Kc L 2 |u(q)|2 q 4 , 2 q
(A2)
and using the equipartition theorem then gives the wellknown result, |u(q)|2 =
kT 1 1 , Kc (N a)2 q 4
(A3)
where we have also used Na = L, where a is the linear spacing between points and N is the number of points in a linear direction with N2 being the total number of points. The normalized spectrum is then given by Su (qm ) in Eq. (9) with b2 given by Eq. (10). Now let us compare the one-dimensional analogue 2 L 2 d u(x) K E = c1 dx, (A4) 2 0 d 2x which immediately shows that Kc1 has the units of energy times distance, unlike Kc in Eq. (A1) which has units of energy. Again Fourier transform of u(x)
Kc1 L 2 iqm x 2 −iqn x E= −u(qm )qm e −u(qn )qn e dx 2 0 m n (A5) =
L Kc1 u(qm )u(qn )qm2 qn2 ei(qm −qn )x dx 2 m n 0
=
Kc1 L |u(qm )|2 qm4 . 2 m
(A6)
To conform to the two-dimensional case, we now set Kc1 = Kc L. Then, E=
Kc L 2 |u(qm )|2 qm4 . 2 m
(A7)
Applying the equipartition theorem to each mode gives the same result as Eq. (A3) and the subsequent renormalization gives the same Su (qm ) in Eq. (9) with b2 given by Eq. (10). An alternative way to derive the one-dimensional analogue is to assign an infinitely stiff bending modulus Kcy to curvature in the y direction of a two-dimensional membrane while retaining a finite Kcx in the x direction. Then, 2 L
2 2 1 1/2 d u(x, y) 1/2 d u(x, y) E= Kcx + Kcy dxdy, 2 0 d 2x d 2y (A8)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
Any mode with non-zero qn has infinite energy and therefore zero amplitude, so the equipartition theorem directly gives |u(qm )|2 =
1 1 kT , 2 Kcx (N a) qm4
(A10)
equivalent to Eq. (A3). These similarities support using the d = 1 model for our study.
APPENDIX B: COMPARISON OF DISCRETE AND CONTINUUM MODELS FOR FLUCTUATIONS
3
10
1
10
-1
10
0.1
-1
(B1)
−2
As Figure 9 shows, this leads to a q contribution when analyzed for power law behavior. However, the coefficient of the q−2 contribution is only 1/6 as large as the coefficient of the q−4 term, so it only exhibits noticeable devia-
Sρ Continuum Phonons DF Su
4
10
a ~ 0.8 nm
tions from q−4 behavior when q exceeds 1 nm−1 as shown in Figure 9. APPENDIX C: MEAN SQUARE DEVIATIONS IN d = 1 and d = 2
The mean square deviation u2 plays a role in comparing Sρ (q) with Su (q). It is useful to know how this quantity depends upon the size of the system, in both d = 1 and d = 2 dimensions. In one dimension, using Eq. (5) and b = b/N 1 L u2 ≡ |u(x)|2 dx L 0 N/2 b2 1 2b2 N 2 a 4 1 = ∼ N 2. 2 q4 4 4 N (2π ) m m m m=1
(C1)
The final sum depends only weakly upon N, going from 1 at N = 2 to ζ (4) ≈ 1.0823, so u2 scales as N2 . In two dimensions L 1 2 u ≡ |u(x, y)|2 dxdy (L)2 0
Paracrystalline DF Su
10
qB=π/a
FIG. 9. Comparison of undulation spectra when a continuum model is used (solid black line) versus a discrete model (blue squares). The difference (open green circles) is well represented over the range available to simulations by a q−2 term (dashed red line).
=
6
1
q (nm )
Figure 8 addresses the generation of lateral fluctuations with the paracrystalline method versus the phonon method. When Eq. (14) was used for the phonon method, the results were too close to those for the paracrystalline method to be visually distinguishable. Instead, the comparison in Figure 8 is to the continuum phonon case which replaces sin(qa/2) in Eq. (14) with the long wavelength limiting qa/2 behavior. Using the continuum phonon spectrum gives the artifact of a discontinuous derivative at the Brillouin zone boundary compared to the discrete phonon model and the paracrystalline model. The comparison of discrete and continuum phonon generation raises a similar issue regarding transverse fluctuations. If one uses the discrete version for the bending energy, then Suexact (q) = b2 (a/2sin(qa/2))4 .
-4
continuum q discrete discrete - continuum -2 q
5
10 2
and the Fourier transformation gives 1 1/2 1/2 |u(qm qn )|2 (Kcx qm2 + Kcy qn2 )2 (L)2 . (A9) E= 2 m n
undulation spectrum |u(q)|
064114-7
Sρ
=
2
10
N/2 1 4b2 N 2 a 4 ∼ N 2, (2π )4 n,m=1 (n2 + m2 )2
(C2)
so the mean square amplitude diverges as the square of the linear size N and L for both the d = 1 and d = 2 systems.
0
10
-2
10
0.01
0.1
-1
q (nm )
1
10
FIG. 8. Comparison of results obtained with paracrystalline generation and discrete phonon generation of lateral fluctuations are indistinguishable in this figure. The phonon result shown is for phonon generation using a continuum model. Input parameters are the same as in Figure 4 except that σ lat = 0.0018 nm for the phonon case.
APPENDIX D: DERIVATION OF CONDITIONS FOR VALIDITY OF THE DIRECT FOURIER METHOD
In this appendix, we first show that the direct Fourier spectrum SuDF (q) is equal to the true fluctuation spectrum uexact (q) convolved with the lateral density structure factor ρ(q) and multiplied by a constant. We further show that in relevant limits this expression reduces to the additive relation
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
064114-8
Albert, Ray, and Nagle
J. Chem. Phys. 141, 064114 (2014)
(2) posited by Brandt et al. when ρ(q) is uncorrelated with uexact (q). SuDF is defined by 2 N DF 2 1 −iqxn u(xn )e Su (q) ≡ N (D1) . N n=1 This can be written as an integral of the summand times N delta functions
2 N ∞ u(x)e−iqx δ(x − xn ) dx . SuDF (q) = N 2 x=0 n=1 (D2) By the convolution theorem, the Fourier transform of this product is the convolution of the individual transforms of u(x) and the sum of delta functions, denoted w(x). Consequently, 2 ∞ DF 2 u(q )w(q − q )dq , (D3) Su (q) = N −∞ where w(q) is w(q) =
L
x=0
1 N
N
Separating the terms with ρ(0) factors, noting that ρ(0) = 1, and recalling that |u(q)| = |u( − q)|, ⎛ 2|u(qp )|2 SuDF (qm ) = N 2 ⎝|u(qm )|2 + qp >0,qp =qm
×
(|ρ(qm + qp )|2 + |ρ(qm − qp )|2 ) 2
. (D7)
In q regions where ρ does not have significant curvature and can be approximated as the mean of ρ(q ± ) for small , one obtains ⎞ ⎛ 2|u(qp )|2 ⎠ SuDF (qm ) ≈ N 2 ⎝|u(qm )|2 + |ρ(qm )|2 qp >0 =q
≈ Suexact (qm )(1 − |ρ(qm )|2 ) + Sρ (qm ).
(D8)
In our emulations and likely for actual simulations, |ρ(q)|2 is small compared to unity, so this reduces to Eq. (2).
δ(x − xn ) e−iqx dx
1 Y.
n=1
N 1 −iqx n = ρ(q). = e N n=1
(D4)
The ρ(q) in Eq. (D4) is exactly the lateral density fluctuation components in Eq. (11). Thus, the direct Fourier spectrum SuDF in Eq. (D1) is equal to the square of the convolution of the true undulation spectrum u(q) with the lateral density spectrum ρ(q), multiplied by a factor of N2 . The analysis so far has been identical to that done for non-equispaced time series by Ref. 23. Unfortunately, as noted in Ref. 23, deconvolving Eq. (D3) to get the u is not always feasible when working with actual, noisy data. We next demonstrate that Eq. (D3) reduces to Eq. (2) proposed in Ref. 7, if lateral fluctuations are sufficiently small and statistically uncorrelated with undulations. Crucially, if either of these conditions is violated, the additive relation Eq. (2) will not necessarily hold. Substituting Eq. (D4) into (D3) and using the assumption of periodic boundary conditions on u(q) gives 2 DF 2 u(qp )ρ(qm − qp ) . (D5) Su (qm ) = N qp Assuming the u and ρ modes are uncorrelated, then SuDF (qm ) = N 2 |u(qp )|2 |ρ(qm − qp )|2 .
(D6)
Shibata, J. J. Hu, M. M. Kozlov, and T. A. Rapoport, Annu. Rev. Cell Dev. Biol. 25, 329 (2009). 2 L. V. Chernomordik and M. M. Kozlov, Annu. Rev. Biochem. 72, 175 (2003). 3 Y. Liu and J. F. Nagle, Phys. Rev. E 69, 040901 (2004). 4 E. Lindahl and O. Edholm, Biophys. J. 79, 426 (2000). 5 R. Goetz, G. Gompper, and R. Lipowsky, Phys. Rev. Lett. 82, 221 (1999). 6 J. F. Nagle, M. S. Jablin, S. Tristram-Nagle, and K. Akabori, “What are the true values of the bending modulus of simple lipid bilayers?,” Chem. Phys. Lipids (in press). 7 E. G. Brandt, A. R. Braun, J. N. Sachs, J. F. Nagle, and O. Edholm, Biophys. J. 100, 2104 (2011). 8 W. K. den Otter and W. J. Briels, J. Chem. Phys. 118, 4712 (2003). 9 J. Stecki, J. Chem. Phys. 120, 3508 (2004). 10 J. Stecki, Adv. Chem. Phys. 144, 157 (2010). 11 J. Stecki, J. Chem. Phys. 137, 116102 (2012). 12 E. R. May, A. Narang, and D. I. Kopelevich, Phys. Rev. E 76, 021913 (2007). 13 M. C. Watson, E. G. Brandt, A. J. Welch, and F. L. H. Brown, Phys. Rev. Lett. 109, 028102 (2012). 14 W. Helfrich, Z. Naturforsch., C 28, 693 (1973). 15 W. Helfrich and R. M. Servuss, Il Nuovo Cimento D 3, 137–151 (1984). 16 G. Brannigan and F. Brown, Biophys. J. 90, 1501 (2006). 17 M. C. Watson, E. S. Penev, P. M. Welch, and F. L. H. Brown, J. Chem. Phys. 135, 244701 (2011). 18 F. N. Fritsch and R. E. Carlson, SIAM J. Numer. Anal. 17, 238 (1980). 19 R. Hosemann and S. N. Bagchi, Direct Analysis of Diffraction by Matter (North-Holland, Amsterdam, 1962). 20 A. Guinier, X-Ray Diffraction (W. H. Freeman, San Francisco, 1963). 21 J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. Acta, Rev. Biomembr. 1469, 159 (2000). 22 A. R. Braun, E. G. Brandt, O. Edholm, J. F. Nagle, and J. N. Sachs, Biophys. J. 100, 2112 (2011). 23 T. J. Deeming, Astrophys. Space Sci. 36, 137 (1975).
qp
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.24.51.181 On: Sat, 22 Nov 2014 19:19:10
