A Computational Approach for Modeling Neutron ... - ACS Publications

Subscriber access provided by University of Newcastle, Australia

Article

A Computational Approach for Modeling Neutron Scattering Data from Lipid Bilayers Jan-Michael Y. Carrillo, John Katsaras, Bobby G. Sumpter, and Rana Ashkar J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00968 • Publication Date (Web): 12 Jan 2017 Downloaded from http://pubs.acs.org on January 13, 2017

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.

Journal of Chemical Theory and Computation 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 31

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

Journal of Chemical Theory and Computation

A Computational Approach for Modeling Neutron Scattering Data from Lipid Bilayers Jan-Michael Y. Carrillo,∗,†,‡ John Katsaras,¶,§,k Bobby G. Sumpter,†,‡ and Rana Ashkar∗,¶,k †Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States ‡Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States ¶Biology & Soft Matter Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States §Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, United States kShull Wollan Center: a Joint Institute for Neutron Science, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment


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

Abstract Biological cell membranes are responsible for a range of structural and dynamical phenomena crucial to a cell’s well-being and its associated functions. Due to the complexity of cell membranes, lipid bilayer systems are often used as biomimetic models. These systems have led to significant insights into vital membrane phenomena such as domain formation, passive permeation and protein insertion. Experimental observations of membrane structure and dynamics are, however, limited in resolution, both spatial and temporal. Importantly, computer simulations are starting to play a more prominent role in interpreting experimental results, enabling a molecular understanding of lipid membranes. In particular, the synergy between scattering experiments and simulations offers opportunities for new discoveries in membrane physics, as the length and time scales probed by molecular dynamics (MD) simulations parallel those of experiments. Here, we describe a coarse-grained MD simulation approach that mimics neutron scattering data from large unilamellar lipid vesicles over a range of bilayer rigidities. Specifically, we simulate vesicle form factors and membrane thickness fluctuations determined from small angle neutron scattering (SANS) and neutron spin echo (NSE) experiments, respectively. Our simulations accurately reproduce trends from experiments and lay the groundwork for studies of more complex membrane systems.

1

INTRODUCTION

The lipid bilayer is the underlying structure in biological cell membranes. Membranes are responsible for a range of functions, some of which include the regulation of different protein functions, and the exchange of nutrients in and out of the cell. Such processes take place over a range of length and time scales; i.e., from membrane ionic transport on the sub-nanoscale, to vesicular budding and cell motility on the micron scale. Research on the physical properties of lipid membranes has shown that many of these phenomena, including the transport of small molecules and the insertion of membrane proteins, are strongly de-

2


Page 2 of 31

Page 3 of 31

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


pendent on lipid-membrane dynamics. 1 The expectation is that the length and time scales over which these phenomena occur ought to be coupled to the motions of the lipid bilayer matrix, necessitating a full understanding of the hierarchy of dynamical modes present in lipid bilayers. Commonly used experimental techniques are capable of probing the fastest and slowest regimes known to exist in membranes. For example, picosecond conformational transitions in lipid acyl tails can be directly examined by NMR and other spectroscopic techniques, 2 whereas micron-scale membrane reshaping can be easily accessed by optical microscopy. 3 Computationally, atomistic molecular dynamics (MD) simulations have been used to interpret low resolution neutron scattering data. 4,5 In addition, atomistic MD simulations have contributed significant insights into fast lipid dynamic processes using small bilayer patches (computationally efficient), bolstering the NMR and fast diffusion measurements. 6 The mesoscopic or intermediate dynamic regime, on the other hand, has proven to be far more challenging to access both experimentally and computationally. However, the emergence of new experimental methods and recent advances in computational capabilities have started to shed new light on important mesoscale membrane phenomena that are intimately linked to collective membrane motions such as protein insertion, 7 lateral phase separation, 8 and protein localization during vesiculation 9 or sporulation. 10 The primary collective dynamic modes in the mesoscopic regime are thermally-excited undulations and peristaltic motions, usually referred to as bending fluctuations and thickness fluctuations, respectively. 11 The theoretical realization of collective motions in lipid membranes dates back to the early 70’s, namely to Helfrich’s continuum theory of structureless elastic membranes. 12 In it, the equilibrium properties of undulating membranes are determined by the energetic penalties of membrane motions incurred from bending and stretching, for example. However, Helfrich’s theory fails to provide a detailed description of the molecular arrangement of the lipids within bilayers, and the effect it has on membrane properties and dynamics. Fortunately, the development of MD simulations over the last few decades has allowed us

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

to bridge the gap between atomistic and collective membrane phenomena. The possibility of varying the extent of atomic grouping, or coarse-graining, in simulations has also enabled virtual measurements of membrane dynamics to determine specific membrane properties at an affordable computational cost. For instance, modeling mesoscale fluctuations requires systems that span length and time scales of hundreds of nm and ns, respectively. Although these conditions are extremely challenging for all-atom simulations, they are well-suited for coarse-grained MD simulations. Indeed, coarse-grained (CG) simulations have proven to be tremendously useful in predicting mesoscale membrane dynamics 11,13 and reproducing experimental results. 14–16 However, in order to stabilize the membrane, most of these models include explicit solvent molecules, which generally account for the majority of the simulated particles and which in turn account for most of the computational cost. Recently, Cooke et al. 17,18 suggested a solvent-free model which is far less-computationally expensive than the standard CG models of lipid membranes, allowing for simulations of even larger membrane systems that better represent those studied experimentally — i.e., microscale flat bilayers and large unilamellar vesicles. Experimentally, membrane bending fluctuations have been studied by a number of different techniques including fluctuation spectroscopy, scattering, and mechanical deformation measurements. 19 Membrane thickness fluctuations, on the other hand, have proven experimentally more challenging since they can only be observed over a narrow window of length and time scales; i.e. length scales that are comparable to the membrane thickness (few nm) and timescales of a few ns to a hundreds of 100 ns. In fact, the first experimental observation of membrane thickness fluctuations was reported recently in the pioneering neutron spin echo (NSE) measurements by Farago and colleagues. 20,21 Later, Nagao and coworkers refined this approach to study thickness fluctuations in oil-swollen surfactant bilayers, 22 as well as single-component 23 and binary 24 lipid vesicles. To our knowledge, NSE spectroscopy is the only technique that has the required patial and temporal resolution for measuring thickness fluctuations in lipid membranes. NSE has the added advantage of simultaneously access-

4


Page 4 of 31

Page 5 of 31

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


