Molecular Dynamics Simulations and Neutron ... - ACS Publications

Subscriber access provided by UNIV MASSACHUSETTS BOSTON

Article

Molecular dynamics simulations and neutron reflectivity as an effective approach to characterize biological membranes and related macromolecular assemblies Leonardo Darré, Javier Iglesias-Fernández, Axel Kolhmeyer, Hanna Wacklin, and Carmen Domene J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b00635 • Publication Date (Web): 17 Sep 2015 Downloaded from http://pubs.acs.org on September 30, 2015

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 35

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

1

Molecular dynamics simulations and neutron

2

reflectivity as an effective approach to characterize

3

biological membranes and related macromolecular

4

assemblies

5

L. Darré,1 J. Iglesias-Fernandez,1 Axel Kohlmeyer,2 H. Wacklin,3,4 and C. Domene1,5

6 7

1

8 9

2

Department of Chemistry, King’s College London, Britannia House, 7 Trinity Street, London SE1 1DB, UK Institute for Computational Molecular Science (035-07), College of Science and Technology, Temple University, 1901 N. 13th Street, Philadelphia, PA 19122, USA 3

10 4

11 12 13

5

European Spallation Source ESS AB, P.O.Box. 176, 221 00 Lund, Sweden

Physical Chemistry, Department of Chemistry, Lund University, P.O. Box 124, 22100 Lund, Sweden

Chemistry Research Laboratory, Mansfield Road, University of Oxford, Oxford OX1 3TA, UK

14

1

15

KEYWORDS: Atomistic; Coarse-Grained; NSLD; Voronoi; Abeles Matrix; NeutronRefTools

Corresponding authors: [email protected] & [email protected]

16

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 35

1

ABSTRACT

2

In combination with other spectroscopy, microscopy, and scattering techniques, neutron

3

reflectivity is a powerful tool to characterize biological systems. Specular reflection of neutrons

4

provides structural information at the nano- and subnano-meter length scales, probing the

5

composition and organization of layered materials. Currently, analysis of neutron reflectivity

6

data involves several simplifying assumptions about the structure of the sample under study,

7

affecting the extraction and interpretation of information from the experimental data. Computer

8

simulations can be used as a source of structural and dynamic data with atomic resolution. We

9

present a novel tool to compare the structural properties determined by neutron reflectivity

10

experiments with those obtained from molecular simulations. This tool allows benchmarking the

11

ability of molecular dynamics simulations to reproduce experimental data, but also promotes

12

unbiased interpretation of experimentally determined quantities. Two application examples are

13

presented to illustrate the capabilities of the new tool. The first example is the generation of

14

reflectivity profiles for a DMPC lipid bilayer from molecular dynamics simulations using data

15

from both atomistic and coarse-grained models, and comparison with experimentally measured

16

data. The second example is the calculation of lipid volume changes with temperature and

17

composition from all atoms simulations of single and mixed DOPC and DPPC bilayers.

18


2

Page 3 of 35

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


TOC

2 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

1 2

Page 4 of 35

Introduction Surfaces, interfaces and thin films, ubiquitous in nature, are planar structures defining the

3

spatial boundary between two different media. In biological systems, boundaries between many

4

phases are defined by soft interfaces, (i.e. membranes and biopolymers) embedded in

5

physiological electrolyte solutions. For example, biological membranes are vital components

6

delimiting the outer boundary of the cell from the environment as well as that of intracellular

7

compartments within the cytoplasm. They serve as smart filters allowing for functional

8

compartmentalization, protection from toxic substances and transport of metabolites.

9

Additionally, they provide a two-dimensional physico-chemical framework enhancing

10

biochemical processes, such as diffusion-limited reactions.1 Consequently, understanding the

11

structure of this type of systems constitutes a cornerstone of multiple biological and

12

technological applications. These complex systems are studied by using both integrated and

13

reconstituted model systems that closely mimic natural membranes.

14

Neutron reflectometry is a powerful tool for the characterization of thin films and interfaces in

15

the nano- and subnano-meter length scales2, 3 and is particularly advantageous due to its ability to

