Determining the Lipid Tilt Modulus by Simulating Membrane Buckles

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Determining the Lipid Tilt Modulus by Simulating Membrane Buckles Xin Wang, and Markus Deserno J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b02016 • Publication Date (Web): 12 Apr 2016 Downloaded from http://pubs.acs.org on April 17, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Determining the Lipid Tilt Modulus by Simulating Membrane Buckles Xin Wang and Markus Deserno∗ Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA E-mail: [email protected] Phone: +1 412 268 4401

Abstract

of the pivotal plane for Berger DMPC lies just below the glycerol backbone—unlike for MARTINI DMPC, where this plane lies closer to the middle of the lipid.

Lipid tilt affects the energy of membrane deformations on scales comparable to the membrane’s thickness. The surface divergence of the tilt field creates a local spontaneous curvature, while tilt itself is quadratically penalized with a strength given by the tilt modulus. Traditionally, this modulus is determined by measuring the power spectrum of lipid orientation fluctuations. Here we present a novel approach which does not rely on fluctuations but instead exploits the fact that curvature gradients induce a tilt field. Its implementation extends a technique previously developed by us for localizing the position of the pivotal plane in buckling simulations, which quantifies the lipid imbalance across segments cut out from a complete buckle. Lipid tilt affects this count in a predictable way, and the signal can be quantified well enough to back out the tilt modulus—at no additional cost and with about 5% precision for not too coarse models. We apply our technique to three lipid models of very different resolution: the highly coarse grained Cooke model, and two versions of DMPC, using both the (less highly coarse grained) MARTINI and the (united atom) Berger force field. For Cooke we find an effective bilayer tilt modulus of (29 ± 9) pN/nm, and for the less generic DMPC lipid we find (115±6) pN/nm for MARTINI and (39 ± 2) pN/nm for Berger, both in reasonable agreement with existing studies for these models. We also show that the position

Introduction On length scales not much larger than molecular, lipid membranes can be described remarkably well by simple continuum theories that contain only a small number of elastic parameters. At sufficiently large dimensions a membrane’s curvature deformation and de facto inextensibility are all that remains relevant; this is the regime of classical Helfrich theory. 1–3 But on scales not much greater than membrane thickness, say O(10 nm), other collective fields matter as well. One of them is lipid tilt, the deviation of a lipid’s orientation away from the local membrane normal (or a suitably defined local spatial average thereof). For any set of order parameters (geometry, tilt, or others), elastic Hamiltonians can be deduced based on symmetry principles (“topdown”). For instance, even though Helfrich motivated the discussion in his seminal paper by drawing analogies to liquid crystal theory, he argued for the final curvature elastic energy functional now associated with his name largely on the basis of symmetry. 2 Gompper and Klein have constructed a Ginzburg-Landau theory for lipid membranes, 4 combining a scalar (density) and a vectorial (tilt) order parameter. And Nelson and coworkers have discussed tilt (and more

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 22

plementary method is to impose an active deformation and measure the force required to do so. For example, the bending rigidity of membranes can be deduced from the force needed to pull or hold thin membrane tethers, a method that was first proposed experimentally 51–53 and later recreated in a simulation framework. 50,54–56 Unfortunately, simulating equilibrated tethers is challenging, because both the chemical potential of lipids between the two leaflets, as well as the chemical potential of water between the inside and the outside of the tether is hard to balance; in the first case due to the low lipid flipflop rate, in the second case due to the low permeability of the bilayer to water. However, an alternative active deformation approach, avoiding both problems, has recently been proposed: membrane buckling. 57–59 The investigations in the present paper also rely on buckling, but our starting point is a technique that aims to measure a parameter that is not energetic but geometric in nature. We have recently shown that the position of the pivotal plane can be measured by quantifying a suitably defined local lipid imbalance across the leaflets of a buckled membrane. 60 Here we first of all extend this analysis to an essentially all-atom model of the lipid dimyristoylsn-glycerophosphocholine (DMPC). More importantly, though, we show that a simple nocost extension of our previous analysis also gives access to the lipid tilt modulus. This connection arises because tilt is created by gradients in membrane curvature, which in turn affect the precise lipid counting across leaflets. Since the latter can be done with very high precision, the extent of lipid tilt can be quantified. But the latter can also be predicted from the average membrane shape and the value of the tilt modulus. Comparing measurements and prediction hence gives access to that modulus.

complicated order parameter fields) based on symmetry, adding additional field-theoretical insights about the relevance of such degrees of freedom. 5–8 Alternatively, one can obtain larger-scale effective theories from underlying smaller-scale models (“bottom-up”). For instance, meanfield techniques adapted from polymer theory have been used to describe the lipid tails, and from there derive a larger-scale curvatureelastic Hamiltonian for which the values of the moduli are predicted based on more microscopic observables. 9–14 Hamm and Kozlov have derived an energy functional for lipid monolayers, which accounts for both curvature elasticity and lipid tilt, by a suitable dimensional reduction, on the basis of a quadratic elastic theory that describes thin in-plane fluid sheets subject to internal lateral pre-stresses. 15,16 Their Hamiltonian has been developed further by May et al. 17,18 and more recently by Brown and coworkers. 19–23 In such approaches, the elastic moduli (say, for bending or tilting) arise either as free parameters (top-down), or are linked to lower-level elastic properties, such as the elastic tensors of the fluid-elastic thin sheet or its underlying pre-stresses (bottom-up). Either way, their numerical values are not directly predicted and hence must be measured—both in experiment and in computer simulations. In the top-down case this is clear; in the bottom-up case this happens because the lower-level observables to which the theory relates (say, the trans-bilayer pre-stresses) are often not themselves predicted by it or otherwise known. The aim of the present paper is to propose a novel method for measuring the tilt modulus in simulations by means of active deformations, and apply it to three different lipid models of rather different resolution. Thermal fluctuations of both bending and tilt degrees of freedom are usually large enough so that measuring the resulting fluctuation spectra can be used as a way to deduce the moduli. For instance, monitoring the spectrum of membrane undulations is a standard procedure for obtaining the bending modulus in experiment 24–34 and in simulation. 17–20,35–50 A com-

Theoretical Framework Helfrich Hamiltonian The classical model of curvature elasticity expresses the energy of a membrane as an integral


2

Page 3 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hamiltonian becomes 16 Z 2 1 H = dA κm K + ∇ · T − K0,m 2 1 2 ˜ . (4) +¯ κm KG + κt,m T 2

over its area, with an energy density quadratic in the curvatures. For each leaflet (i.e., a single monolayer) we can write down the energy Z 1 2 H = dA κm (K − K0,m ) + κ ¯ m KG . 2 (1) Here, K = Tr(K) = c1 + c2 and KG = det(K) = c1 c2 are the total and Gaussian curvature— trace and determinant of the curvature tensor K, whose eigenvalues ci are the local principal curvatures; 61 K0,m is the monolayer’s spontaneous curvature, and the two monolayer moduli κm and κ ¯ m multiply the two independent quadratic curvature invariants. Essentially the same functional form also holds for bilayers. The reason we start out with a single leaflet is that the generalization to tilt is more natural on the monolayer level. We will later discuss a simplified effective bilayer version.

˜ G is Since the difference between KG and K of higher order (it contains terms quadratic in both tilt and curvature), we can ignore the tilde to lowest order. In fact, since we will subsequently be interested in planar buckles, for which KG = 0, the Gaussian term will vanish identically. A functional variation of Eqn. (4) with respect to T gives the Euler-Lagrange equation satisfied by the tilt: q ∇ ∇·T −ℓ−2 T = −∇K with ℓ = κm /κt,m , (5) where the characteristic length ℓ describes the range over which a tilt excitation decays. Since ℓ is on the order of a nanometer, tilt-physics is indeed very local. Eqn. (5) shows that if we impose curvature gradients on a membrane, we will induce a tilt field.

Lipid tilt within the model of Hamm and Kozlov The Helfrich Hamiltonian (1) does not contain lipid tilt, i.e. the fact that a lipid’s orientation n might deviate from the local membrane normal N . This degree of freedom can be described by an additional tangential vector field T , defined as n −N , (2) T = n·N such that the tilt angle θ is given by tan θ = |T | .

Lipid tilt in the geometry of a buckled membrane An undulatory membrane buckle varies only along one spatial direction (say, the xdirection), while it is flat along the other one (say, y). Hence, the problem becomes effectively one-dimensional. At a given shape, the Euler-Lagrange equation for the tilt is then an ordinary differential equation for the tangential tilt field T (s), measured along and as a function of arc length s:

(3)