ing and resolving both bending and thickness fluctuation modes. For such measurements, the NSE instrument is configured to detect length scales that correspond to the membrane thickness. Thus, NSE measurements are usually coupled with small angle neutron scattering (SANS) for structural characterization of the membranes, usually in the form of vesicles, prior to experimentation. Here we demonstrate the feasibility of CG MD simulations in reproducing neutron scattering results from large unilamellar lipid vesicles. We specifically focus on SANS and NSE measurements, which can access unique structural and dynamical lipid membrane information on length and time scales compatible with simulations. The present simulations are built on a coarse-grained description of lipid molecules 17,18 used by Reynwar et al.. 9 We describe the MD simulation protocol that we used to reproduce the scattering data, and present the method used for obtaining the energy spectra. Specifically, we use a discrete Fourier transform that grows on the order of O(ng log ng ) as a function of the number of grid points ng . In contrast to the direct calculation of the spectra, which grows as O(n2p ), this method is appealing when modeling systems with features spanning ∼ 0.1 nm to 1 µm, and comprised of a large number of particles, np . We then compare the spectra obtained from simulations with known methods for extracting information from scattering experiments, and self consistently check the results with the corresponding real-space calculations obtained from particle trajectories. We want to emphasize that the aim of this manuscript is to validate the efficacy of our computational approach in reproducing available scattering patterns and to set the stage for its later use as an analysis tool for structural and dynamical scattering studies of lipid membranes. The results show that the CG model discussed here is ideal for simulating real size membranes measured experimentally and that are currently challenging to model via allatom simulations. We also found that the model is highly efficient in computing membrane structures and collective dynamics observed in scattering measurements, and can be reliably used to simulate more complex membranes. In the future, this computational approach will

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

Page 6 of 31

be extended to incorporate inputs from finer grain simulations (e.g., atomistic or united atom models), in consonance with the “computational microscopy” concept described in recent reviews by Marrink et al. 25 and Schulten et al. 26

2 2.1

METHODS Coarse-Grained Lipid Model

We use a coarse-grained model for lipid bilayers based on the solvent-free approach developed by Cooke et al. 17,18 In this model, each lipid molecule is represented by three connected beads — a “head” and two “tail” beads — with the solvent being implicitly treated. The pairwise interaction potential between beads consists of a purely repulsive potential for head-head and head-tail interactions, and an attractive potential for tail-tail interactions. The repulsive portion is described by the Weeks-Chandler-Andersen (WCA) potential, and is given by,

UW CA (r) =

   4ε (b/r)12 − (b/r)6 + 1 , r ≤ rc 4   0,

(1)

r > rc ,

where r is the radial distance between the centers of two beads in units of σ, rc is the cutoff given by rc = 21/6 b, ε = kB T ( where kB T is thermal energy), and b is respectively set to 0.95 σ, 0.95 σ, and 1.0 σ for head-to-head, head-to-tail and tail-to-tail distances for which the inter-particle potential is zero at r > rc . The attractive portion of the pair-wise potential describing tail-tail interactions is given by,     −ε, r < rc     Ucos (r) = −ε cos2 [π(r − rc )/2wc ], rc ≤ r ≤ rc + wc       0, r > rc + w c

(2)

which has a well depth of ε = kB T , and for r > rc smoothly approaches zero with a decay 6


Page 7 of 31

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


of wc . Cooke et al. showed that wc = 1.3 σ is at the boundary of the unstable phase, while wc = 1.7 σ is at the boundary of the gel phase. 18 In this study we constrain wc in a range from 1.3 σ to 1.7 σ for fluid membrane systems. The three beads are connected by two FENE bonds described by, " 2 # 1 r UF EN E (r) = − kbond ro2 ln 1 − 2 ro

(3)

with kbond = 30 kB T /σ 2 and ro = 1.5 σ. The beads are straightened by applying a bending potential

Ubend (θ) = kbend [1 + cos(θ)] ,

(4)

where θ is the angle between the unit vectors defined by the middle tail-bead to head-bead and the middle tail-bead to end tail-bead. The minimum of eq.4 is at θ = π rad, which favors a colinear configuration among the 3 lipid beads. A kbend of 20 kB T was used, which results in an approximate lipid persistence length of 19 σ. 27

2.2

Planar Bilayer

Preassembled planar bilayer membranes consisting of N = 18720 or N = 168480 lipids arranged in a hexagonal lattice with an approximate area per lipid of 1 σ 2 (see Figure 1(a)) were equilibrated using MD simulations. The simulation box was subjected to periodic boundary conditions in the x, y and z directions. Equilibration was performed using the LAMMPS 28,29 software, where the equations of motion were integrated at a time step of 0.01 τ , and the lipids were in an isobaric-isothermal (NPT) ensemble so as to render a laterally tensionless bilayer (in the xy-plane). A spring force with restoring energy Uspring,z = 1 k z2 2 s cm

and ks = 2000 kB T /σ 2 was applied to all beads, pegging the z coordinate of the

center-of-mass, zcm , to zero. In the equilibration simulation, the temperature was controlled by a Nosé-Hoover thermostat 30 with a damping parameter of 1 τ , while the pressures in 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

the x and y directions were independently controlled by a Nosé-Hoover barostat 30 with a damping constant of 5 τ . The equilibration run proceeded for up to 5 × 103 τ . The size of the x and y dimensions of the simulation box (Lx and Ly ) saturated at ≈ 3×103 τ . The final 1 × 103 τ frames of the simulation were used for averaging, thus enabling one to determine the size of the simulation box in the succeeding production runs. During the production runs, the x and y dimensions of the simulation box were adjusted to match the average values of the NPT ensemble run for a time interval of 50 τ , after which the canonical ensemble (NVT) production run proceeded for up to 2 × 104 τ . Data with elapsed times greater than 1×104 τ were used for data analysis. Concurrently, we performed 5 independent production simulation runs for different wc of the system with 18720 lipids. The thermostats for these simulation runs were replaced by a Langevin thermostat 31 with a friction coefficient of τ −1 , matching the original parameters used by Cooke et al. 17,18 The results of these simulations were then compared with the Helfrich theory 12 to extract the estimates of the bending modulus κ, presented in the “Results and Discussion” section.

2.3

Vesicles

Simulations on vesicles were performed using a protocol analogous to that for planar bilayer systems. The initial vesicle configuration consisted of N = 18996 lipids preassembled with an approximate radius of 30 σ and an area per lipid of 1 σ 2 for both inner and outer leaflets (see Figure 1(b)). NVT simulations of lipid systems having wc parameters ranging from 1.3 σ to 1.7 σ were performed for a duration of 5 × 104 τ . The temperature was controlled by a Langevin thermostat 31 with a friction coefficient of τ −1 . A spring force with energy defined 2 as Uspring,r = 12 ks rcm with ks = 2000 kB T /σ 2 was applied to all beads, fixing the vesicle

center-of-mass, rcm , to coordinates {0, 0, 0}. The application of a restoring spring introduces effects similar to thermal fluctuations or protrusions of individual lipid molecules but has little effect on the collective dynamics or structural properties, such as lipid bilayer thickness and fluctuations. The simulated systems equilibrated at ≈ 1 × 103 τ for the less rigid 8


Page 8 of 31

Page 9 of 31

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