16

detect the dimensions, chemical composition and relative location of specific molecular species,

17

in non-crystalline, disordered thin films. This ability relies on several factors. The first of them is

18

the sensitivity of neutrons to the coherent nuclear scattering length, an intrinsic property of the

19

elements, which displays a particularly large difference between hydrogen (-3.74 fm) and its

20

heavy isotope deuterium (+6.67 fm) allowing isotopic substitution to be exploited to highlight

21

specific parts of the system. Other beneficial factors are their nondestructive nature and their

22

long penetration ability. In addition, light elements such as hydrogen, carbon, nitrogen, and


4

Page 5 of 35

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

oxygen strongly scatter neutrons, thus making neutron reflectivity suitable for biological

2

samples.

3

In a specular neutron reflection experiment, reflectivity, i.e. the ratio between the reflected and

4

incident intensities (Equation 1) of a collimated neutron beam directed at an interface at a

5

grazing angle (Figure 1B) is measured as a function of the momentum transfer vector Qz

6

(Equation 2), which depends on the incident neutron beam angle (θ) and wavelength (λ). I I0

7

R=

Equation 1

8

Qz =

9

The molecular organization in the direction orthogonal to the planar interface, i.e. the z axis,

10

determines the density of atoms at every value of z and, as the neutron scattering length depends

11

on the atomic nuclei, this defines the Neutron Scattering Length Density (NSLD) profile in that

12

direction. The NSLD profile determines the reflectivity, thus the latter contains the information

13

about the molecular organization of the system. However, such information is not directly

14

available from the reflectivity data, and typically, an initial structural guess to model the system

15

is required. From this initial guess, the reflectivity can be calculated, using, for example, the

16

Abeles optical matrix method4 (see Methods), and compared to the experimentally measured

17

data. This procedure is iterated to reach an agreement between the modeled and the experimental

18

reflectivity. The final outcome is a NSLD profile associated with the model generated for the

19

molecular organization of the system. One drawback of the method is that, like in all other

20

diffraction-based methods, phase information is lost on measuring the reflected intensity, making

21

it difficult to determine a unique model. To overcome this limitation, multiple neutron contrasts,

22

i.e. solute/solvent deuteration schemes, are experimentally measured and simultaneously fitted

4 ⋅ π ⋅ sin ( θ )

Equation 2

λ


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 35

1

against different models. This allows incorrect models to be filtered out, often aided by

2

additional constraints like molecular connectivity or physical properties known about the system

3

(e.g. physical density).

4

Model NSLD profiles can be also obtained from computer simulations. In molecular dynamics

5

(MD) simulations, the ensemble average structure of a system at equilibrium is readily available,

6

and thus, a direct model of the NSLD (NSLDMD) is calculable. This information can be used to

7

generate a reflectivity profile from the simulated trajectory, which can be directly compared to

8

the experimental data or used to build an MD-model-based initial guess for the fitting procedure

9

mentioned above. This is particularly appealing considering the ever-growing performance of

10

computers, maturity and efficiency of MD simulation software, and the availability of coarse-

11

grained (CG) models that allow overcoming length and time scales limitations of atomistic MD

12

simulations. Similar approaches have already been applied, for example, to fit neutron

13

reflectivity (NR) data of the α-hemolysin ion channel reconstituted in tethered lipid bilayers by

14

generating a model of the contribution of the protein to the NSLD profile from a crystal

15

structure.5 Another example is the software SIMtoEXP

16