Hamm and Kozlov have argued that including T in the elastic description will introduce two major changes. 16 First, tilt itself will be quadratically penalized by a tilt modulus κt.m . And second, the curvature tensor K gets re˜ = placed by a “dressed” curvature tensor K K + ∇ ⊗ T (where ∇ is the vector differential operator on the surface). This changes the total ˜ = Tr(K + ∇ ⊗ T ) = ˜ = Tr(K) curvature into K K + ∇ · T (meaning, the tilt divergence adds to the curvature), and the Gaussian curvature into ˜ Hence, their joint curvature-tilt ˜ G = det(K). K

T ′′ − ℓ−2 T = −K ′ ,

(6)

showing that tilt is sourced by curvature gradients. The Green function of the left hand side is GT (s) = − 12 ℓ e−|s|/ℓ , showing that, indeed, tilt decays on the length-scale ℓ, and that the exci−1/2 tation magnitude is proportional to ℓ ∼ κt,m . For a buckle in angle-arclength parametrization K(s) = −ψ ′ (s), where ψ(s) is the angle


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

which the buckle locally makes with the horizontal. If we ignore tilt, the Euler Lagrange equation for the geometry of the buckle is given by 57,58 p ψ ′′ + λ−2 sin ψ = 0 with λ = κm /fx , (7)

where λ is the stress-dependent characteristic length from Eqn. (7) and m is the elliptic parameter; both are fixed by the boundary conditions. The function sn[x, m] is a Jacobi elliptic function. 63 Using some elementary identities for such functions, we can rewrite the sine of the angle as

where fx is the buckling stress (per monolayer). Combining this with the Euler-Lagrange equation for the tilt, we get T ′′ − ℓ−2 T = −λ−2 sin ψ ,

Page 4 of 22

√ d sin ψ(s) = −2 m λ cn[s/λ, m] . ds

(10)

Our aim is to expand the right hand side in a Fourier series in s. To do so, let us recall that elliptic functions are doubly periodic—in both the real and the imaginary direction. The quarter period in the real direction is given by K = K[m], and the one in the imaginary direction is K′ = K[1 − m], where K[m] is the complete elliptic integral of the first kind. 63 From these periods we can define the so-called “nome” q = exp{−πK′ /K}. It now turns out that all elliptic functions can be written as Fourier expansions involving powers of the nome. Specifically, we are interested in cn[x, m], for which the expansion is 63

(8)

an equation which can be solved for the tilt along the buckle. Had we not ignored tilt in the Euler-Lagrange equation of the shape, we would not have ended up with a stand-alone equation for ψ(s), and we would instead have to solve for ψ(s) and T (s) simultaneously. Removing the influence of tilt on shape clearly makes for easier equations— but why would this be allowed? The answer is that tilt is such a stiff degree of freedom (as evidenced by the smallness of ℓ compared to the geometric length L of the whole buckle) that it only modifies the shape equation very minimally. We will provide a more careful argument for this below. Since we know the solution ψ(s) for the buckle, 57,58,62 we could simply convolve it with the Green function of the tilt equation and arrive at T (s). Yet, proceeding this way is cumbersome, as the right hand side of Eqn. (8) is a fairly complicated expression involving Jacobi elliptic functions. But we know that it has periodicity L, where L is the contour length of the buckle, and so we can expand it as a Fourier series, for which there turns out to exist an analytical expression. This is not only easier, it will also suggest some convenient approximations.

∞ 2π X q n/2 πnx cos . (11) cn[x, m] = √ n 2K K m n=1 1 + q n odd

Using this in Eqn. (10), we get the series expansion ∞ 2π 2 X n q n/2 2πns sin ψ(s) = 2 sin . n K n=1 1 + q 4Kλ

(12)

n odd

Since we furthermore know that this function has period L, it must be true that L = 4Kλ ,

(13)

and with this we can eliminate the occurrence of K in the argument of the sine, rewriting the Fourier series expansion in such a way that the period L becomes manifest:

Fourier expansion of a buckle The exact analytical solution for the bilayer buckling problem is well known, 57,58,62 and the angle as a function of arc length can be written as n√ o ψ(s) = 2 arcsin m sn[s/λ, m] , (9)

∞ 2πns 2π 2 X n q n/2 sin . sin ψ(s) = 2 n K n=1 1 + q L n odd


4

(14)

Page 5 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Coefficients An (γ) for the series expansions introduced in Eqns. (16) and (17). n

will search for the solution in the form T (s) =

n=1 n odd

An (γ)

1

1+

3

1+

19 2 7 2867 1 γ + 512 γ + 512 γ 3 + 524 γ4 8 288 9 159 2 1411 3 γ + 512 γ + 8182 γ + ··· 16

5

1+

15 γ 16

7

1+

21 γ 16

9

1 + ···

355 2 γ 512

+

+ ···

+ ···

Cn (γ) =

+ ···

2π L

2 X ∞

n=1 n odd

Bn (γ) sin

2πns , (16) L

γ n−1 2 , Bn (γ) = 2 γ An (γ) n 16

n2

Bn (γ) . + (L/2πℓ)2

B1 ( 12 ) ≈ 1.52 B3 ( 12 ) ≈ 0.184 B5 ( 12 ) ≈ 0.0119 B7 ( 12 ) ≈ 0.000629 B9 ( 12 ) ≈ 0.0000281 .

where the coefficients Bn are defined by √

2πns . L

(18)

(19)

Besides the (good) approximation that we did not let the tilt field affect the shape in the first place, this solution is exact. However, the desire to reap a very substantial technical convenience suggests one further approximation: notice that in the denominator of Eqn. (19) the term (L/2πℓ)2 is large if the tilt modulus is large—recall that ℓ is microscopic but L is macroscopic. Hence, for the first two modes n = 1 and n = 3, and maybe even for the third mode n = 5, the competing term n2 can be neglected. And we may hope that for even higher modes the coefficients Bn (γ) become so small, mostly because of the exponentially decreasing term (γ/16)(n−1)/2 , that we can ignore the mode (and all subsequent ones) altogether. For instance, at γ = 12 , the first few values of the Bn are

1 1 11 4 17 5 m(γ) = γ− γ 2 − γ 3 − γ − γ −· · · . 8 32 1024 4096 (15) From this we obtain the coefficients of the Fourier expansion as a Taylor series in γ:

Cn sin

Inserting this into Eqn. (8) and comparing terms with Eqn. (16), we find the coefficients

This expansion still contains the elliptic quarter periods K and K′ , which in turn depend on the elliptic parameter m. We have previously shown that the latter can be written as a rapidly converging series expansion in the buckling strain γ = (L−Lx )/L. 58 The first few terms are

sin ψ(s) = λ2

∞ X

(17)

(20a) (20b) (20c) (20d) (20e)

This is interesting because if we scratch the term n2 in the denominator, what remains are—up to some overall factor—the Fourier coefficients of sin ψ. This means that to some decent (but admittedly hard to control) approximation, we have the following result: 2 2 2K 2πℓ T (s) ≈ sin ψ(s) = ε(γ) sin ψ , π L (21)

and the Taylor series An (γ) for n ≤ 9 (i.e., the first 5 modes) are listed in Table 1 (the notation is chosen such that An (0) = 1 for all n).

Determining the tilt field We can now come back to the differential equation (8). Since we have the right hand side as a Fourier series, and the equation is linear, we


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

were we defined 2 2πℓ 1 9 2 ε(γ) = 1 + γ + γ + ··· . L 2 32 (22) This expression for ε(γ) results from inserting the series expansion of m(γ) from Eqn. (15) into K = K[m(γ)] and expanding again. Being the amplitude in front of sin ψ(s), ε(γ) is the largest possible value of the tilt; the largest one that is realized occurs at the inflection points of the buckle, where ψ(s) = ψi takes its maximum value. From Eqn. (9) or Eqn. (10), combined with Eqn. (13), we can readily deduce that p sin ψi = sin ψ(L/4) = 2 m(1 − m) , (23)

Page 6 of 22

so we return to the condition ℓ ≪ L—which is precisely what motivated us to drop the term n2 in the denominator of Eqn. (19).

Tilt affects the energy only on small scales Let us finally return to the question why we can in good approximation neglect the back effect of tilt on membrane shape. The answer hinges on the fact that for curvatures small compared to ℓ−1 the contribution of tilt to the membrane energy is strongly suppressed, and consequently tilt affects the equilibrium large scale shape only very little. To see this, consider a membrane buckle with a single mode of wave vector k = 2π/L and small amplitude h0 ,

and inserting the series expansion for m(γ) from Eqn. (15), we arrive at 9 25 2 √ sin ψi = 2 γ 1 − γ − γ − ··· . 16 512 (24) Together with the series expansion for ε(γ) from Eqn. (22), this gives a maximum tilt amplitude of 1 Tmax (γ) 25 2 √ =2 γ 1− γ− γ − ··· . (2πℓ/L)2 16 512 (25) Since the right hand side becomes 1 for γ ≈ 0.26, this result shows that (2πℓ/L)2 can be viewed as the maximum tilt which lipids assume in a buckle strained by about a quarter of its length. From Eqn. (3), the maximum tilt angle is hence θmax = arctan(Tmax ). There is another useful way to look at the approximation (21): combining it with the EulerLagrange equation for the geometry (7) and using the identity (13), we see that it is equivalent to T (s) = ℓ2 K ′ (s) . (26)