bilayers, which is relatively fast because of the high flip-flop rate in the Cooke model. 17,18 It is therefore, reasonable to assume that the number ratio between the inner and outer leaflets is at equilibrium. The equilibrium structure is then reached, as indicated by the saturation of the average value of the outer-leaflet radius, hRo i shown in Figure 2. Data analysis was performed after an elapsed time greater than 3 × 104 τ , or when structural properties such as hRo i reached equilibrium. The equilibrium structure of the vesicle and the collective motion of the lipids in the simulations were analyzed in reciprocal space, simulating typical SANS and NSE experiments. The collective motions describing the thermal fluctuations of the lipid bilayers were compared with the Zilman-Granek model, 13 which is based on the Helfrich free energy model. 12

(a) Plane

(b) Vesicle

Figure 1: CG simulations of planar (a) and vesicular (b) lipid bilayers. Each lipid molecule consists of 3 connected beads – 1 head bead (black bead) and 2 tail beads (gray). The image in (b) is a cross-section of a vesicle.

9



1.03

hRo (t)i(θ,φ) / hRo i

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 10 of 31

κ κ κ κ κ

1.02

= = = = =

2.1 kB T 4.7 kB T 10.4 kB T 23.4 kB T 52.3 kB T

1.01 1.00 0.99 0.98

0

1

2

3

4

5

t [104 τ ] Figure 2: Evolution of the outer leaflet radius of the vesicle at different bending moduli, κ. The instantaneous outer leaflet radius (radially averaged in θ and φ), hRo (t)i(θ,φ) is normalized by the equilibrium outer leaflet radius, hRo (t)i.

2.4

Static and Dynamic Scattering Functions

We used a discrete Fourier transform to calculate the scattering function of the real space density distribution P(x, y, z, t) of all lipid beads, where {x, y, z} are the wrapped coordinates of a bead in a simulation box under periodic boundary conditions. {0, 0, 0} is the center of the simulation box with dimensions L = Lx = Ly = Lz , and t is the elapsed time. No distinction was made between the head and the tail beads. Each simulation snapshot — representing the positions of lipid beads at elapsed time t — was mapped on a three dimensional grid X[i, j, k] with a grid size ∆g,x = ∆g,y = ∆g,z = ∆g = 0.25 σ, number of grid points ng = L/∆g , and indices i, j, k = ⌊ x+L/2 ⌋, ⌊ y+L/2 ⌋, ⌊ z+L/2 ⌋. The complex function ρ(l, m, n) = ρ(~q, t) was ∆g ∆g ∆g calculated as, ng −1 ng −1 ng −1

ρ[l, m, n] =

X X X i=0

X[i, j, k] e(−2πi/ng )li e(−2πi/ng )mj e(−2πi/ng )nk

(5)

j=0 k=0

where the wave vector ~q has coordinates qx = 2πl/L, qy = 2πm/L and qz = 2πn/L, and q = (qx2 + qy2 + qz2 )1/2 . The static scattering function S(~q) was calculated as 1/N hρq~ ρ−~qi for correlations within the same snapshot (at the same t), whereas the dynamic scattering func10


Page 11 of 31

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


tion S(~q, t) was calculated as 1/N hρq~(∆t) ρ−~qi, assuming correlations between a snapshot captured at a reference time t and another captured after an elapsed time t+∆t. All discrete Fourier transforms were performed using the FFTW library. 32 For the S(q, t) measurements, different sampling time intervals (0.01 τ , 0.1 τ ,1 τ and 50 τ ) were taken in order to observe the decay over both short and long time scales.

3 3.1

RESULTS AND DISCUSSIONS Bending Modulus

The classical continuum description for the bending energy of a quasi-flat membrane is given by the Helfrich theory 12 as,

EBend

1 = 2

Z

dxdy κ(∆h)2 + Σ(∇h)2 ,

(6)

where κ is the bending modulus, Σ is the lateral tension, and h(x, y) is the height of the membrane relative to a reference plane. In describing h(x, y), we used the plane at z = zcm as the reference plane, and lipid molecules were sorted as either belonging to the upper or lower bilayer leaflet. To improve statistics, we independently computed h(x, y) for both the upper (hu (x, y)) and lower (hl (x, y)) leaflets such that hu (x, y) = zhead − zcm and hl (x, y) = zcm − zhead . To estimate κ we followed the procedure described by Cooke et al. 17,18 along with a two-dimensional discrete Fourier transform of the arrays, hu (x, y) and hl (x, y). By P expanding the membrane height into its Fourier modes, h(r) = h(x, y) = q hq eiq.r , and

averaging the norms of the resulting complex number arrays, we obtain the expectation

value |h2q | of the membrane height, where q = (qx2 + qy2 )1/2 . According to the linearized

Helfrich theory, the |h2q | spectrum has the following expression,

|h2q | =

kB T , Lx Ly [κq 4 + Σq 2 ] 11


(7)


where κ and Σ are fitting parameters. Although a barostat was applied in the x and y directions, the diagonal pressure tensor components pxx and pyy are not exactly zero, and typical of any MD simulation, pressure fluctuations were significant. Hence an estimate of Σ – calculated in the simulations as −Lz (pxx + pyy )/2 – was used as an initial value for the

fitting of eq.7. A typical fit of |h2q | is presented in Figure 3(a). The relationship of κ as a function of wc , obtained with different numbers of beads, N , and different thermostats in

the production runs, resulted in a best fit of κ = e 8.051 wc −9.730 for planar bilayer membranes (See Figure 3(b)). It is worth noting that the bending rigidity from simulations can be estimated by using alternative methods, which involve measuring the tensile force required for stretching a cylindrical vesicle 33 or for actively buckling a membrane. 34,35

(a)

(b) 105

102 N = 18720, Langevin N = 18720, Nos` e-Hoover N = 168480, Nos` e-Hoover

104

κ [kB T ]

q −2

103

E

Lx Ly |h2q | [σ 4 ]

102

q −4

D

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 31

101 100

101

wc = 1.6 σ ε = 1.0 kB T

10−1 −4 10

10

−3

10

10

−2

−1

10

100 1.2

0

ε = 1.0 kB T 1.3

q 2 [σ 2 ]

1.4

1.5

1.6

1.7

1.8

wc [σ]

Figure 3: Dependence of Lx Ly |h2q | on q 2 for a planar lipid bilayer system with wc = 1.6 σ and ε = 1 kB T (a). The q −2 and q −4 segments show the asymptotic tension and bending regimes, respectively. 17,18 Bending modulus κ as a function of wc for ε = 1.0kB T (b). The dashed line in (b) represents the best fit.

3.2

Small Angle Neutron Scattering of Vesicles

The vesicle’s structural parameters were determined by calculating the average distances of the lipid heads in the outer, hRo i, and inner, hRi i, bilayer leaflets from the center-of-mass of the vesicle. The brackets in Ro and Ri denote ensemble averages in both orientation and time. The vesicle bilayer thickness, hti, was calculated as the ensemble average of the 12


Page 13 of 31

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