(http://sourceforge.net/projects/simtoexp/),6 which combines MD data with the analysis of small

17

angle neutron and X-ray scattering (SANS/SAXS) data.

18

NeutronRefTools is a plugin for the VMD molecular visualization and analysis software,7 and

19

allows extraction of the NSLDMD from either all-atom (AA) or coarse-grained (CG) simulations,

20

as well as the calculation of the neutron reflectivity profile. This can be extended to the

21

corresponding electron density and X-ray reflectivity profiles. The corresponding experimental

22

data can also be loaded and directly compared to the simulations. Additionally, NeutronRefTools

23

is able to calculate volume fractions and absolute volumes of the components of the system,


6

Page 7 of 35

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

providing further structural detail, and basis for molecular scattering length densities models.

2

The combination of these features and the integration into the VMD software package facilitates

3

direct comparison of data, without the need to transfer or convert data sets between MD and

4

neutron scattering software. Two application examples are presented to illustrate the capabilities

5

of the tool. In test case 1, reflectivity profiles obtained from AA and CG MD simulations are

6

compared with experimental data of a DMPC lipid bilayer measured in different solvent

7

contrasts. Additionally, the molecular volumes, both volume fraction profiles and absolute

8

values, are computed and compared to previous literature reports. In test case 2, the molecular

9

volumes in pure and binary DOPC/DPPC lipid bilayers are obtained from AA simulations. This

10

example illustrates the ability of the approach to obtain molecular insight into more

11

experimentally challenging systems.

12

Methods

13

The NeutronRefTools plugin computes four properties of a system: i) the NSLD profile, ii) the

14

volume fraction (Φ) profile of the components, iii) absolute molecular volumes, and iv) the

15

reflectivity profile of the system at an interface that can be defined by the user.

16

Both, the NSLD and volume fraction profiles are one-dimensional properties computed along

17

the direction orthogonal to the surfaces of the layers of the system. As an example, a DMPC lipid

18

bilayer/water system with the membrane aligned to the x-y plane will be considered. The

19

dimension in which the NSLD and volume fractions are measured corresponds to the z axis (see

20

Figure 1A).

21 22 23

Once the direction is defined, the system is discretized in layers (from i=3 to N) along that direction, and the NSLD is calculated in each volume element according to Equation 3: NSLD MOL = j

nk ⋅ bk Vlayer k = atomtypes MOL ∈ layeri i j

∑


Equation 3

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

1

Page 8 of 35

where nk is the number of atoms of type k corresponding to a given molecular species MOLj, in

2

a given layer layeri, bk is the coherent neutron scattering length for the atom type, and Vlayer,i is

3

the volume of the layer defined by the x and y dimensions of the box, and the selected bin size in

4

the z-direction. The total NSLD for a given layer is computed by addition of the NSLD

5

contribution of all molecular species present in the layer.

6

7

The volume fraction (Φ) is calculated using Equation 4:

Φ MOL = j

k=

∑Vk atoms MOL j ∈ layeri

Equation 4

Vlayer

i

8

where Vk is the volume of an atom from a given molecular species MOLj located in a given

9

layer Layeri. The calculation of each Vk is performed by using a 3D Voronoi tessellation of the

10

system (NeutronRefTools uses the Voro++ library8 for this purpose) per layer or for the

11

complete system box, where the former is slower than the latter.

12

These calculations can be iterated over selected MD trajectory frames providing average

13

values per layer, which are finally used to generate the NSLDMD(z) and ΦMD(z) profiles. The

14

NSLDMD(z) is the quantity compared to experimentally fitted NSLD profiles, and its

15

decomposition in molecular contributions, as well as that of the ΦMD(z), allows interpreting the

16

structure of the system in terms of the spatial distribution of its components. NeutronRefTools

17

allows the user to decompose the NSLDMD(z) and ΦMD(z) in a very flexible manner by creating

18

as many atomic groups as desired/required.

19

It is worth emphasizing that the calculation of the NSLDMD in the present approach is

20

independent of the volume fraction calculation. As indicated by Equation 3, the NSLD is

21

obtained using the volume of each layer (Vlayer,i) in which the system has been discretized.


8

Page 9 of 35

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

The calculation of absolute molecular volumes is also based on the construction of 3D Voronoi

2

cells for each atom of the system. However, in this case, no partition over layers is required and

3

the reported value corresponds to the volume (in units of Å3) of the selected group of atoms (see

4

Equation 5) from each frame of the MD trajectory.

5

VolMOL = j

∑ Vk atomsMOL j

Equation 5

k=

6

To account for the boundaries of the system in the tessellation process, the dimensions of the box

7

are passed to Voro++ to define the planes to cut the voronoi cells of those particles close to the