h(x) = h0 sin(kx) .

(27)

Within linearized Monge gauge, its curvature is given by K(x) = −h′′ (x) = h0 k 2 sin(kx) .

(28)

Since the tilt field T (x) satisfies the equation T ′′ − ℓ−2 T = −K ′ = −h0 k 3 cos(kx) ,

(29)

we immediately see that it is given by h0 k 3 . k 2 + ℓ−2 (30) A simple calculation now leads to the energy density of this system averaged over one period: Z o 1 L nκ κt 2 dx e= [K(x) + T ′ (x)] + T 2 (x) L 0 2 2 1 (31) = κeff (k) hK 2 i , 2 T (x) = T0 cos(kx) with T0 =

with the effective scale-dependent rigidity

This is just the Euler-Lagrange equation for the tilt, (6), but with the second derivative term removed. Now, T ′′ ≪ T /ℓ2 holds if the tilt varies slowly compared to the decay length ℓ, and this in turn happens when the function sourcing the tilt, K ′ , varies sufficiently slowly. But the latter rate is set by the period of the buckle, and

κeff (k) =

κ κ = . 2 2 1+k ℓ 1 + (2πℓ/L)2

(32)

Since ℓ is a microscopic length, and L is macroscopic, this softens the bare curvature modulus only very little. Notice that (up to minor cor-


6

Page 7 of 22

20 !

rections) (2πℓ/L)2 = ε(0) is the slope of z0 (ξ), which we use to extract the tilt modulus, and which we found to be a few percent. Curvature softening due to tilt is hence an effect on the order of a few percent. Notice, though, that it will be larger if the tilt correlation length ℓ is larger. Could this happen, for instance in buckles made from more densely packed gel phases? 59 This is difficult to say, because ℓ2 = κm /κt,m depends on both the bending and tilt modulus. While we know that the bending modulus increases by about an order of magnitude in gel phases, 64–66 it is unclear what the tilt modulus does. One could either imagine it growing as well, given that the system becomes denser and the area per lipid shrinks; but one could also imagine it staying the same or even decreasing, since tilt measures transverse shear (i.e., the vorticity direction is tangential to the membrane), and this shear modulus need not increase if the Young modulus does—in fact, it could decrease due to the higher tail order of lipids that are sheared against one another. The fate of the length ℓ is hence not obvious and would be worth studying further. Should it increase significantly, the approximation that tilt does not affect the shape would have to be revisited. As an interesting side-note: observe that 1 1 1 = + , 4 4 κeff (k) k κk κt k 2

x(s), z(s)

15

" ψ(s)