difference between Ro (θi , φj ) and Ri (θi , φj ), where θi and φj are the polar and azimuthal angles of the radial bins, respectively. (See Figure 4.) The fluctuations or error bars in 1/2 Ro and Ri shown in Figure 4(b) is δR = hR2 i − hRi2 , similarly the bilayer thickness 1/2 fluctuation is δt = ht2 i − hti2 . The average area per lipid in the simulated vesicles,

hAl,sim i, and bilayer thickness, hti, as a function of κ are shown in Figures 4(c) and (d). The

parameter trends indicate a clear decrease in hAl,sim i and an increase in hti as κ increases, a result which is in qualitative agreement with neutron and X-ray scattering data of fluid lipid membranes. 36 The equation hAl,sim i = e(−0.104 κ−0.821) +1.06 is used to interpolate hAl,sim i for different values of κ (see dashed line in Figure 4(c)). This allows us to use the experimental κ values of the different lipid bilayer membranes to interpolate the corresponding hAl,sim i values – assuming that the experimental values corrrespond to those from simulations. The interpolated hAl,sim i provides the conversion factor from real units (˚ A) to reduced units (σ), as given in Table 1. Next, we converted published head-to-head thicknesses hti to reduced units and compared these values against the simulated results for both vesicular and planar bilayers (see Figure 4(d) and Table 1). Using this approach, we find that published bilayer thicknesses for DMPC, DOPC, POPC, DPhyPC and DLPC fall within the error bars of the simulated thicknesses, hti (See Figure 3d), validating the notion of mapping a particular lipid bilayer system to its simulation. However, it should be noted that the area per lipid, bending rigidity and membrane thickness are not necessarily coupled experimentally as they are in the present simulations. Experimentally, parameters such as the length of the acyl chains, their degree of saturation, and the relative size of the head-to-tail units can individually vary among the different lipid species, and can change the membrane thickness and its bending rigidity independently. In the current simulation approach, these membrane parameters are coupled and are primarily determined by the size of the beads and by the chosen potential, which are fixed. Hence, precise mapping of the experimental results to the simulation may require tuning ε and b in eq.2, kbend in eq.4 or changing the number of coarse-grained beads in the model lipid, items which are beyond the scope of this paper but will be addressed

13



in future publications. Tuning those parameters would allow the decoupling of membrane parameters such as the layer thickness, area per lipid and bending rigidity, permitting these variables to vary independently and consequently allowing the opportunity to simulate a wide range of lipid species as well as mixed lipid bilayers.

(a)

(b) 40.0 40

hRo i , hRi i [σ]

t

y [σ]

20 Ro

Ri

0 −20

35.0

30.0

25.0

−40 −40

−20

0

20

20.0 0 10

40

x [σ]

101

102

κ [kB T ]

(c)

(d) 1.5 5.0

hti [σ]

1.4

hAl,sim i [σ 2 ]

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 31

1.3 1.2

4.0

DMPC DHPC DOPC SOPC Plane

3.0 1.1 1.0 0 10

101

2.0 0 10

102

κ [kB T ]

101

DLPC DPPC POPC DPhyPC Vesicle

102

κ [kB T ]

Figure 4: Distances of outer, Ro , and inner, Ri , leaflets from the vesicle center-of-mass and the bilayer thickness, t (a). hRo i and hRi i as a function of bending modulus κ (b). Images in (b) are snapshots of the lipid vesicle at different values of κ. Average area per lipid, hAl,sim i of the vesicle outer leaflet (c). The curve in (c) is the best fit. Average bilayer thicknesses, hti based on the lipid head-to-head distance between leaflets (d). The scattered data points in (d) are bilayer thicknesses of common lipids mapped to the simulation results. The lines in (b) and (d) are guides to the eye.

The static scattering function, S(q), of the simulated vesicle is calculated from the density-density correlation function by taking the discrete Fourier transform of the density distribution of all lipid beads, such that S(~q) = 1/N hρq~ ρ−~qi, where ~q is the scattering 14


Page 15 of 31

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: Mapping of the different lipid membranes to simulations using experimentally determined values for the bending modulus κ, area per lipid Al , and average membrane thickness hti.

T [o C] DMPC 30 30 DLPC 48 DHPC DPPC 50 30 DOPC 30 POPC SOPC 30 30 DPhyPC

κ [kB T ] 15.8 37,38 13.2 37 9.4 39 15.0 38,39 18.2 38 20.3 41 21.5 42 12.4 43