8

box edge for every MD snapshot. This process could incur volume underestimation but will

9

likely be minimized by averaging over MD frames.

10

A challenging aspect of the present approach in the calculation of volumes is the presence of

11

void space in the simulation box, which induces an overestimation of the volume for those

12

particles surrounding the ‘hole’. This is certainly a limitation to the general applicability of this

13

approach. However, when the void space corresponds to a slab of vacuum, e.g. water/vacuum

14

interface, the overestimation is minimized when the system is suitably partitioned instead of

15

using the whole box (tessellation per layer mentioned above).

16

Conversely, 3D Voronoi tessellation provides a consistent measurement of effective molecular

17

volumes that responds to environmental effects and phase changes, while relying only on the

18

position of the particles during the simulation.

19

Finally, the reflectivity profile is calculated using the Abeles matrix method.4 This approach is

20

commonly used in the fitting of experimental neutron scattering data, in packages such as

21

MOTOFIT,9 Reflpak/Refl1D (http://www.ncnr.nist.gov/reflpak/), GenX

22

(http://genx.sourceforge.net), Parratt (http://www.physics.brocku.ca/~tharroun/parratt/), and


9


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 35

1

many others (http://material.fysik.uu.se/Group_members/adrian/reflect.htm#Analysis).

2

Typically, the thickness and NSLD of each layer in a model are parameters optimized during the

3

fitting procedure, starting from an initial guess defined by the user. In NeutronRefTools, this

4

somewhat arbitrary initial guess is replaced by the MD-based structure of the system (Figure

5

1B), from which the reflectivity profiles for different neutron contrasts can be generated by

6

defining the deuteration level of e.g. the solvent or other molecular constituents.

7

To describe the reflectivity calculation, let us consider an incident neutron beam that is

8

refracted by each of the layers considered in Figure 1B. The value of the wavevector, k, in layer

9

n, is given by Equation 6: 2

10

Q k n =   - 4 ⋅ π ⋅ (NSLDn - NSLD0 ) 2

11

where Q is the momentum transfer (Equation 2). From the wavevectors corresponding to two

Equation 6

12

adjacent layers, for example n and n+1, the Fresnel reflection coefficient between those layers

13

can be calculated using Equation 7: k n - k n+1 k n + k n+1

(

exp - 2 ⋅ k n ⋅ k n+1 ⋅ σ n,n+1

2

)

14

rn,n+1 =

Equation 7

15

where a Gaussian modification is applied to account for the roughness/diffuseness of each

16

interface, as described by Nevot and Croce.10 Using the Fresnel coefficient for each interface and

17

a phase factor, βn, which accounts for the thickness of each layer (see Equation 8), a

18

characteristic matrix Cn, is constructed for each layer n, as shown in Equation 9.

19

20

0 ↔ n = 0 βn =   k n ⋅ thickness n ↔ 1 ≤ n ≤ N

Cn = exp(i ⋅ βn 1 )  rn ⋅ exp(- i ⋅ βn

1

)

rn ⋅ exp(i ⋅ βn 1 )  exp(- i ⋅ βn 1 ) 


Equation 8

Equation 9

10

Page 11 of 35

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 2


Considering the product of the characteristic matrices for all the layers, a resultant matrix is obtained (Equation 10): N

3

M = ∏ Cn

Equation 10

n=0

4 5

whose M11 and M21 elements give the reflectivity corresponding to a given Q value (Equation 11):

M 21 ⋅ M *21 M 11 ⋅ M *11

6

R(Q ) =

7

In addition to the graphical comparison between the experimental and theoretical reflectivity

Equation 11

8

profiles, NeutronRefTools also provides a more quantitative comparison by means of the

9

reduced χ2 metric, calculated as indicated in Equation 12:

∑

Nq i =1

( Rexp (Q) − RMD (Q))2 sd exp Nq − 1

2

10

χ =

11

where Nq refers to the number of q values for which the experimental and MD-based

2

Equation 12

12

reflectivities (Rexp(q) and RMD(q), respectively) were measured, and the sdexp corresponds to the

13

experimental standard deviation of each data point. The usage of the reduced χ2 value may be

14

useful in the selection of the parameters during the fitting procedure.

15

To account for the experimental resolution, that is, the relative error associated to the Q values

16

(dq/Q) NeutronRefTools also offers the option to smear the reflectivity profile by means of

17

discrete convolution with Gaussian functions with user-defined FWHM (Full Width at Half

18

Maximum) and number of points for discretizing the Gaussian . In the present approach, the

19

Gaussian FWHM equals the absolute error dq, obtained from dq/Q for each Q value in the


11


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 35

1

reflectivity profile. The value of dq/Q in turn can be a user defined constant, or read directly for

2

each point from the experimental data file.

3

NeutronRefTools implementation of the methods described above allows the user to choose

4

different computing options and to input specific parameters depending on the desired property

5

to be calculated. In the case of the NSLD and molecular volume calculations some of the

6

available options are : i) the model type: AA or CG; ii) the calculation type: NSLD and/or

7

molecular volumes; iii) trajectory-related options i.e. range of frames to be analyzed and/or how