z [nm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10

s

s

N+

5

ψ(s)

N−

0 Lx -5

-5

0

5

10

15 20 x [nm]

25

30

35

40

Figure 1: Snapshot of a MARTINI buckle simulation viewed in the xz–plane. The solid curve follows the bilayer’s midplane. Along its arc length s we pick two points—one at s and one at s—that cut out a segment for which we count the number of lipids N± on the upper and lower leaflet. The boundary angles ψ(s) and ψ(s) are also indicated. Reprinted from Wang and Deserno 60 with the permission of AIP Publishing.

first, we ignore tilt (since it is small) and calculate the shape. Next, we calculate the tilt that is induced by this shape, which is of order ε ∼ (2πℓ/L)2 . So far this is exactly what we do. Next, correct the shape by solving the shape equation while using the small tilt field just calculated as a fixed input. This will give rise to a small first order correction. This correction in turn induces a correction to the previously calculated tilt field, but clearly this is now higher order, ε2 , and so we can ignore it at the level of accuracy that we need.

(33)

showing that the term on the left hand side, which typically occurs as the amplitude of a bilayer’s undulatory power spectrum, gets an additional tilt contribution at large wave-vectors if the effect of tilt on bending is considered. Our ground-state analysis hence provides an independent illustration of an effect well-discussed in the literature of fluctuating bilayers. 22,23,34 One might worry that the whole point of our analysis was to look at very small deviations of the area difference across the two leaflets, in order to extract z0 and the tilt modulus. Would this not require us to solve for the shape more carefully, for self-consistency reasons? The answer is no, because the small effect we are after is the slight tilt due to an imposed curvature. Imagine solving the problem iteratively:

From the pivotal plane to the tilt modulus In a previous publication we have shown how the simulation of a buckle can be used to extract the pivotal plane of the monolayer. 60 Here we show how an extension of our reasoning can not only explain a curious finding from our previous analysis, but in fact give access to the monolayer’s tilt modulus κt,m .


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 22

Pivotal plane—without tilt The basic idea of finding the position of the pivotal plane via a buckled membrane configuration is illustrated in Fig. 1: If we cut a section out of a buckle, then the number of lipids in the upper and lower leaflet differ, because the total area of these two leaflets is different by an amount that is proportional to both the integrated curvature and the distance of the pivotal plane to the bilayer’s midplane. If a± and N± are the actual area per lipid and the number of lipids in the upper (+) and lower (−) leaflet, respectively, then we have previously shown that∗ a± N± = ∆s ∓ z0 ∆ψ , Ly

T

n e N

(34) Figure 2: Illustration of lipid tilt fields in a buckled lipid bilayer. Solid black curves are the outer surfaces of the bilayer, the dashed curve is the bilayer’s midplane, and the two dotted curves are the pivotal planes of each leaflet. Local normals are indicated by solid lines perpendicular to the midplane. Tilted lipids are illustrated as red and blue lines in the upper and lower leaflet, respectively; the extent of tilt is strongly exaggerated. The green outline is an example of a cut-out segment. Notice that for this particular cut the head groups of the lipids in the upper leaflet tilt into the cut-out segment, while the head groups of the lipids in the lower leaflet tilt out of the cut. The small inset illustrates the tilt definition from Eqn. (2).

where ∆s is the arc length between the two cuts, ∆ψ is the angle by which the buckle changes along the cut, z0 is the distance of the pivotal plane from the bilayer’s midplane, and Ly is the width of the buckle. We previously used this equation to extract the pivotal plane position z0 . 60

Pivotal plane—with tilt The key insight—as far as tilt is concerned—is that tilt further increases or decreases the available area. As Fig. 2 illustrates, different regions along the height of a lipid may tilt into or out of any chosen cut, thus changing the count of lipids in that region. Moreover, the effect does not cancel but amplify across the two leaflets, thus affecting the overall lipid imbalance in a tilt-dependent way. Considering the geometry outlined in Fig. 2, we arrive at the following tilt-correction:

where ξ is the distance away from the bilayer’s midplane which we use as a reference for lipid counting. For instance, if we use a particular bead as the reference bead, then ξ is the (average) distance of that bead from the midplane. As it stands, this result is not yet very useful. We do not measure the tilt in our simulation, and in fact it is very hard to do so; we come back to this in the Discussion section. Hence, we will use Eqn. (26) to link the tilt to the local curvature gradient. However, there is yet one more complication we need to consider: the curvature we speak of is the one of the bilayer midplane, but the monolayer tilt fields couple to the curvatures Kp,± at the two pivotal planes,

a± N± = ∆s ∓ z0 ∆ψ + (ξ − z0 )[T± (s) − T± (s)] , Ly (35) ∗

Actually, the equation in Ref. 60 contains a “±” instead of the “∓”. This derives from a different sign convention for the relation between curvature and an˙ “∓” is the right gle. Using the convention K = −ψ, sign. It should be noted, though, that as long as one is consistent, this will at the end only result in one overall minus sign, which causes no trouble.


8

Page 9 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


not the curvature K at the midplane. However, these curvatures are connected by the parallel surface relation 60,67 K ± 2KG z0 1 ± Kz0 + KG z02 K KG =0 = 1 ± Kz0

Kp,± =

Kz0 ≪1

≈

K ∓ K 2 z0 ,

form ∆(Kψ) yp = z0 1 + ε 1 − 2(ξ − z0 ) − εξ . xp ∆ψ (41) This leads to the following analysis prescription: create a large number of cut-out segments for one individual buckling snapshot and calculate the variables xp and yp for each segment. This yields a collection of points {xp , yp } that define a line through the origin, whose slope will lead to z0 . Averaging over many snapshots will increase precision and provide error bars. The only remaining nuisance is the correction term ∆(Kψ)/∆ψ. Fortunately, it turns out that we can ignore it. Notice that this term emerges from the last contribution to Eqn. (39), which (with our choice of signs in this equation) (p) gets multiplied by M− , the difference in lipid number across the whole buckle. Since we al(p) ways set up initial buckles such that M− = 0, it can only ever become nonzero due to lipid flip-flop, but for many lipid models—especially the more highly resolved ones (such as MARTINI and Berger)—the flip-flop rate is so small that the initial lipid symmetry is preserved throughout the entire simulation, and thus the nuisance term vanishes exactly. If the flip flop rate does not vanish (such as in the Cooke (p) model), then M− will not remain zero during the simulation—but its average will. Hence, we can remove the nuisance term from Eqn. (41), and arrive at the very simple relation

(36a) (36b) (36c)

showing that the tilt fields should in fact be written as T± = ±ℓ2 K ′ (1 ∓ 2Kz0 ) .

(37)

Inserting this into Eqn. (35) and using Eqns. (21,26), we arrive at a± N± = ∆s ∓ z0 ∆ψ ± ε(γ)(ξ − z0 )∆ψ Ly − 2ε(γ)(ξ − z0 )z0 ∆(Kψ) . (38) Notice that we used the further approximation sin ψ ≈ ψ, in order to subsequently be able to combine the tilt correction with the one that stems from pure geometry. To extract the two lengths z0 and ℓ (the latter enters via ε) from Eqn. (38), it helps to define two further quantities: first, M± = N+ ± N− , the sum/difference of the lipid numbers across (p) the cut-out segment, and M± , this quantity evaluated not for the cut-out segment but for the whole buckle. With these definitions, we can rewrite Eqn. (38) as (p) M± ∆s

− M± L

z0 =

(p) = M∓ ∆ψ[z0 (1 + ε) − εξ] (p) − 2M± ǫz0 (ξ − z0 )∆(Kψ)

. (39)

(p)

(p)

(42)

This tells us that the pivotal plane position z0 is almost equal to the slope yp /xp , but that it will very slightly depend on the choice of reference bead, since that bead’s average position ξ away from the bilayer midplane enters the value for z0 : the farther away the reference bead is, the larger the value of z0 we will deduce. The self-consistent choice ξ = z0 leads to z0 = yp /xp , and this is indeed the prescription we followed in our previous work. 60 But there is more information that can be extracted from the function z0 (ξ). In our earlier work we had not yet understood the origin of its non-

(p)

Since Ntotal = M+ ≫ M− ≈ 0, we will now choose the lower tier signs in order to preserve the information leading to z0 and ε. Introducing the two variables (p)

and yp = M− ∆s − M− L , (40) we can turn Eqn. (39) into the (almost) final xp = M+ ∆ψ

ε yp /xp + ξ. 1+ε 1+ε


9

The Journal of Physical Chemistry choline

model, a 3 beads per lipid implicit solvent representation, which chiefly aims to capture the phenomenology of amphipathic surfactant assembly into fluid bilayer membranes. 44,45,68 Our medium-resolution model is a 10 bead representation of DMPC using the MARTINI model, which aims to capture the rough geometry of a lipid molecule as well as the solvation properties of its components. 69,70 And at the finest level, we use a 46 united atom representation of DMPC using the Berger force field. 36,71 All three models are extensively described in the literature, and we have discussed their use in buckling simulations previously. 58 In fact, the simulation trajectories we analyze are the same as the ones used in Ref. 58 , and so we refer the reader to that reference for all technical details. Here we merely point out a few key aspects of these models, relevant for understanding their behavior under the present conditions. Fig. 3 provides a quick visual overview of the three models. The Cooke model is an implicit solvent model, so assembly is driven by a suitably chosen attraction between the tails. The depth of that interaction potential, ǫ, is the basic unit of energy; we choose it such that kB T = 1.1 ǫ. The range of the potential features a length parameter wc , which we choose to be wc = 1.6 σ, where σ is the excluded volume parameter of the tail beads’ Weeks-Chandler-Andersen potential. Due to its strongly coarse-grained nature, the Cooke model has a lipid flip-flop rate that is large enough to equilibrate the lipid distribution across leaflets during a simulation. We ran the simulations using ESPResSo. 72 The MARTINI model maps on average four heavy atoms to one CG bead and can describe many different lipids by suitably choosing tail lengths and interaction potentials; an explicit CG solvent combines four atomistic water models into one effective solvent bead. We represent DMPC by 10 beads (3 beads per tail) and ran the simulations with GROMACS 4.5, 73 employing a Berendsen thermostat 74 and a time constant τT = 1 ps at a reference temperature T = 300 K. The Berger force field is a united atom atomistic force field for lipids in which nonpolar hy-

head

choline phosphate glycerol backbone

phosphate glycerol backbone

tail

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sn 2 chain sn 2 chain

sn1 chain

sn1 chain

Figure 3: Illustration of the three different lipid models used in this study, from left to right: Cooke, MARTINI DMPC, Berger DMPC. The somewhat non-standard colors are chosen to help identify structural parts. vanishing slope z0′ (ξ)—hence the need for the self-consistency argument. The analysis above now shows that this slope derives from ε, which in turn is directly related to the tilt modulus. Using Eqn. (22), and simplifying slightly (because ǫ ≪ 1), we therefore get κt,m = κm

2π L

2

Page 10 of 22

9 2 γ + ··· 1 + 12 γ + 32 . (43) ′ z0 (ξ)

This equation permits us to determine the tilt modulus from buckling simulations—at no additional cost, since we get z0′ (ξ) for free in the analysis leading to z0 , 60 and even κm = κ/2 can be extracted from the buckle. 58 In the next section we will apply this analysis to three different lipid models, relying entirely on trajectories previously calculated. 58

Results Description of the lipid models Let us first summarize the essential features of the three different lipid models we use in this study—Cooke, MARTINI, and Berger. The results of all simulations conducted with them are summarized in Table. 2. We use three different lipid models in this paper, which differ markedly in their resolution— from highly coarse grained (CG) to essentially atomistic. At the CG level, we have the Cooke


10

Page 11 of 22

The Journal of Physical Chemistry 1.4

tions, from which the value of z0 can be determined as a function of reference bead position ξ (data not shown). Since DMPC in the Berger force field contains 46 heavy atoms per lipid, this gives 46 points in a plot of z0 (ξ). Fig. 4 shows such a plot for the buckling strain γ = 0.24. Unlike for the MARTINI model, these points do not all fall onto a straight line. Instead, the tail region shows a linear relation with fairly good accuracy, while beyond the glycerol backbone the slope noticeably decreases and ultimately even changes sign. A self-consistent determination of z0 (i.e., intersecting the measured data with the line z0 = ξ) finds the pivotal plane at the end of the linear region, close to the glycerol backbone. A more careful analysis, averaging over all buckles and all values of γ, gives z0 = 1.3225(44) nm, placing the pivotal plane roughly at the carbonyl carbon in the sn1 chain. The error is probably underestimated, given that the goodness of fit is fairly small, but even the outlier data points are not farther away than a single chemical bond length. Besides the change in pivotal plane position, the noticeably nonlinear shape of z0 (ξ) has further interesting physical implications—not just for measuring the tilt modulus, but for defining what exactly we mean by lipid tilt. We will come back to this in the Discussion section.

1.35

z0 [nm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.3 1.25 1.2 1.15

0.5

1 !ξi " [nm]

1.5

Figure 4: Position of pivotal plane z0 for the Berger DMPC model, plotted against the average position hξi i of the ith reference atom in the lipid, for the strain γ = 0.24. The solid line is a linear fit to the 28 carbon atoms in the two myristoyl tails, flanked by the 1 σ and 2 σ confidence bands. The dashed line is the identity z0 = hξi i, which intersects the data at the self-consistently determined value of z0 . drogen atoms are grouped together with carbons into CH, CH2 , CH3 beads. 36,71 As water model we used SPC. The isothermal-isobaric ensemble was achieved by combining a NoseHoover thermostat 75,76 (τT = 1.0 ps at T = 300 K) and a Parrinello-Rahman barostat 77,78 (τP = 50.0 ps, κT,x = κT,y = 0, κT,z = 5 × 10−5 Pa, at 1 bar). The initial buckled configurations were backmapped from MARTINIDMPC simulations using a special version of GROMACS 79

Monolayer tilt modulus To calculate the monolayer tilt modulus from Eqn. (43), we need to know (a) the monolayer bending modulus κm and (b) the slope of the line z0 (ξ). For all three models we have previously calculated the bilayer bending modulus κ via the buckling protocol, 58 and since κm = κ/2 we are done in terms of moduli. For the calculation of z0 (ξ), we re-use the same trajectories created during the buckling studies, 58 and subject them to the novel analysis described in this paper.

Pivotal plane We have previously determined the pivotal plane position z0 both for the Cooke model and for MARTINI-DMPC. 60 For Cooke, we find z0 = 1.4685(32) σ, close to the middle of the lipid. For MARTINI-DMPC we find z0 = 0.850(11) nm, a value which is about 0.4 nm closer towards the bilayer midplane than the glycerol backbone. Applying the analysis of the previous section to the buckling simulations from Ref., 58 we find that the {xp , yp } pairs indeed form linear rela-

Cooke model The Cooke model only has three beads per lipid, and hence the function z0 (ξ) only has three


11


Table 2: Results for all simulations conducted in this work. To ease notation, we choose the Cooke-model generic parameters to be σ = 1 nm and kB T = 1.1 ǫ = 4.1 pN nm. Values with an asterisk have been previously determined in Ref., 58 those with a dagger in Ref. 60 Model

κm [kB T ]

z0 [nm]

ℓ [nm]

Cooke MARTINI Berger

6.4(2)∗ 14.5(5)∗ 12.4(5)∗

1.4685(32)† 0.850(11)† 1.3225(44)

1.32(22) 1.008(34) 1.609(43)

50

30 20 10 0 −10 0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

2πℓ 2 L

0.0156(51) 0.0184(12) 0.0665(36)

θmax (γ = 14 )

pN κt,m [ nm ]

0.85(28)◦ 1.008(66)◦ 3.64(20)◦

14.4±4.6 57.4±3.0 19.4±1.0

simulated strains γ. The result averaged over all simulations gives κt,m = (3.5 ± 1.1) kB T /σ 2 . With the mapping σ = 1 nm and kB T = 1.1 ǫ = 4.1 pN nm, this gives κt,m = (14.4±4.6) pN/nm. The error is fairly large, since the two available points do not constrain the slope of z0 (ξ) very precisely. The tilt decay length is found to be ℓ = 1.3 nm, and is hence clearly microscopic. For a strain of 25% the inferred maximum tilt angle is 0.85◦ and thus very small. Somewhat disconcertingly, the scatter of data points in Fig. 5 is substantially smaller than what the size of the error bars would suggest. We suspect that this originates from an overestimation of the errors of z0 from individual simulations, which are determined via Monte Carlo resampling of the raw data. One conceivable way how this could result in an overestimation is if some underlying systematic effect is mistaken for statistics. For instance, if the {xp , yp } plots result on average in a line, but with a systematic small “waviness” of integer period,† then this will not affect the precision with which the average slope can be determined, but the increased residuals will nevertheless suggest a higher variance. At this point we have not traced down the origin of this effect.

40 κt,m [pN nm−1 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 22

0.45

γ

Figure 5: Value of the tilt modulus κt,m for the Cooke model against strain γ. The shaded regions are the 1 σ and 2 σ confidence bands when fitting to a constant. points for any value of γ. Moreover, as we have found previously, 60 these three points do generally not lie on a common line—in the following sense: for each of the three points we can calculate an average and an error bar, and a linear fit is generally a statistically implausible description of these three points. However, a line drawn through the two points belonging to the two tail beads yields consistent values for the slope. This might sound trivial, since one can always draw a line through two points, but “consistency” here means that the slope is consistent for all values of the strain γ, a finding that would not have to hold, and which in fact does not hold if we were to fit a line through all three points. Fig. 5 shows the value of the tilt modulus κt inferred via Eqn. (43), plotted against all

MARTINI model The MARTINI model for DMPC contains 10 CG beads, and this offers a broader basis for †

Imagine, for instance, data points which scatter according to xp = sin(t) + δx sin(nx t) and yp = z0 [sin(t) + δy sin(ny t)] with t ∈ [0, 2π], δx , δy ≪ 1, and small integer coefficients nx and ny .


12

Page 13 of 22

24

80 75

22 κt,m [pN nm−1 ]

70 κt,m [pN nm−1 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


65 60 55 50 45 40 35

20 18 16 14

0.05

0.1

0.15

0.2

0.25 γ

0.3

0.35

0.4

0.45

0.1

0.15

0.2

0.25

0.3

0.35

γ

Figure 6: Value of the tilt modulus κt,m for the MARTINI model against strain γ. The filled symbols are determined from Eqn. (43); the open symbols are calculated without the strain 9 2 correction factor of 1 + 21 γ + 32 γ + · · · , showing that the latter is needed to arrive at a result that is no longer strain dependent. The shaded regions are the 1 σ and 2 σ confidence bands when fitting the filled data points to a constant.

Figure 7: Value of the tilt modulus κt,m for the Berger DMPC model against strain γ. The shaded regions are the 1 σ and 2 σ confidence bands when fitting to a constant. connecting tilt and shape and the replacement sin ψ → ψ in Eqn. (38) to ease including tilt in the lipid-imbalance equation. Again, the error bars in Fig. 6 appear too large if viewed in reference to the actual data scatter, and we again suspect that this is not due to a highly fortuitous absence of scatter but an overestimation of the error bars, but we have not traced down the origin of this phenomenon.

estimating the slope of z0 (ξ). Moreover, unlike in the Cooke model case, all points lie plausibly on a single line (see Fig. 8 in our previous publication 60 ). There are still small systematic effects, this time related to the fact that the inferred z0 values belonging to the beads of the sn1-chain are slightly smaller (by about 0.01 nm) than those inferred from the sn2-chain. Fig. 6 shows the values of the tilt modulus inferred from the slope in such graphs for all values of γ studied. Averaging the data, we find κt,m = (57.4 ± 3.0)pN/nm, and a tilt decay length of ℓ = 1.0 nm. The maximum inferred tilt angle at γ = 25% is found to be 1◦ , again very small. Unlike for the Cooke data, the MARTINI results for κt,m (γ) are precise enough to see that the γ-dependent correction in Eqn. (43) is crucial in order to get strain-independent values of the tilt modulus. This nontrivial correction illustrates that our theory describes the data well, despite several approximations, of which the two most important ones were the use of a simplified Euler Lagrange equation (26) for

Berger model While determining the pivotal plane position of the Berger model, we have seen that the function z0 (hξi i) is no longer a simple line—see again Fig. 4. This does not necessarily matter for a self-consistently determined z0 value, but for the tilt modulus we need the slope of z0 (hξi i), and this is now an ill-defined quantity, or at least one that depends on position. In the Discussion section we will reflect on some further implications of this result. For now, a pragmatic but plausible way forward is to observe that within the tail-region of the lipid the relation between z0 and the distance hξi i of a reference atom away from the bilayer’s midplane is still linear, and so we will define the slope of z0 (ξ) only with this subset of atoms. This means that, as far as tilt is concerned, we restrict our analysis to the tails.


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 22

tics. 60 The in principle curvature-dependent entropy embodied in the tail statistics is replaced by a curvature independent energy of a CG interaction potential. This reduces the extent to which tails can respond to curvature changes, which in turn pulls the pivotal plane closer to the leaflet’s midplane (i.e., the trivial position one would expect if the leaflet were elastically uniform along its entire width). If so, our explanation would also predict that the essentially atomistically resolved Berger force field should not suffer from this problem and instead place the pivotal plane closer to the glycerol backbone. This is indeed what it does. While it is dangerous to read too much into the outcome of a single investigation, the underlying CG approximation is well known and uncontested: CG potentials invariably inherit aspects of the state point at which they are constructed or for which they have been developed. Lipid models are generally optimized to reproduce the Lα phase of membranes, not highly curved lipid leaflets. Notice, however, that MARTINI is not a structurally coarsegrained model, and so it should suffer less from any all-too-eager attempt to optimize the structural aspects of a flat bilayer. And still, its very construction (equal bead size all the way down the tail) is at least inspired by the geometry of a bilayer phase. Had the developers aimed at surfactants that assemble into spherical micelles, a more “tapered” design might have suggested itself. The surprising lesson is that such subtle aspects in an otherwise fairly well-resolved model can show up so noticeably even in relatively benignly deformed bilayers. We hasten to add that this does of course does not invalidate the use of CG models to study membrane phenomena that involve membrane bending. The curvature modulus could be well reproduced, and this is often all that matters at larger scales. But we believe that CG studies which explicitly aim to explain the value of the bending modulus or its dependence on lipid architecture—say, how it depends on tail length and saturation—ought to be viewed with some skepticism. If the chain conformations in curved coarse-grained membranes do not— and can not—capture the physics of free fatty

Fig. 7 shows the values of the thus-inferred tilt modulus against all investigated values of the strain γ. Averaging over all strains, we get κt,m = (19.4 ± 1.0) pN/nm, which is about three times smaller than the value inferred from the CG MARTINI version of DMPC. In contrast, the tilt decay length is fairly similar, ℓ = 1.3 nm. Notice that the small value of the modulus also renders the maximum tilt angle bigger, we find θmax (γ = 41 ) = 3.67◦ , almost 4 times the maximum value for MARTINI under the same conditions.

Discussion Pivotal plane for Berger DMPC Before we go over the tilt moduli-measurements, let us first reflect on the pivotal plane position we measured for the united atom Berger model of DMPC, especially on the point that our result, z0,Berger = 1.3225(44) nm, is noticeably larger than our earlier MARTINI-derived value z0,MARTINI = 0.850(11) nm. A different way to state this is by comparing the position of the pivotal plane to that of the Luzzati plane, the membrane’s dividing surface towards the water phase, for which the phosphate position can be used as a simple proxy. Its distance from the bilayer’s midplane is in fact very similar for the two models: hξPh i = 1.74 nm for MARTINI and hξPh i = 1.80 nm for Berger. This leads to z0 /hξPh i = 0.49 for MARTINI, almost exactly in the middle of the leaflet, and z0 /hξPh i = 0.73 for Berger, substantially higher up towards the head region. Both general expectation and experiment suggest that z0 is just barely below the glycerol backbone region, where the lipid tails get joined to the more rigid head group region, giving a ratio z0 /hξPh i not in the middle of the bilayer but about 32 up towards the head group region. 80 Why would a coarse-grained but still fairly well resolved lipid model get that number wrong? In our previous paper we attributed this discrepancy to an invariable shortcoming of the lipid tails’ coarse-grained representation to capture the change of chain configurational statis-


14

Page 15 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


acid tails, and hence end up bending the bilayer around a noticeably shifted reference surface, the relative contributions of head- and tail-regions to the overall rigidity are misbalanced, rendering claims about the relation between bending modulus and specific tail properties somewhat suspect.

of the monolayer, κt = 2 κt,m ,

(46)

and this is the tilt modulus that the publications referenced above refer to. Strictly speaking, a consistent way to construct a bilayer Hamiltonian from two monolayer Hamiltonians is more subtle. For once, there is also a tilt sum field 12 (T↑ + T↓ ) which should enter the physics, and the leaflets need not follow each other perfectly (i.e., there might be breathing modes). But even in the absence of breathing, the translation of two leaflets which have their reference surfaces not at the bilayers midplane into one central reference surface gives rise to some interesting new terms, details of which we will discuss in another publication. Fortunately, for the present analysis none of these subtleties matter, and we may simply take the monolayer tilt moduli measured in this paper, multiply them by two, and arrive at effective bilayer tilt moduli. Watson et al. have determined the tilt modulus of Berger-DMPC at 303 K by measuring transverse lipid orientation fluctuations, finding κt = 56 pN/nm. 21 Using the same technique, but a CHARMM 36 model of DMPC, Venable et al. have found κt = 40.2 pN/nm. 23 Our own value of κt = (38.8 ± 2) pN/nm is of quite comparable magnitude. In fact, Venable et al. studied twelve common lipids and found values in the range 40 . . . 85 pN/nm. The only experimentally available tilt modulus so far has been determined by Jablin et al., 34 who find (95 ± 7) pN/nm, albeit for DOPC. Previous measurements of the tilt modulus of MARTINI DMPC do not yet appear to exist, but Watson et al. have measured the tilt modulus of MARTINI DPPC, 69 which differs by having one extra CG bead per chain. They find κt = 110 pN/nm, about a factor of 2 larger than their Berger DMPC value. Our own value for MARTINI DMPC, κt = (115 ± 6) pN/nm is about a factor 3 larger that the united atom Berger value. Both measurements suggest that MARTINI lipids tilt less easily than the more highly resolved united atom models. In their investigation of bending and tilt

The tilt moduli The tilt moduli of several different atomistic lipid models have recently been determined in simulations, relying on a spectral analysis of lipid orientation in a (possibly quite small) flat bilayer. 22,23 Also, a first experimental measurement has been presented, which analyzed membrane undulations at large wave vectors. 34 However, in both cases the authors quote an effective bilayer tilt modulus, which is the modulus that couples to a tilt that is defined as a suitable average of the tilt fields in the two apposing monolayers. With the sign conventions from Watson et al., 21 we can define a bilayer tilt field as 1 T (b) = (44) T↑ − T↓ , 2 where T↑ and T↓ are the tilt fields in the upper and lower monolayer, respectively. The minus sign derives from the fact that if the tilted lipid orientations are collinear across the two leaflets (as in Fig. 2), then the tangential tilt field points in the opposite direction in the lower leaflet. The simplest possible tilt theory one can now assume for the bilayer has a form that mimics Eqn. (4): Z 2 1 H = dA κ K + ∇ · T (b) 2 1 (b) 2 ˜ (45) +¯ κKG + κt (T ) , 2 where the moduli now refer to the bilayer (and the spontaneous bilayer curvature vanishes for symmetry reasons). Since tilting two leaflets costs twice the energy as tilting a single leaflet, the modulus for the bilayer has to be twice that


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

moduli through lipid orientation fluctuations, Levine et al. point out that the value for the tilt modulus may shift by up to ±20% depending on how they define the lipid normal. 22 In the supplementary material to their publication they compare three different definitions. While they all use the midpoint between the terminal methyls of the two tails as one anchor point, they differ in the choice of the other anchor point: either the phosphate atom, or the glycerol C2 atom, or the midpoint between those two (the latter is the authors’ preferred choice). They find that the inferred tilt modulus increases as one moves the upper reference point further out. Our own studies suggest an explanation of this finding. Recall from Fig. 4 that z0 (ξ) is not simply a straight line, the shape which our theory predicts under the assumption that the lipids tilt uniformly. Instead, the slope decreases beyond the glycerol backbone and even becomes negative. This clearly shows that the head of the lipid does not participate in tilting to the same extent as the tails do. Had we described the lipid by a simple line between two reference points, then the resulting slope in z0 (ξ) would be smaller the higher up the headgroup choice is taken. But from Eqn. (43), smaller slopes result in bigger moduli, and so we recover the same trend as Levine et al. do. The non-uniform tilting of lipids is more than a nuisance for defining the orientation. It reminds us that the notion of lipid tilt is more subtle than simply a rotation of the lipid without a change of its shape. Interestingly, our own analysis never requires us to define an orientation vector, and so we can empirically test whether or not lipids tilt uniformly— concluding that they do not. Moreover, the data suggest that the tail region of lipids shows a clean linear trend in z0 (ξ), suggesting that the tails do tilt uniformly, while the head group shows a more complex behavior. One might be tempted to define a local tilt modulus along the height of the lipid, via the local slope of z0 (ξ), but notice that this would give an infinite modulus at the point where z0 (ξ) has its maximum (roughly at ξ = 1.6 nm, where the C3 atom of the glycerol backbone sits). We instead prefer to define a tilt modulus for the tail region alone

Page 16 of 22

(since it is well-defined). This also means that our effective lipid orientation is more tilted than the one which Levine et al. infer from their preferred choice, 22 which is one possible reason for why our tilt modulus is slightly smaller than theirs. Let us finally address the question why we do not determine the tilt modulus more directly. Eqn. (26) clearly states that κt,m = κm

K ′ (s) , T (s)

(47)

suggesting that we could measure κt,m by separately measuring tilt and curvature gradients. Unfortunately, this is an extremely noisy measurement. It is difficult enough to measure the tilt—for which the hard part is not to find the lipid orientation n, but to fit the membrane shape and determine the direction of the local membrane normal N , so that we can calculate T essentially as the difference between those two—see Eqn. (2). The latter ends up calculating a small number (recall that lipids tilt by not more than a few degrees) as the difference between two large numbers (at γ ≈ 0.543 we have ψmax = 90◦ ). But even worse, the curvature gradient is K ′ = −ψ ′′ (s), showing that we need to differentiate the angle two more times to get to the second variable in Eqn. (47), and this deteriorates the precision. These noise problems are compounded by the fact that we cannot trivially average over different simulation frames, since bilayers fluctuate strongly enough that every shape needs to be individually fitted, as previously described. 60 In contrast, the more indirect method we propose here never needs to calculate more than angles, resulting in manageable noise—especially since the outputs from individual frames, the slopes in the {xp , yp } plots, have been fully disentangled from the individual buckle they derive from, so that they can be trivially averaged. All this shows that while our proposed method might appear somewhat circuitous, it thereby avoids noisy variables and hence is more accurate than the seemingly more direct way. At this point we do not suggest that our method is necessarily more efficient or more


16

Page 17 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


low the glycerol backbone. This supports the idea that lipid head and lipid tail differ noticeably in their ability to respond to local lateral stress. The tails are more pliable, and so the pivotal plane divides a bilayer leaflet in an uneven way, with the larger fraction devoted to the soft tails and the smaller fraction to the stiffer heads, thereby lowering the bending energy compared to an even division. In turn, our previous observation 60 of a very even division in MARTINI DMPC suggests that this difference in stiffness is not well captured on the coarse-grained level, suggesting that attempts to explain on the CG level how changes in lipid topology affect a membrane’s bending rigidity might not be very reliable.

precise than traditional fluctuation approaches, especially since very small systems appear to suffice to get power spectra of orientation fluctuations over the necessary range of wave vectors. 21,22 However, given the persistent difficulties with pinpointing elastic moduli of membranes, 81 it seems advantageous to have alternative means available for measuring the same observable. This not only enables independent measurements that might suffer from different statistics or systematics, but it also permits one to test how well the microscopic theories on which all analysis techniques ultimately rely really work.

Conclusions We have shown that our previously proposed method for extracting the pivotal plane of a lipid membrane can also give access to the lipid’s tilt modulus at no additional cost. We have applied our analysis to three coarse grained models of very different resolution and showed that it works well in all three cases. For MARTINI DMPC and Berger DMPC we find values for the moduli which are compatible (within the scatter still common for simulations of moduli) with other simulations that use very different techniques relying on a fluctuation analysis. For the highly coarse-grained Cooke model no numbers for comparison existed, but the results for maximum tilt angle and tilt decay length are very similar to the other two models, especially MARTINI. We find that MARTINI DMPC is harder to tilt than a more finely resolved united atom model of DMPC. Moreover, analysis of the latter model suggests that these lipids do not tilt uniformly, only the tails do. Hence, the definition of tilt becomes less obvious, as previously noted, 22 but our analysis strategy does not have to commit to a specific definition of the tilt direction, since it instead backs out the extent of local tilt from the data. Finally, our new determination of the pivotal plane for Berger-DMPC shows that, in contrast to the CG MARTINI DMPC counterpart, the pivotal plane indeed lies very slightly be-

Acknowledgement We thank Mert Terzi for stimulating discussions. Partial support from NSF grant CHE-1464926 is also gratefully acknowledged.

References (1) Canham, P. B. The Minimum Energy of Bending as a Possible Explanation of the Biconcave Shape of the Human Red Blood Cell. J. Theoret. Biol. 1970, 26, 61–81. (2) Helfrich, W. Elastic Properties of Lipid Bilayers—Theory and Possible Experiments. Z. Naturforsch. C 1973, 28, 693– 703. (3) Evans, E. A. Bending Resistance and Chemically-Induced Moments in Membrane Bilayers. Biophys. J. 1974, 14, 923– 931. (4) Gompper, G.; Klein, S. Ginzburg-Landau Theory of Aqueous Surfactant Solutions. J. Phys. II (France) 1992, 2, 1725–1744. (5) Nelson, P.; Powers, T. Rigid Chiral Membranes. Phys. Rev. Lett. 1992, 69, 3409– 3412.


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(6) Nelson, P. N.; Powers, T. Renormalization of Chiral Couplings in Tilted Bilayer Membranes. J. Phys. II (France) 1993, 3, 1535–1569.

Page 18 of 22

(16) Hamm, M.; Kozlov, M. M. Elastic Energy of Tilt and Bending of Fluid Membranes. Eur. Phys. J. E 2000, 3, 323–335.

(7) Powers, T.; Nelson, P. Fluctuating Membranes with Tilt Order. J. Phys. II (France) 1995, 5, 1671–1678.

(17) May, E. R.; Narang, A.; Kopelevich, D. I. Role of Molecular Tilt in Thermal Fluctuations of Lipid Membranes. Phys. Rev. E 2007, 76, 021913.

(8) Seifert, U.; Shillcock, J.; Nelson, P. Role of Bilayer Tilt Difference in Equilibrium Membrane Shapes. Phys. Rev. Lett. 1996, 77, 5237–5240.

(18) May, E. R.; Narang, A.; Kopelevich, D. I. Molecular Modeling of Key Elastic Properties for Inhomogeneous Lipid Bilayers. Molecular Simulation 2007, 33, 787–797.

(9) Ben Shaul, A.; Szleifer, I.; Gelbart, W. M. Chain Organization and Thermodynamics in Micelles and Bilayers. 1. Theory. J. Chem. Phys. 1985, 83, 3597–3611.

(19) Watson, M. C.; Brown, F. L. Interpreting Membrane Scattering Experiments at the Mesoscale: The Contribution of Dissipation within the Bilayer. Biophys. J. 2010, 98, L9–L11.

(10) Szleifer, I.; Ben Shaul, A.; Gelbart, W. M. Chain Organization and Thermodynamics in Micelles and Bilayers. 2. Model Calculations. J. Chem. Phys. 1985, 83, 3612– 3620.

(20) Watson, M. C.; Penev, E. S.; Welch, P. M.; Brown, F. L. H. Thermal Fluctuations in Shape, Thickness, and Molecular Orientation in Lipid Bilayers. J. Chem. Phys. 2011, 135, 244701.

(11) Szleifer, I.; Kramer, D.; Ben-Shaul, A.; Roux, D.; Gelbart, W. M. Curvature Elasticity of Pure and Mixed Surfactant Films. Phys. Rev. Lett. 1988, 60, 1966–1969.

(21) Watson, M. C.; Brandt, E. G.; Welch, P. M.; Brown, F. L. H. Determining Biomembrane Bending Rigidities from Simulations of Modest Size. Phys. Rev. Lett. 2012, 109, 028102.

(12) Szleifer, I.; Kramer, D.; Ben-Shaul, A.; Gelbart, W. M.; Safran, S. A. Molecular Theory of Curvature Elasticity in Surfactant Films. J. Chem. Phys. 1990, 92, 6800–6817.

(22) Levine, Z. A.; Venable, R. M.; Watson, M. C.; Lerner, M. G.; Shea, J.-E.; Pastor, R. W.; Brown, F. L. H. Determination of Biomembrane Bending Moduli in Fully Atomistic Simulations. J. Am. Chem. Soc. 2014, 136, 13582–13585.

(13) Rawicz, W.; Olbrich, K.; McIntosh, T.; Needham, D.; Evans, E. Effect of Chain Length and Unsaturation on Elasticity of Lipid Bilayers. Biophys. J. 2000, 79, 328– 339.

(23) Venable, R. M.; Brown, F. L.; Pastor, R. W. Mechanical Properties of Lipid Bilayers from Molecular Dynamics Simulation. Chem. Phys. Lipids 2015, 192, 60– 74.

(14) Uline, M.; Szleifer, I. Mode Specific Elastic Constants for the Gel, LiquidOrdered, and Liquid-Disordered Phases of DPPC/DOPC/Cholesterol Model Lipid Bilayers. Faraday Discuss. 2013, 161, 177–191.

(24) Brochard, F.; Lennon, J. F. Frequency Spectrum of Flicker Phenomenon in Erythrocytes. J. de Physique 1975, 36, 1035– 1047.

(15) Hamm, M.; Kozlov, M. Tilt Model of Inverted Amphiphilic Mesophases. Eur. Phys. J. B 1998, 6, 519–528.

(25) Brochard, F.; de Gennes, P.-G.; Pfeuty, P. Surface-Tension and Deformations of


18

Page 19 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Theory of Biomembrane Mechanics. Phys. Rev. Lett. 2014, 113, 248102.

Membrane Structures—Relation to 2Dimensional Phase-Transitions. J. de Physique 1976, 37, 1099–1104.

(35) Goetz, R.; Gompper, G.; Lipowsky, R. Mobility and Elasticity of Self-Assembled Membranes. Phys. Rev. Lett. 1999, 82, 221–224.

(26) Schneider, M. B.; Jenkins, J. T.; Webb, W. W. Thermal Fluctuations of Large Cylindrical Phospholipid-Vesicles. Biophys. J. 1984, 45, 891–899.

(36) Lindahl, E.; Edholm, O. Mesoscopic Undulations and Thickness Fluctuations in Lipid Bilayers from Molecular Dynamics Simulations. Biophys. J. 2000, 79, 426– 433.

(27) Schneider, M. B.; Jenkins, J. T.; Webb, W. W. Thermal Fluctuations of Large Quasi-Spherical Bimolecular Phospholipid-Vesicles. Biophys. J. 1984, 45, 1457–1472.

(37) Marrink, S. J.; Mark, A. E. Effect of Undulations on Surface Tension in Simulated Bilayers. J. Pjys. Chem. B 2001, 105, 6122–6127.

(28) Faucon, J. F.; Mitov, M. D.; Méléard, P.; Bivas, I.; Bothorel, P. Bending Elasticity and Thermal Fluctuations of LipidMembranes—Theoretical and Experimental Requirements. J. de Physique 1989, 50, 2389–2414.

(38) Ayton, G.; Voth, G. A. Bridging Microscopic and Mesoscopic Simulations of Lipid Bilayers. Biophys. J. 2002, 83, 3357–3370.

(29) Henriksen, J.; Rowat, A. C.; Ipsen, J. H. Vesicle Fluctuation Analysis of the Effects of Sterols on Membrane Bending Rigidity. Eur. Biophys. J. 2004, 33, 732–741.

(39) Hofsäß, C.; Lindahl, E.; Edholm, O. Molecular Dynamics Simulations of Phospholipid Bilayers with Cholesterol. Biophys. J. 2003, 84, 2192–2206.

(30) Liu, Y. F.; Nagle, J. F. Diffuse Scattering Provides Material Parameters and Electron Density Profiles of Biomembranes. Phys. Rev. E 2004, 69, 040901.

(40) Farago, O. “Water-free” Computer Model for Fluid Bilayer Membranes. J. Chem. Phys. 2003, 119, 596–605.

(31) Pan, J.; Tristram-Nagle, S.; Ku¸cerka, S.; Nagle, J. F. Temperature Dependence of Structure, Bending Rigidity, and Bilayer Interactions of Dioleoylphosphatidylcholine Bilayers. Biophys. J. 2008, 94, 117–124.

(41) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750–760. (42) Brannigan, G.; Tamboli, A. C.; Brown, F. L. H. The Role of Molecular Shape in Bilayer Elasticity and Phase Behavior. J. Chem. Phys. 2004, 121, 3259–3271.

(32) Pan, J.; Mills, T. T.; Tristram-Nagle, S.; Nagle, J. F. Cholesterol Perturbs Lipid Bilayers Nonuniversally. Phys. Rev. Lett. 2008, 100, 198103.

(43) Wang, Z.-J.; Frenkel, D. Modeling Flexible Amphiphilic Bilayers: A Solvent-Free Off-Lattice Monte Carlo Study. J. Chem. Phys. 2005, 122, 234711.

(33) Pan, J.; Tristram-Nagle, S.; Nagle, J. F. Effect of Cholesterol on Structural and Mechanical Properties of Membranes Depends on Lipid Chain Saturation. Phys. Rev. E 2009, 80, 021931.

(44) Cooke, I. R.; Kremer, K.; Deserno, M. Tunable Generic Model for Fluid Bilayer Membranes. Phys. Rev. E 2005, 72, 011506.

(34) Jablin, M. S.; Akabori, K.; Nagle, J. F. Experimental Support for Tilt-Dependent


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 22

(54) Harmandaris, V. A.; Deserno, M. A Novel Method for Measuring the Bending Rigidity of Model Lipid Membranes by Simulating Tethers. J. Chem. Phys. 2006, 125, 204905.

(45) Cooke, I. R.; Deserno, M. Solvent-Free Model for Self-Assembling Fluid Bilayer Membranes: Stabilization of the Fluid Phase Based on Broad Attractive Tail Potentials. J. Chem. Phys. 2005, 123, 224710.

(55) Arkhipov, A.; Yin, Y.; Schulten, K. FourScale Description of Membrane Sculpting by BAR Domains. Biophys. J. 2008, 95, 2806–2821.

(46) Brannigan, G.; Philips, P. F.; Brown, F. L. H. Flexible Lipid Bilayers in Implicit Solvent. Phys. Rev. E 2005, 72, 011915.

(56) Baoukina, S.; Marrink, S. J.; Tieleman, D. P. Molecular Structure of Membrane Tethers. Biophys, J. 2012, 102, 1866–1871.

(47) Wang, Z.-J.; Deserno, M. A Systematically Coarse-Grained Solvent-Free Model for Quantitative Phospholipid Bilayer Simulation. J. Phys. Chem. B 2010, 114, 11207–11220.

(57) Noguchi, H. Anisotropic Surface Tension of Buckled Fluid Membranes. Phys. Rev. E 2011, 83, 061919.

(48) Wang, Z.-J.; Deserno, M. Systematic Implicit Solvent Coarse-Graining of Bilayer Membranes: Lipid and Phase Transferability of the Force Field. New J. Phys. 2010, 12, 095004.

(58) Hu, M.; Diggins, P.; Deserno, M. Determining the Bending Modulus of a Lipid Membrane by Simulating Buckling. J. Chem. Phys. 2013, 138, 214110.

(49) Brandt, E. G.; Braun, A. R.; Sachs, J. N.; Nagle, J. F.; Edholm, O. Interpretation of Fluctuation Spectra in Lipid Bilayer Simulations. Biophys. J. 2011, 100, 2104– 2111.

(59) Diggins, P.; McDargh, Z. A.; Deserno, M. Curvature Softening and Negative Compressibility of Gel-Phase Lipid Membranes. J. Am. Chem. Soc. 2015, 137, 12752–12755.

(50) Shiba, H.; Noguchi, H. Estimation of the Bending Rigidity and Spontaneous Curvature of Fluid Membranes in Simulations. Phys. Rev. E 2011, 84, 031926.

(60) Wang, X.; Deserno, M. Determining the Pivotal Plane of Fluid Lipid Membranes in Simulations. J. Chem. Phys. 2015, 143, 164109.

(51) Bo, L.; Waugh, R. E. Determination of Bilayer Membrane Bending Stiffness by Tether Formation from Giant, ThinWalled Vesicles. Biophys. J. 1989, 55, 509–517.

(61) Kreyszig, E. Differential Dover: New York, 1991.

Geometry;

(62) Levien, R. The Elastica: A Mathematical History; 2008.

(52) Cuvelier, D.; Derényi, I.; Bassereau, P.; Nassoy, P. Coalescence of Membrane Tethers: Experiments, Theory, and Applications. Biophys. J. 2005, 88, 2714–2726.

(63) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, 9th ed.; Dover: New York, 1972.

(53) Tian, A.; Baumgart, T. Sorting of Lipids and Proteins in Membrane Curvature Gradients. Biophys. J. 2008, 96, 2676– 2688.

(64) Dimova, R.; Pouligny, B.; Dietrich, C. Pretransitional Effects in Dimyristoylphosphatidylcholine Vesicle Membranes: Optical Dynamometry Study. Biophys. J. 2000, 79, 340–356.


20

Page 21 of 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(65) Lee, C.-H.; Lin, W.-C.; Wang, J. All-Optical Measurements of the Bending Rigidity of Lipid-Vesicle Membranes Across Structural Phase Transitions. Phys. Rev. E 2001, 64, 020901.

(74) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690.

(66) Steltenkamp, S.; M¨ uller, M. M.; Deserno, M.; Hennesthal, C.; Steinem, C.; Janshoff, A. Mechanical Properties of Pore-Spanning Lipid Bilayers Probed by Atomic Force Microscopy. Biophys. J. 2006, 91, 217–226.

(75) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (76) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697.

(67) Deserno, M. Fluid Lipid Membranes: From Differential Geometry to Curvature Stresses. Chem. Phys. Lipids 2015, 185, 11–45.

(77) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190.

(68) Deserno, M. Mesoscopic Membrane Physics: Concepts, Simulations, and Selected Applications. Macromol. Rapid Comm. 2009, 30, 752–771.

(78) Nosé, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055–1076. (79) Rzepiela, A. J.; Schäfer, L. V.; Goga, N.; Risselada, H. J.; De Vries, A. H.; Marrink, S. J. Reconstruction of Atomistic Details from Coarse-Grained Structures. J. Comput. Chem. 2010, 31, 1333–1343.

(69) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812– 7824.

(80) Rand, R. P.; Fuller, N. L.; Gruner, S. M.; Parsegian, V. A. Membrane Curvature, Lipid Segregation, and Structural Transitions for Phospholipids under DualSolvent Stress. Biochem. 1990, 29, 76–87.

(70) Marrink, S. J.; Tieleman, D. P. Perspective on the MARTINI Model. Chem. Soc. Rev. 2013, 42, 6801–6822. (71) Berger, O.; Edholm, O.; Jähnig, F. Molecular Dynamics Simulations of a Fluid Bilayer of Dipalmitoylphosphatidylcholine at Full Hydration, Constant Pressure, and Constant Temperature. Biophys. J. 1997, 72, 2002–2013.

(81) Nagle, J. F.; Jablin, M. S.; TristramNagle, S.; Akabori, K. What are the True Values of the Bending Modulus of Simple Lipid Bilayers? Chem. Phys. Lipids 2015, 185, 3–10.

(72) Limbach, H.-J.; Arnold, A.; Mann, B. A.; Holm, C. ESPResSo – An Extensible Simulation Package for Research on Soft Matter Systems. Comput. Phys. Commun. 2006, 174, 704–727. (73) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. ACS Paragon Plus Environment

21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Buckling a lipid membrane induces the lipids in both leaflets to tilt, and comparing the extent of tilt with the prediction from continuum theory leads to the tilt modulus. We use this technique to measure the tilt modulus for three lipid models—Cook, MARTINI, and Berger—which differ in their chemical resolution. 61x39mm (300 x 300 DPI)


Page 22 of 22

Determining the Lipid Tilt Modulus by Simulating Membrane Buckles

Recommend Documents