Real Units ˚] Al [˚ A2 ] hti [A 37 60.6 35.3 37 63.2 37 30.8 37 65.1 39 38.2 39 63.1 40 38.0 40 67.4 40 36.7 40 68.3 41 37.0 41 67.0 42 39.1 42 80.5 43 36.4 43

Reduced Units Al [σ 2 ] hti [σ] 1.15 4.86 1.17 4.20 1.23 5.25 1.15 5.14 1.13 4.75 1.11 4.73 1.11 5.03 1.18 4.41

CF* [˚ A/σ] 7.27 7.34 7.28 7.39 7.73 7.83 7.77 8.25

*CF is the conversion factor used to convert real to reduced units.

vector. The static scattering function in this case, is analogous to the form factor obtained by SANS measurements, the data of which can be fitted using a core-shell model, 44

S(q) =

1 Vshell

3Vo J1 (qRo ) 3Vi J1 (qRi ) − qRo qRi

2

,

(8)

where J1 (x) = [sin(x) − x cos(x)]/x2 , Vo = 4/3 πRo3 , Vi = 4/3 πRi3 and Vshell = Vo − Vi . Fits to the scattering functions using eq.8 implemented in SasView 45 are shown in Figures 5(a), and an expanded view of the high q region is shown in Figure 5(b). The fitting using SasView takes into consideration polydispersity of the vesicle radius and shell thickness, in this case described by a Gaussian distribution. Linear correlations between the directly computed average bilayer thickness and thickness fluctuations with the fitting parameters given by SasView are shown in Figures 5(c-d). A systematic shift between the values of hti and tf it , and of δt and δt,f it can be attributed to the different binning procedures used in obtaining these values, where the real-space values were binned in a spherical coordinate system while the Fourier spectra were obtained in the Cartesian coordinate system. The Cartesian binning, in this case, is much finer than the spherical binning, such that each Cartesian bin is a cube with a dimension of 0.25 whereas the binning in the spherical coordinates is performed based

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

on θ and φ slicing, yielding 3 to 5 head-beads per bin (depending on the inner or outer leaflet). In the latter case, window averaging over each bin yields a local membrane thickness which is later used in the ensemble average over the entire vesicle. The two approaches generate instrumental errors that can be thought of in terms of the resolution of the computational methodology, which is analogous to the instrumental resolution of detectors in scattering experiments. Such resolution differences generally result in systematic shifts in measured quantities, as we observe here. The errors on membrane thickness values associated with the binning procedures are estimated to be 0.25 σ in Cartesian coordinates and ∼ 0.3 σ to 0.5 σ (depending on bending rigidity) in spherical coordinates. In addition, the more pronounced fringes in the spectra, due to the lower radial polydispersity in the more rigid systems (κ = 23.4 kB T and κ = 52.3 kB T ), render a higher uncertainty in the fits of the minimum whose position defines the membrane thickness. All these factors ought to be taken into account when comparing the results of the two approaches. Similar trends are observed in SANS experiments on lipid vesicles, where vesicle radius is determined from the low-q fringe patterns and bilayer thickness is determined from the high-q data. If the background level is low enough, as is the case with acyl chain deuterated vesicles in a deuterated buffer, a minimum in the form factor can also be observed in the data at high q. 23 In experiments, this intensity minimum shifts to lower q values as the lipids transition from a fluid to a gel state; i.e., as the bending rigidity increases, indicating an increase in bilayer thickness (e.g., see Fig. 1 in ref. 23 and Fig. S4 in ref. 24 ). Likewise, the simulations show that as the bending modulus of the lipids comprising the vesicle increases, the minimum in S(q), i.e., q0 , moves toward lower q values, in accord with the thickening of the bilayer. That is to say that q0 decreases with increasing hti and κ. We note that the present calculations of the form factors do not consider the contrast between the membrane and the surrounding medium, as is the case in experiments, and are solely based on the structural correlations between the lipids. Modifications of such calculations are underway for direct comparison with SANS data in order to capture and interpret variations that

16


Page 16 of 31

Page 17 of 31

result from compositional asymmetry or lipid imbalance between the inner and outer bilayer leaflets.

(a)

(b) = = = = =


S(q)

S(q)

κ κ κ κ κ

10−2

100

10−1

101

0.8

1.0

1.2

−1

q [σ ]

1.4 1.6 −1

1.8

2.0

q [σ ]

(c)

(d) 4.8 0.6

2.1

52.3

4.6

0.5

4.4

δt [σ]

23.4

hti [σ]

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.4

0.4 4.7 0.3

4.7 0.2

4.2 2.1

23.4

10.4

52.3

0.1 4.0 4.8

5.0

5.2

5.4

5.6

0.0 0.0

5.8

tf it [σ]

0.1

0.2

0.3

0.4

0.5

0.6

0.7

δt,f it [σ]

Figure 5: Static scattering function S(q) of vesicles with lipids having different bending moduli, κ. Solid curves are fits to the data using eq.8 obtained from fits using SasView. The spectra are shifted in the vertical axis for clarity. The spectra in the high q region are shown in (b). The vertical solid arrows in (b) are the minima of the spectra, q0 . Comparison between the average bilayer thickness hti and the fitting parameter for the thickness obtained from SasView, tf it (c). Comparison between the average thickness fluctuation δt and the standard deviation of the thickness distribution obtained from SasView, δt,f it (d). The dashed line in (c) is the best fit. In (d) the values in the vicinity of the dashed line, f (x) = x, are of the bending modulus κ.

3.3

Neutron Spin Echo

NSE spectroscopy is a powerful technique for investigating lipid bilayer bending and thickness fluctuations as it directly accesses the length- and time-scales over which those fluctuations 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

Page 18 of 31

take place. 20,23,46 An NSE experiment’s output is the intermediate dynamic scattering function, which can be calculated in simulations as S(q, t)/S(q, 0), where the scattering function S(q) at time t is correlated with S(q) after an elapsed time, ∆t as 1/N hρq~(∆t) ρ−~qi. The function is normalized with S(q) such that the S(q, t)/S(q, 0) at ∆t = 0 is unity. S(q, t)/S(q, 0) provides the temporal decay of collective dynamic modes of the membrane. Figure 6(a) shows S(q, t)/S(q, 0) for q values betweenn 0.85 σ −1 and 1.70 σ −1 , a range relevant to membrane fluctuations and complementary to the S(q) range shown in Figure 5(b). The lines in Figure 6(a) are fits to the equation,

β

S(q, t)/S(q, 0) ≈ (1 − A) e−ΓT rans t + A e−(Γ t) ,

(9)

where ΓT rans is the relaxation time of the collective motion or relative translation of individual lipids as a result of the center-of-mass motion of the vesicle; this mode is represented by a single exponential decay. The bending modes follow a stretched exponential with β = 2/3, as predicted by Zilman and Granek. 13 Ideally S(q, t)/S(q, 0) is described by the coupled motion of translation and bending undulation modes, such that β

S(q, t)/S(q, 0) = e−ΓT rans t [(1 − A) + A e−(Γ t) )]. 13,47 However, this equation is not easily applied. Instead the linear combination of modes described by eq.9 is used under the assumption that at the simulation time window the translational motion is non-relaxed, thus exp(−ΓT rans t) ≈ 1 − o(ΓT rans t). 20,46 The fits of the data to eq.9 show a variety of dynamic modes at different q values, as seen in Figure 6(b). We should first clarify that A = 1 implies that the dynamics are purely bending fluctuations and that deviations of A from unity signify the onset of other dynamic modes in the system. Figure 6(b) shows that at small q values, the coefficient of bending motions, A ≈ 1, indicating that bending modes are dominant at these length scales and the bilayer behaves like a continuum surface with negligible translational motion. For large q’s, on the other hand, the coefficient A plateaus at values less than 1, with the biggest deviation observed for the softest vesicle (κ = 2.1 kB T ). These dynamics are attributed to translational 18


Page 19 of 31

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


motions, which are expected to be more noticeable at these length scales. Keep in mind that the translational motions result from the restoring spring forces that are used to peg the vesicle to the center of the simulation box, and that the spring corrects the trajectory of individual beads over very short distances (i.e., large q) in a fashion reminiscent of thermal fluctuations or protrusions of individual lipid molecules. In this case, softer membranes are expected to experience more significant translational dynamics, which is indeed what we observe in Figure 6(b). The most peculiar observation, however, is the dip in the A values at the intermediate q range, indicating excess dynamics beyond bending fluctuations and translational motions. Examination of the dip positions shows that they occurs at q values that correspond to the membrane thickness, q0 , in excellent agreement with NSE data from lipid membranes (see Fig. S3 in ref. 24 ). Since these dynamics occur exactly at length scales of the membrane thickness, they are attributed to membrane thickness fluctuations. At these length scales, the bilayer is seen to behave like two separate surfaces, where the out-of-phase surface leaflet fluctuations (thickness fluctuation modes) are significantly pronounced. If we subtract the translational motions, we can then safely assume that for q values far from q0 bending modes are dominant and are adequately described by ΓZG ∼ q 3 as predicted by Zilman and Granek 13 (see Figure 6(d)), while thickness and bending fluctuations modes are coupled at values near q0 . The patterns in Figure 6(c) are extracted using the approach by Nagao et al., 22,23,48 in which the bending and thickness fluctuation relaxation modes are empirically expressed as the sum of two terms, Γ 1 ΓZG ΓT F = 3 + 3 , 3 q q q 1 + (q − q0 )2 ξ −2

(10)