8

many are selected (e.g. one every ten frames), iv) the thickness of the layers in which the

9

simulation box will be split into; v) the percentage of solvent deuteration; vi) the selection of

10

atoms for the decomposition of the NSLD and volume fraction profiles, with their respective

11

deuteration; vii) the selection of atoms to compute the absolute volumes.

12

Regarding the reflectivity calculation, the main options include: i) the Q range to evaluate the

13

reflectivity (input manually or read from experimental data file); ii) parameters such as the

14

incoherent background noise in the experimental data ; iii) the selection of the NSLDMD region

15

to be used, and the parameters of any additional layers, such as substrate material for supported

16

bilayers, to generate an MD-based NSLD of the full system (e.g. the substrate-lipid bilayer-water

17

interface) that represents the experimental setup. For example, to compare a MD simulation of a

18

lipid bilayer in water with a supported bilayer experiment, the NSLDMD is used to generate

19

layers 3 to N, representing the solvated bilayer and, layers 0, 1 and 2 are used to model the

20

supporting surface (i.e. layer 0: Silicon wafer, layer1: SiO2 surface and layer 2: optional layer of

21

solvent between the surface and membrane), and layer N+1 is used to complete the system with

22

an infinitely thick layer of bulk water (see Figure 1). In addition, by removing a suitable number


12

Page 13 of 35

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

of points from the NSLDMD, the amount of water between the bilayer and the supporting surface

2

can be adjusted to match the experimental system.

3

In neutron reflectivity experiments, the method of contrast variation allows to selectively

4

highlight or hide different molecular components of the system by partial or complete

5

deuteration of the components and/or the solvent. Due to the notable difference in the coherent

6

scattering length between hydrogen and deuterium, this approach leads to multiple reflectivity

7

profiles with different intensity for the same physical system. Accordingly, for a given structural

8

model to be considered correct, it must reproduce the reflectivity profiles for all the

9

experimentally generated contrasts/deuteration schemes.

10

NeutronRefTools allows the selection of a percentage of the solvent and/or specific

11

components of the system, to be deuterated and the generation of the corresponding NSLD and

12

reflectivity profiles, which can be compared to the corresponding experimental measurements. It

13

is important to bear in mind that this approach is a first approximation based on assuming

14

equivalent dynamics for deuterated and non-deuterated molecules. Alternatively, multiple

15

simulations, each with the desired deuteration level can be performed and analyzed

16

independently with NeutronRefTools.

17

Molecular Models. The simulation of processes involving surfaces, interfaces or thin films,

18

such as lipid bilayers, are often computationally expensive because of the large system size and

19

long time-scales required to sample the conformations of interest. For this reason, analysis of

20

MD trajectories generated using CG models, in addition to AA simulations, is available. In

21

particular, the current implementation is compatible with the SDK CG force-field.11-13 In the

22

current version, there is no special treatment for the CG beads in the generation of the volume

23

fraction profiles. As in the AA case, the Voronoi cell is generated for each particle, weighted by


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 35

1

the radius (AA or CG), and the entire volume of the cell is assigned to the layer where the center

2

of the bead is located. This is the reason for the choice of a bin-size of 6 Å in the CG volume

3

fraction profile in test case 1 (see Results). An alternative approach, where the cell volume is

4

split between the layer that contains the center of the particle and the two adjacent layers was

5

also considered. In this case, the proportion of volume assigned to each layer depends on the

6

distance of the center of the particle relative to the boundaries of the layer containing it. From

7

such distance, the fraction of the cell volume corresponding to each layer could be estimated

8

assuming a spherical particle. However, as potential artifacts might arise with such assumption,

9

in addition to a higher computational cost, the simpler approach described before was preferred.

10

Test Case 1. Neutron Reflectivity of a DMPC bilayer

11

The first model example consists of generating the NSLD, volume fraction and reflectivity

12

profiles for a DMPC lipid bilayer from MD simulations at both atomistic and coarse-grained

13

levels, and comparison with experimental reflectivity data. For this purpose, a pre-equilibrated

14

lipid bilayer containing 120 DMPC (1,2-Dimyristoyl-sn-glycero-3-phosphocholine) molecules

15

was placed in a rectangular box with the bilayer plane aligned to the x-y plane. The system was

16

solvated using the Solvate plug-in of VMD7 and water molecules overlapping with the lipids

17

were removed. The CHARMM3614 force-field was used for the lipids and the TIP3P model15 for

18

water. Initial steric clashes were removed by 5000 steps of minimization followed by 500 ps of

19

equilibration with restraints on the phosphate groups to maintain the thickness of the bilayer

20

during the relaxation/thermalization of the system. MD simulations were performed in the NpT

21

ensemble, with a semi-isotropic pressure coupling at 1 atm using the Nose-Hoover Langevin

22

piston,16 while the temperature was controlled at 303 K by means of the Langevin thermostat.17

23

Long-range electrostatic interactions were treated using the particle mesh Ewald algorithm,18 and


14

Page 15 of 35

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

van der Waals forces were smoothly switched off, using the vdwForceSwitching option, between

2

10-12 Å. The RATTLE algorithm19 was used to constrain all bonds involving hydrogen atoms, in

3

order to use a 2-fs time-step. The multi-time step algorithm r-RESPA20 was used to integrate the

4

equations of motion. Non-bonded short- and long-range forces were computed every time step.

5

A 120-ns AA simulation was generated using NAMD2.9

6