where ΓZG indicates the decay rate due to the bending fluctuations, as predicted by Zilman and Granek, 13 and ΓT F pertains to the decay rate attributed to thickness fluctuations that manifest themselves as excess dynamics observed at q = q0 . The second term in the r.h.s of eq.10 is a Lorentzian function, where q0 , ΓT F /q03 and ξ −1 are the peak position, peak height and half width at half maximum (HWHM) of the Lorentzian, respectively. We emphasize 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

that the decay rates in Figure 6(c) closely follow the trends of experimentally observed decay rates from NSE studies of lipid membranes, with shifts of the peak positions towards lower q values for more rigid membranes and a concomitant narrowing of the peak width (see Fig. 2 in ref. 23 and Fig. S2 in ref. 24 ). The width of the peaks in the decay rates in q space can be thought of as an uncertainty on q0 , such that δq ≈ ξ −1 . Using error analysis and knowing that q0 ≈ 1/ hti, one can easily see that δt ≈ ξ −1 . In this framework, the narrowing of the peak is consistent with the notion of smaller thickness fluctuation amplitudes, and the results then become physically intuitive. As the membrane stiffens when it transitions into a gel phase, it becomes thicker and its thickness fluctuations become less prominent, as expected. We checked the validity of the combined approach of Nagao et al. 22,23,48 and Monroy et al. 20,46 to describe the collective motions of the simulated bilayer membrane, by comparing the computed fluctuations, δt , with the fluctuations indirectly obtained from a convolution of a Gaussian (center-of-mass translation) with a Lorentzian (bending motion, and thickness fluctuations obtained from eq.10 ) distributions. (See Figure 7.) If this approach is valid, we expect δt to be comparable with the full width at half maximum (FWHM) of the convolution of the Gaussian with the Lorentzian distribution. Recall that we applied a harmonic spring, Ucm (rcm ) = 1/2 ks (~rcm )2 , which pegs the center-of-mass position of the vesicle ~rcm at {0,0,0} in the Cartesian coordinate system of the simulation box. The fluctuations of ~rcm about the origin can be computed as σcm = 1/2 2 hrcm i − hrcm i2 , and σcm can be analytically estimated from Ucm by using the Boltzmann R∞ 4 R∞ 2 2 factor, where σcm = 0 rcm exp(−Ucm /kB T )dr/ 0 rcm exp(−Ucm /kB T )dr = 3kB T /ks . (See dashed line Figure 7). The estimated σcm clearly describes the vesicles made up of lipids

forming membranes with a larger κ. The FWHM of the center-of-mass fluctuations, assuming p a Gaussian distribution, is αG = 8 ln(2)σcm . Next, the thickness fluctuation amplitude of the membrane obtained using eq.10 is estimated as ∆t = hti ξ −1 /qo , hence the FWHM of

the Lorentzian function is αL = 2∆t . An estimate of the FWHM of the convolution of a

20


Page 20 of 31

Page 21 of 31

(a)

(b) κ = 2.1 kB T 1.0 0.8

q = 0.85 σ −1

0.8

q = 0.90 σ −1 q = 1.00 σ −1

0.6

q = 1.10 σ −1

A

S(q, t)/S(q, 0)

1.0

q = 1.25 σ −1

0.4

q = 1.30 σ

0.6

−1

0.4

q = 1.35 σ −1

0.2

q = 1.45 σ −1

0.2

q = 1.55 σ −1 q = 1.70 σ −1

0.0 −2 10

10−1

100

101

0.0 0.8

102

t [τ ]

κ κ κ κ κ

= = = = =


1.0

1.2

1.4 −1

1.6

1.8

q [σ ]

(c)

(d) 10−2

10

−1 κ κ κ κ κ

10−2

= = = = =


ΓZG ∼ κ−1/2

ΓZG [τ −1 ]

Γ/q 3 [σ 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


10−3

10−4 0.8

1.0

1.2

1.4 −1

1.6

10−3

10−4 0 10

1.8

q [σ ]

101

102

κ [kB T ]

Figure 6: S(q, t)/S(q, 0) at different values of q near q0 for lipids that form membranes with κ = 2.1 kB T (a). The bending prefactor, A from the S(q, t)/S(q, 0) fits using eq.9 at different values of κ (b). Decay rate, Γ from fits using eq.9 at different values of κ (c). The lines in (c) are Lorentzian fits using eq.10. Dependence of ΓZG as a function of κ are obtained from fits using eq.10 (d).

21



Gaussian with a Lorentzian function is given as, 49

αv =

1 2 1/2 , C1 αL + (C2 αL2 + 4C3 αG ) 2

(11)

where C1 = 1.0692, C2 = 0.86639 and C3 = 1.0. In terms of the quantities ∆t and σcm , eq.11 2 1/2 becomes αv = 1.0692∆t + (0.86639∆2t + 8 ln(2)σcm ) . Figure 7 shows excellent agreement

between δt and αv , validating our combined approach.

Fluctuations [σ]

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 22 of 31

δt σcm ∆t αv

0.6

0.4

0.2

0.0

100

101

102

κ [kB T ] Figure 7: Direct calculation of the average thickness fluctuation amplitude, δt , and average vesicle center-of-mass fluctuations, σcm , along with the estimated average thickness fluctuation amplitudes, ∆t , (based on eq.10) and the average effective fluctuation amplitude, αv (based on eq.11). The dashed line is the estimated average center-of-mass of the fluctuations (3kB T /ks )1/2 .

Note that the system with the smallest bending rigidity (κ = 2.1 kB T ) exhibits the largest center-of-mass fluctuations. This is because the spring applied to all lipids fixes the unrwapped center-of-mass coordinates to the origin. Since this system coexists with the unstable phase – where there are many stray lipids that move to other periodic images – the unwrapped center-of-mass does not correspond to the center-of-mass of the original simulation box. To remedy this and to prevent the stray lipids from transferring to other periodic box images, we applied a repulsive WCA wall on the faces of the simulation box. With this adjustment, the center-of-mass fluctuations of the system became comparable to 22


Page 23 of 31

the analytical estimate, i.e., (3kB T /ks )1/2 . Also, calculating the intermediate scattering function with the system having a repulsive wall broadens the Lorentzian fit in the gamma Γ/q 3 vs. q plot. (See Figure 8(a)) This is attributed to a decrease in the center-of-mass fluctuations, which enhance the thickness fluctuations calculated from S(q, t) . However, the broadening of the Lorentzian compensates the decrease in the center-of-mass fluctuations keeping αv the same for both systems. The total thickness flucations, δt , is also approximately equal, which is expected as the application of a global spring force should not affect bilayer properties, such as δt . (See Figure 8(b))

(a)

(b) κ = 2.1 kB T 10−1

0.6 w/o wall

Fluctuations [σ]

w/ wall

Γ/q 3 [σ 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


10−2

0.4

0.2

0 10−3 0.8

1.0

1.2

1.4 −1

1.6

1.8

δt

q [σ ]

δcm w/o wall

∆t

αv

w/ wall

Figure 8: Decay rate, Γ from fits using eq.9 for lipids with κ = 2.1 kB T which are free or enclosed in a WCA wall. (a) The wall is described by eq.1 with b = σ. Comparison of δt , δcm, ∆t and αv between systems with lipids that are free of or enclosed by a WCA wall (b).

4

CONCLUSIONS

We performed CG MD simulations on experimentally relevant planar and vesicular lipid bilayers with different bending rigidities, similar of those used in scattering measurements. Importantly, we demonstrated that the adopted computational approach can efficiently reproduce neutron scattering data by computing the static, S(q), and normalized dynamic, S(q, t)/S(q, 0), structure factors, analogous to SANS and NSE data, respectively. The struc23



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

tural S(q) of the simulated vesicles was well-described by a core-shell model with a Gaussian distribution of shell thicknesses and vesicle radii. The validation of this approach using simple, single-component vesicles is encouraging for its future use as a tool for the direct analysis of SANS data on more complex systems, such as asymmetric vesicles or multicomponent lipid membranes, in which subtle deviations in the scattering signal can be indicative of important structural rearrangements within the bilayer. The simulations were also utilized to explore collective membrane dynamics that are challenging for all-atom simulations, and are of the type generally accessed by NSE, namely bending undulations and thickness fluctuations. Bending moduli were estimated using the Helfrich Hamiltonian 12 and the procedure described by Cooke et al. 17,18 for planar bilayer systems. The simulations show excellent agreement with the scaling predictions of Zilman and Granek, 13 ΓZG ∼ κ−1/2 q 3 , as observed in NSE and dynamic light scattering experiments. The simulations also reproduced the excess dynamics, beyond bending fluctuations, uniquely observed in NSE measurements at length scales commensurate with the bilayer thickness. These excess dynamics are attributed to membrane thickness fluctuations and closely follow experimentally observed trends. This required a decoupling of the collective motions captured in the simulations and described by S(q, t), namely translational motions, bending undulations and thickness fluctuations. To this end, we used an approach incorporating the methods of Nagao et. al. 22,23,48 — which separates thickness fluctuations from the bending motion — and that of Monroy et al. 20,46 — which couples the bending motion and thickness fluctuations with translational motion. We traced the translational motion as originating from the fluctuations of the center-of-mass of the vesicles that manifests itself as translational motion in the individual lipids. This enabled us to describe the total fluctuations of the vesicle thickness as the convolution of fluctuations from the center-of-mass of the vesicle, and the combined bending and thickness fluctuations of the membrane. More notably, the confirmation of the presence of an additional dynamical mode corresponding to membrane thickness fluctuations and the remarkable resemblance of simulation results to experimental

24


Page 24 of 31

Page 25 of 31

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


data trends of bending and thickness fluctuations is most promising about our computational approach, as it opens new possibilities for combining simulations and scattering experiments to explore long-standing problems in membrane physics, such as the mechanical coupling of membrane leaflets and the effect of phase separation on thickness fluctuation modes in mixed lipid bilayer systems.

Acknowledgement This work was conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. J. K. is supported by DOE’s Scientific User Facilities Division. R. A. is supported by ORNL’s Neutron Sciences Directorate Clifford G. Shull Fellowship. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

References (1) Phillips, R.; Ursell, T.; Wiggins, P.; Sens, P. Emerging Roles for Lipids in Shaping Membrane-Protein Function. Nature 2009, 459, 379–385. (2) König, S.; Sackmann, E. Molecular and Collective Dynamics of Lipid Bilayers. Curr. Opin. Colloid & Interface Sci. 1996, 1, 78–82. (3) Evans, E.; Rawicz, W. Entropy-Driven Tension and Bending Elasticity in CondensedFluid Membranes. Phys. Rev. Lett. 1990, 64, 2094. (4) Vogel, A.; Roark, M.; Feller, S. E. A Reinterpretation of Neutron Scattering Experiments on a Lipidated Ras Peptide Using Replica Exchange Molecular Dynamics. BBABiomembranes 2012, 1818, 219–224.

25



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

(5) Engelbrecht, T.; Hauß, T.; S¨ uβ, K.; Vogel, A.; Roark, M.; Feller, S. E.; Neubert, R. H.; Dobner, B. Characterisation of a New Ceramide Eos Species: Synthesis and Investigation of the Thermotropic Phase Behaviour and Influence on the Bilayer Architecture of Stratum Corneum Lipid Model Membranes. Soft Matter 2011, 7, 8998–9011. (6) Wohlert, J.; Edholm, O. Dynamics in Atomistic Simulations of Phospholipid Membranes: Nuclear Magnetic Resonance Relaxation Rates and Lateral Diffusion. J. Chem. Phys. 2006, 125, 204703. (7) Kim, T.; Lee, K. I.; Morris, P.; Pastor, R. W.; Andersen, O. S.; Im, W. Influence of Hydrophobic Mismatch on Structures and Dynamics of Gramicidin A and Lipid Bilayers. Biophys. J. 2012, 102, 1551–1560. (8) Nickels, J. D.; Cheng, X.; Mostofian, B.; Stanley, C.; Lindner, B.; Heberle, F. A.; Perticaroli, S.; Feygenson, M.; Egami, T.; Standaert, R. F.; Smith, J. C.; Myles, D. A. A.; Ohl, M.; Katsaras, J. Mechanical Properties of Nanoscopic Lipid Domains. J. Am. Chem. Soc. 2015, 137, 15772–15780. (9) Reynwar, B. J.; Illya, G.; Harmandaris, V. A.; M¨ uller, M. M.; Kremer, K.; Deserno, M. Aggregation and Vesiculation of Membrane Proteins by Curvature-Mediated Interactions. Nature 2007, 447, 461–464. (10) Huang, K. C.; Mukhopadhyay, R.; Wingreen, N. S. A Curvature-Mediated Mechanism for Localization of Lipids to Bacterial Poles. PLoS Comput. Biol. 2006, 2, e151. (11) Lindahl, E.; Edholm, O. Mesoscopic Undulations and Thickness Fluctuations in Lipid Bilayers from Molecular Dynamics Simulations. Biophys. J. 2000, 79, 426–433. (12) Helfrich, W. Elastic Properties of Lipid Bilayers: Theory and Possible Experiments. Zeitschrift f¨ ur Naturforschung C 1973, 28, 693–703.

26


Page 26 of 31

Page 27 of 31

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


(13) Zilman, A.; Granek, R. Undulations and Dynamic Structure Factor of Membranes. Phys. Rev. Lett. 1996, 77, 4788. (14) 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. (15) Watson, M. C.; Penev, E. S.; Welch, P. M.; Brown, F. L. Thermal Fluctuations in Shape, Thickness, and Molecular Orientation in Lipid Bilayers. J. Chem. Phys. 2011, 135, 244701. (16) 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. (17) 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. (18) Cooke, I. R.; Kremer, K.; Deserno, M. Tunable Generic Model for Fluid Bilayer Membranes. Phys. Rev. E 2005, 72, 011506. (19) Dimova, R. Recent Developments in the Field of Bending Rigidity Measurements on Membranes. Adv. Colloid Interface Sci. 2014, 208, 225–234. (20) Arriaga, L. R.; López-Montero, I.; Monroy, F.; Orts-Gil, G.; Farago, B.; Hellweg, T. Stiffening Effect of Cholesterol on Disordered Lipid Phases: a Combined Neutron Spin Echo+ Dynamic Light Scattering Analysis of the Bending Elasticity of Large Unilamellar Vesicles. Biophys. J. 2009, 96, 3629–3637. (21) Farago, B. Recent Developments and Applications of NSE in Soft Matter. Curr. Opin. Colloid & Interface Sci. 2009, 14, 391–395.

27



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

(22) Nagao, M. Observation of Local Thickness Fluctuations in Surfactant Membranes Using Neutron Spin Echo. Phys. Rev. E 2009, 80, 031606. (23) Woodka, A. C.; Butler, P. D.; Porcar, L.; Farago, B.; Nagao, M. Lipid Bilayers and Membrane Dynamics: Insight into Thickness Fluctuations. Phys. Rev. Lett. 2012, 109, 058102. (24) Ashkar, R.; Nagao, M.; Butler, P. D.; Woodka, A. C.; Sen, M. K.; Koga, T. Tuning Membrane Thickness Fluctuations in Model Lipid Bilayers. Biophys. J. 2015, 109, 106–112. (25) Ingólfsson, H. I.; Arnarez, C.; Periole, X.; Marrink, S. J. Computational ‘Microscopy’of Cellular Membranes. J. Cell Sci. 2016, 129, 257–268. (26) Lee, E. H.; Hsin, J.; Sotomayor, M.; Comellas, G.; Schulten, K. Discovery Through the Computational Microscope. Structure 2009, 17, 1295–1306. (27) Dobrynin, A. V.; Carrillo, J.-M. Y. Universality in Nonlinear Elasticity of Biological and Polymeric Networks and Gels. Macromolecules 2010, 44, 140–146. (28) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1–19. (29) Brown, W. M.; Wang, P.; Plimpton, S. J.; Tharrington, A. N. Implementing Molecular Dynamics on Hybrid High Performance Computers–Short Range Forces. Comp. Phys. Comm. 2011, 182, 898–911. (30) Shinoda, W.; Shiga, M.; Mikami, M. Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation under Constant Stress. Phys. Rev. B 2004, 69, 134103. (31) Schneider, T.; Stoll, E. Molecular-Dynamics Study of a Three-Dimensional OneComponent Model for Distortive Phase Transitions. Phys. Rev. B 1978, 17, 1302.

28


Page 28 of 31

Page 29 of 31

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


(32) Frigo, M.; Johnson, S. G. Acoustics, Speech and Signal Processing, 1998. Proceedings of the 1998 IEEE International Conference on; 1998; Vol. 3; pp 1381–1384. (33) 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. (34) Hu, M.; Diggins Patrick, P.; Deserno, M. Determining the Bending Modulus of a Lipid Membrane by Simulating Buckling. J. Chem. Phys. 2013, 138, 214110. (35) Noguchi, H. Anisotropic Surface Tension of Buckled Fluid Membranes. Phys. Rev. E 2011, 83, 061919. (36) Kuˇcerka, N.; Nieh, M.-P.; Katsaras, J. Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochim. Biophys. Acta 2011, 1808, 2761–2771. (37) Kuˇcerka, N.; Liu, Y.; Chu, N.; Petrache, H. I.; Tristram-Nagle, S.; Nagle, J. F. Structure of Fully Hydrated Fluid Phase DMPC and DLPC Lipid Bilayers Using X-ray Scattering from Oriented Multilamellar Arrays and from Unilamellar Vesicles. Biophys. J. 2005, 88, 2626–2637. (38) Nagle, J. F. Introductory Lecture: Basic Quantities in Model Biomembranes. Farad. Discuss. 2013, 161, 11–29. (39) Guler, S. D.; Ghosh, D. D.; Pan, J.; Mathai, J. C.; Zeidel, M. L.; Nagle, J. F.; TristramNagle, S. Effects of Ether versus Ester Linkage on Lipid Bilayer Structure and Water Permeability. Chem. Phys. Lipids 2009, 160, 33–44. (40) Kuˇcerka, N.; Nagle, J. F.; Sachs, J. N.; Feller, S. E.; Pencer, J.; Jackson, A.; Katsaras, J. Lipid Bilayer Structure Determined by the Simultaneous Analysis of Neutron and XRay Scattering Data. Biophys. J. 2008, 95, 2356–2367.

29



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

(41) Kuˇcerka, N.; Tristram-Nagle, S.; Nagle, J. F. Structure of Fully Hydrated Fluid Phase Lipid Bilayers with Monounsaturated Chains. J. Membr. Biol. 2006, 208, 193–202. (42) Greenwood, A. I.; Pan, J.; Mills, T. T.; Nagle, J. F.; Epand, R. M.; Tristram-Nagle, S. CRAC Motif Peptide of the HIV-1 GP41 Protein Thins SOPC Membranes and Interacts with Cholesterol. Biochim. Biophys. Acta 2008, 1778, 1120–1130. (43) Tristram-Nagle, S.; Kim, D. J.; Akhunzada, N.; Kuˇcerka, N.; Mathai, J. C.; Katsaras, J.; Zeidel, M.; Nagle, J. F. Structure and Water Permeability of Fully Hydrated DiphytanoylPC. Chem. Phys. Lipids 2010, 163, 630–637. (44) Guinier, A.; Fournet, G. Small Angle Scattering of X-Rays; John Wiley & Sons: New York, 1955. (45) SasView, 3rd ed.; http://www.sasview.org. (46) Arriaga, L. R.; López-Montero, I.; Orts-Gil, G.; Farago, B.; Hellweg, T.; Monroy, F. Fluctuation Dynamics of Spherical Vesicles: Frustration of Regular Bulk Dissipation into Subdiffusive Relaxation. Phys. Rev. E 2009, 80, 031908. (47) Milner, S. T.; Safran, S. Dynamical Fluctuations of Droplet Microemulsions and Vesicles. Phys. Rev. A 1987, 36, 4371. (48) Nagao, M.; Chawang, S.; Hawa, T. Interlayer Distance Dependence of Thickness Fluctuations in a Swollen Lamellar Phase. Soft Matter 2011, 7, 6598–6605. (49) Olivero, J.; Longbothum, R. Empirical Fits to the Voigt Line Width: A Brief Review. J. Quant. Spectrosc. Radiat. Transfer 1977, 17, 233–236.

30


Page 30 of 31

Page 31 of 31


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 ACS Paragon Plus Environment

A Computational Approach for Modeling Neutron ... - ACS Publications

Recommend Documents