(http://www.ks.uiuc.edu/Research/namd/).21 Subsequently, the all atom system was mapped into

7

the SDK model13 using in-house Tcl scripts. A CG-MD simulation of 1µs was generated using

8

LAMMPS (http://lammps.sandia.gov).22 The system was evolved in the NpT ensemble using

9

Velocity Verlet time integration with Nose-Hoover thermostat and barostat,23 with a semi-

10

isotropic pressure coupling at 1 atm, and temperature controlled at 303K. A cut-off of 15 Å for

11

non-bonded interactions was used in combination with long-range electrostatic interactions

12

evaluated using the PPPM method.24 A time step of 10 fs was used to integrate the equations of

13

motion.

14

For the experimental counterpart, specular neutron reflection of a supported bilayer of DMPC on

15

a 80mm x 50mm x 15mm silicon (111) wafer was measured in H2O and D2O on the D17

16

reflectometer25 at the Institut Laue Langevin (Grenoble), using neutron wavelengths λ of 2-20 Å

17

and two angles of incidence (0.7 º and 3.0 º). The supported bilayer was formed by in-situ vesicle

18

fusion (0.5mg/ml DMPC SUV in D2O, formed by tip-sonication) after measuring the silicon

19

substrate reflectivity in D2O and H2O. Initial optical matrix modelling of the structure was

20

performed using 3 layers (head-groups, tails, head-groups) to describe the lipid bilayer, which

21

was found to have a total thickness of 41±2 Å, with an average area per molecule of 62±5 Å2,

22

and a polar head group thickness of 7±1 Å. The support surface had a SiO2 layer with a thickness


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

Page 16 of 35

1

of 12±2Å, 24 v/v% solvent and interfacial roughness of 2.5Å (toward the Si) and 2.5Å (toward

2

the solution).

3

Test Case 2. Calculation of molecular volumes for DOPC and DPPC lipids in one- and two-

4

component bilayers at two different temperatures

5

The second application considered is the calculation of lipid volume changes with varying

6

temperature and composition from AA simulations of single and mixed DOPC and DPPC

7

bilayers. To this extent, two lipid bilayers containing 108 DOPC (1,2-Dioleoyl-sn-glycero-3-

8

phosphocholine) and 116 DPPC (1,2-dihexadecanoyl-sn-glycero-3-phosphocholine) molecules

9

were placed in two rectangular boxes with the bilayer plane aligned to the x-y plane.

10

Additionally, a mixed DOPC/DPPC bilayer was built using 56 DOPC and 56 DPPC molecules,

11

placed in a rectangular box and aligned to the x-y plane. The three systems were solvated and

12

equilibration and production runs were performed using the same methods and parameters as

13

those used in the AA simulations of test case 1. Two temperature values were considered (298

14

and 323 K). The centers of mass of the phosphate groups in the DPPC bilayer simulation at 298

15

K were subjected to a sequence of restraints with force constants of: 1, 0.1, 0.01 and 0.001

16

Kcal/mol for eight nanoseconds in the first instance and four nanoseconds in each of the

17

subsequent restraints respectively in order to keep the inter-leaflet distance separation at 43 Å,

18

and to achieve an equilibrated bilayer in the gel phase. Such equilibration was followed by 50 ns

19

of production MD. Equilibrated patches of liquid phase DOPC and gel phase DPPC bilayers

20

were used to generate the mixed DOPC/DPPC system for the simulation at 298 K. This

21

configuration was further equilibrated for 20 ns, and a production run of 50 ns was used for

22

analysis. The area per lipid was calculated using MEMBPLUGIN.26

23

Results


16

Page 17 of 35

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


Test Case 1. Neutron Reflectivity of a DMPC bilayer

2

NeutronRefTools was used to calculate the NSLD, reflectivity and volume fraction profiles

3

from AA and CG MD trajectories of a DMPC bilayer/water system. The MD-based reflectivity

4

data was compared with experimental measurements for a supported DMPC bilayer in water and

5

deuterated water. Additionally, the molecular volumes were also calculated and compared to

6

experimental values reported in the literature. Comparison of the total NSLDMD profiles

7

calculated for both AA and CG MD trajectories for the DMPC bilayer (Figure 2A), and their

8

decomposition into molecular component contributions (Figure 2B), illustrates the potential of

9

the tool in the analysis of both molecular representations. Moreover, in this example, a good

10

agreement between AA and CG simulations is also achieved. As mentioned earlier, generation of

11

multiple contrasts is required to attain a reliable structural model to interpret the experimental

12

reflectivity data. Figure 2A shows the NSLDMD for the DMPC bilayer system in both 100% of

13

the H2O and after replacement of all solvent molecules by D2O. The degree of substitution can

14

be chosen by the user by specifying the percentage of water molecules to randomly substitute by

15

deuterated water. Combining the computed NSLDMD with the definition of parameters for the

16

substrate/solvent layers 0, 1, 2 and N+1, and the instrument background (see Methods section,

17

and Table 1), the reflectivity profiles for the DMPC bilayer system were computed in water and

18

deuterated water. One consideration at this point is that the NSLDMD obtained from simulations

19

of a lipid bilayer using periodic boundary conditions by definition corresponds to a solvent

20

volume fraction at the lipid tails level of zero (except in the presence of water pores). In contrast,

21

in supported bilayer experiments, the supporting surface is usually not perfectly covered by the

22

lipids but patches of water can remain in membrane defects affecting the reflectivity profile and

23

corresponding effective NSLD of the bilayer. This situation can be addressed by correcting the


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 35

1

NSLDMD to include a variable level of solvent volume fraction at the lipid tails level, until good

2

fit between MD and experimental reflectivity is obtained. An initial guess to the amount of

3

water at the lipid tails can be estimated from the difference in MD and experimental area per

4

lipid, when the latter corresponds to the area occupied by both solvent and lipids. In the present

5

example, the experimental area per lipid is 62 ± 5 Å2 while the values measured from the MD

6

trajectories are AA: 61 ± 2 Å2 and CG: 60 ± 1 Å2. This indicates that the initial guess would be

7

in the order of the experimental relative error in the area per lipid (~ 9 %).

8 9

In line with the previous consideration, the raw reflectivity profiles generated from the AA MD trajectories (Figure 3A) indicate a disagreement between the experimental and MD data

10

(dotted traces), particularly clear in the D2O contrast. Besides the water content at the lipid tails

11

region, an alternative (or complementary) explanation for this inconsistency could also be a

12

difference in the water content at the head-group level. This could be a consequence of the

13

contact between the bilayer and the supporting solid surface in the experiment, while such an

14

interaction is absent in the MD simulation. To test this hypothesis, the NSLDMD and associated

15

reflectivity profiles were tuned to consider variations in the solvent volume fraction at the head-

16

groups and lipid tails levels (see Table 1). Using this approach, a corrected MD-based reflectivity

17

profile, showing a reasonable fit to the experimental data, was obtained (Figure 3A, continuous

18

traces). This result indicates that the water content in supported and free bilayers differ at the

19

head-group level, supporting the previous hypothesis. However, it also suggests that the MD data

20

of a free bilayer in solution can in fact be used for the fitting of supported bilayer neutron

21

reflectivity data, as long as the water content at the head-groups and lipid tails levels is properly

22

considered. The suitability of this approach is also evidenced when comparing the water content


18

Page 19 of 35

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

corrected reflectivity profiles obtained from CG-MD simulations with the experimental data

2

(Figure 3B).

3

Although reasonable good fits have been obtained using the present approach, some residual

4

disagreement between experimental and MD reflectivity profiles still remains. Besides the

5

aforementioned differences in the lipid bilayer water content, another source of error could reside

6

in the absence of an explicit supporting surface in the molecular simulations. In other words,

7

mixing free bilayer molecular simulations data with artificial layers representing the supporting

8

material, might fail to mimic subtle structural constraints induced on the lipids by the surface,

9

which would translate into small differences in the NSLD and reflectivity profiles. However, the

10

present approach has the advantage of being independent from the availability/development of

11

suitable force field representations of the supporting surface materials.

12

Finally, the results obtained for both AA and CG models not only illustrate the applicability of

13

the method but also highlight its potential use in testing different molecular representations and

14

force fields for molecular simulations. In addition, the results point to a distinct difference in the

15

head-group hydration in free and solid-supported bilayers. It has long been debated whether the

16

solid support imposes constraints on the molecular conformation, mobility or dynamics of the

17

lipids, albeit without conclusive evidence.

18

The volume fraction calculation using NeutronRefTools relies exclusively on the position and

19

radii of the particles. From this data, the radical Voronoi tessellation of Voro++8 is used to

20

construct the Voronoi cell around each atom or CG bead in the system, and obtain the associated

21

particle volume for each simulation frame. Such per particle data is then gathered by molecular

22

component, defined by the user, and plotted along the z axis. Figure 4 shows the output

23

corresponding to the DMPC/water system for both AA and CG MD data. While the calculation


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 35

1

for the AA system was done using layers of 2 Å thicknesses, the higher granularity of the CG

2

model required layer thicknesses of 6 Å to accommodate the bulkier beads. Consequently, in

3

general, fewer data points are generated in the CG volume fraction profiles. For this reason,

4

intermediate data points were interpolated between the CG measurements to facilitate graphical

5

comparison. Despite the different level of detail between the AA and CG representations,

6

generally, good agreement is observed, particularly in the regions containing pure water (z < 7 Å,

7

z > 57 Å) or pure lipids (15 Å

Molecular Dynamics Simulations and Neutron ... - ACS Publications

Recommend Documents