Article pubs.acs.org/JPCB
Combined Coarse-Grained Molecular Dynamics and Neutron Reflectivity Characterization of Supported Lipid Membranes Alexandros Koutsioubas* Jülich Centre for Neutron Science (JCNS) at Heinz Maier-Leibnitz Zentrum (MLZ), Forschungszentrum Jülich GmbH, Lichtenbergstr. 1, 85748 Garching, Germany S Supporting Information *
ABSTRACT: Supported lipid bilayers on planar surfaces constitute an archetypical experimental system for the study of biological membranes. The popularity of these ordered molecular layers in the literature, is on one hand related to the simplicity of their preparation using the method of vesicle fusion and on the other hand to their compatibility with a multitude of surface sensitive experimental probes. Neutron reflectivity has proven as an important experimental method for the investigation of such systems with the ability to provide subnanometer structural information perpendicular to the supporting plane. Traditionally reflectivity data are compared to theoretical curves of simplified models consisting of stratified layers representing the hydrophilic (lipid heads) and hydrophobic (lipid tails) parts of the bilayer. In the present work we explore the combined use of molecular simulations and neutron reflectivity for the characterization of supported membranes. By performing coarse-grained molecular dynamics simulations based on the MARTINI force field of supported 1,2dihexadecanoyl-sn-glycero-3-phosphocholine (DPPC) bilayers close to a hydrophilic substrate, we compared the obtained reflectivity profiles with neutron reflectivity data for this system at a series of temperatures above and below the main phase transition. It is found that the use of an imperfectly smooth substrate in the coarse grained simulation is of vital importance for avoiding the artificial freezing of water that is trapped between the surface and the bilayer. The observed quantitative agreement between simulation and experiment using “rough” supporting surfaces, especially for the liquid lipid phase, exhibits that the presented methodology may serve as a basis for the detailed and assumption-free investigation of more elaborate systems.
■
INTRODUCTION Biological lipid membranes represent an important selfassembly system in living organisms, that perform the function of separation and compartmentalization by acting as selective permeable barriers. In their most abundant form as the primary constituent of the cell membrane, they provide a bilayer fluid matrix for embedded membrane proteins that in turn regulate the cell’s communication with the surrounding environment and its transport of chemical species. Due to the immense complexity of membranes in living cells, model lipid bilayer membranes with well-defined composition serve as relatively simpler counterparts in experimental studies. In particular, supported lipid bilayer membranes separated from planar substrates by nanometer-thick water layers1 have found widespread use in the literature due to their simple and reproducible preparation through lipid-vesicle fusion2 or monolayer transfer (Langmuir−Blodgett technique) and also due to the fact that many thermodynamic and structural properties of free-standing bilayers are maintained thus providing a versatile biomimetic environment. The special nature of supported membranes (i.e., their reconstitution on a planar surface) enables the application of many different surface-sensitive experimental techniques3,4 © 2016 American Chemical Society
(atomic force microscopy, quartz crystal microbalance, evanescent wave methods, fluorescence microscopy, neutron/ X-ray reflectivity) for the study of their structure in relation to various physical parameters and/or under the presence of membrane active biological molecules. Among these methods, neutron reflectivity that measures the thickness and scattering length density (sld) profile of interface structures, presents several advantages for the study of supported membranes. Notably the highly penetrative nature of neutrons that renders the technique ideal for the study of buried interfaces and the exploitation of high contrast between the scattering length of hydrogen and its heavy isotope deuterium that enables the application of contrast matching techniques, have led to structural insights concerning supported lipid membranes at the subnanometer level.5−14 In neutron reflectivity studies of supported bilayers the specular reflection of a solid/liquid interface supporting the membrane is probed and related to the scattering length density profile normal to the interface which contains direct Received: May 30, 2016 Revised: October 14, 2016 Published: October 17, 2016 11474
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B
supramolecular structures. In this respect the present work is focused on the comparison of neutron reflectivity data of DPPC bilayers formed on a hydrophilic surface through vesicle fusion, with coarse-grained MD trajectories of the aforementioned system at temperatures both above and below the main phase transition. In all simulations we have adopted the MARTINI force-field35 and evaluated the effects of using a polarizable or nonpolarizable water model, and most importantly the effects of the structure of the surface. It is shown that by preparing a “rough” hydrophilic substrate and by carefully tuning the thickness of the trapped water layer below the lower membrane leaflet, no artificial crystallization of water takes place while the structure of the supported DPPC membranes in the gel and liquid phase is in very good agreement with the acquired neutron data.
information about the molecular distribution within the membrane. Traditionally the sld profile is modeled as a series of stratified layers representing the contribution of substrate, hydrophilic lipid heads, hydrophobic lipid tails, and aqueous solvent. Then the corresponding reflectivity is calculated using the appropriate mathematical formalism (Parratt recursion15 or Abeles matrix method16). Minimization of the discrepancy between the model calculated reflectivity and experiment, leads to an estimation of the model’s parameters. The acquisition of reflectivity curves at multiple solvent contrasts (H2O/D2O ratios) or the manipulation of the lipid contrast via isotopic substitution gives a data set of reflectivity curves that under simultaneous model fitting leads to an overall decrease of the ambiguity of the determined system’s parameters. In practice this approach is quite efficient in the case where the interaction of supported membranes with biological molecules (proteins and small peptides)17−19 or nanoparticles20,21 is studied. However, as complexity increases or subtle structural features need to be resolved, the comparison of experimental reflectivity with molecular simulation data might prove beneficial for filtering out candidate models.22 Only recently attempts to explicitly couple neutron reflectivity data with molecular dynamics simulations of free-standing lipid membranes appeared in the literature,23 where it is shown that except for some residual discrepancies (probably due to the absence of a supporting surface), experimental reflectivity curves of supported DMPC membranes are in reasonable good agreement with those corresponding to simulated free-standing membranes. Comparing to free-standing bilayers, supported ones are known to present a large suppression of their characteristic undulations24 and are also characterized by a certain amount of asymmetry between their two leaflets,25,26 due to the interaction of the bottom one with the substrate. While there exists a wide range of molecular dynamics simulation studies of free-standing lipid membranes,27 there are far fewer allatom28,29 or coarse-grained30−34 simulation works dedicated to supported lipid bilayers, a fact related to the added complexity in the presence of a solid support in close proximity to the membrane. An important insight obtained from these studies is related to the substrate induced asymmetry between the proximal and distal leaflets that is manifested both in the membrane’s structural (density,ordering) and dynamical (lipid diffusion) properties. In cases where a coarse-grained representation of the system is considered, water models (as in the MARTINI force field35) tend to artificially freeze and produce layering effects near solid substrates, thus introducing important artifacts in the simulation. Efforts to avoid such problems by modifying the thermodynamic parameters of the model has been criticized recently,36 since they lead to an appreciable modification of the equilibrium area per lipid and thus affecting the whole set of membrane properties. Given that comparing to all-atom simulations, coarse-grained models are able to vastly reduce the number of degrees of freedom by grouping many atoms into effective interaction sites and consequently access wider length and time scales, it would be interesting to investigate the conditions under which coarsegrained simulations of lipid membranes near a solid support could quantitatively reproduce relevant neutron reflectivity experimental data. Success in this direction would permit the combined use of neutron reflectivity and MD simulations for detailed studies of more complicated systems involving the interaction of membranes with biological molecules and
■
MATERIALS AND METHODS
Neutron Reflectivity. Neutrons are characterized by a weak interaction with matter and thus have the ability to penetrate deeply into many materials, a fact that makes them ideal probes for buried interfaces. In particular neutron reflection where the neutron beam penetrates through a solid support and then encounters a solid/liquid interface, has been used extensively for structural studies of various types of supported bilayer membranes.5−12 Specular neutron reflectivity measures the thickness and scattering length (sld) density profile of interface structures along the surface normal (z). Scattering length density ρ(z) depends on the chemical composition of materials and is defined as the sum of the coherent nuclear scattering lengths (bi) of its constituent atoms and their number density (ni(z)) so that ρ(z) = ∑ibini(z). Reflectivity data are gathered by recording the intensity of reflected neutrons relative to the incident beam as a function of the momentum transfer vector (q = 4πsin θ/λ), where θ is the incidence angle and λ the wavelength of incident neutrons. The variation of reflectivity as a function of momentum transfer R(q) is related to the square modulus of the one-dimensional Fourier transform of the sld profile (ρ(q)) through the relation P(q) = (16π2/q2)|ρ(q)|2. Recovery of the sld profile from P(q) is performed usually by model fitting using thin slab models. Materials and Experimental Procedures. Lipid molecules 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) were purchased from Avanti Polar Lipids in the form of lyophilized powder and were used without any further purification. Deuterium oxide (D2O) of 99.8% purity, was purchased from ARMAR (Europa) GmbH. Unilamellar vesicles (ULVs) were prepared by sonication as follows. Weighted amounts of DPPC powder were dissolved in small volumes (2−3 mL) of chloroform and the solution was transferred to a small round-bottom flask. Chloroform was removed using a rotary vacuum evaporator, followed by overnight drying in vacuum. The dry lipid films at the wall of the flask were hydrated with a D2O buffer (150 mM: NaCl, 5 mM: HEPES, pD = 7.4) preheated to 50 °C, followed by vigorous magnetic stirring for 1 h at 50 °C. Then the resulting multilamellar vesicle suspension was subjected to probe tip sonication (Sonics Vibra-Cell VCX130) for several minutes, until the initial turbid solution became almost clear. Final solution concentration was in the range of 2 mg/mL. All ULV solutions were used for measurement within 1 h of preparation. The size of the ULVs was characterized by dynamic light 11475
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B
For the model fitting procedure, each phospholipid bilayer is represented by a 3-layer stack corresponding to their lipid heads/lipid tails/lipid heads structure, where the two lipid head layers are identical. The molecular sld values that were used for each layer are summarized in Table 1. Layer roughness is
scattering (DLS) using a MALVERN Zetasizer Nano-S instrument. Neutron reflectivity data were acquired on the MARIA vertical reflectometer37 of the JCNS outstation at MLZ in Garching (Germany), using temperature-regulated (through a connected Julabo F12-ED circulator) custom liquid cells. The design of the liquid cells is similar as the one used in a previous work,38 with the difference that the main body of the cells was made from boron-glass in order to minimize background. The measurements were performed using two different wavelengths 12 Å for the low-q region and 6 Å for the high-q region up to 0.25 Å−1, with a wavelength spread Δλ/λ = 0.1. Ultrapolished Si blocks (r.m.s. roughness of 1−2 Å, dimensions 150 × 50 × 20 mm3) purchased from Andrea Holm GmbH were used as substrates and were cleaned before each experiment by soaking them in a 5:1:1 deionized water/ 29%, NH4OH/30%, H2O2 solution maintained at 75 °C for 15 min, followed by a rinse with milli-Q water. Since the Si/liquid interface is probed, the collimated neutron beam enters through the one side of the Si block, is then reflected at the solid/liquid interface and finally exits through the other face of the Si block toward the detector. The combination of high-flux at the MARIA instrument and low-background cells, permitted the detection of reflectivities down to 10−6−10−7 with total acquisition times of about 2 h for a single curve. Data Analysis. The acquired neutron reflectivity profiles were analyzed either by classic model fitting using a custom simulated annealing minimization algorithm written in Fortran 90 or by comparison with theoretical reflectivity curves corresponding to the sld profiles obtained by the coarse grained simulations. Theoretical curves were calculated using the Abeles matrix formalism16 where each layer n is defined by three parameters, i.e., scattering length density ρn, thickness dn, and roughness σn. In more detail, the Fresnel reflection coefficient between layer n and n + 1 is given by rn,n + 1 =
kn − kn+1 2 exp( − 2k nk n + 1σn,n + 1) kn + kn+1
Table 1. Molecular Scattering Length Density Values Used in This Study
(1)
1/2
⎡ exp(k ndn) rn exp(k ndn) ⎤ ⎥ cn = ⎢ ⎢⎣ rn exp( −k ndn) exp( −k ndn)⎥⎦
cn
n=0
(3)
(4)
Finite instrumental resolution δq/q is taken into account by assuming a Gaussian resolution function that is convoluted to the theoretical curves using the following equation for performing a p point average p ⎛ a δq ⎞ waR ⎜q + ⎟ / ∑ wa p q ⎠ a =−p ⎝ a =−p p
R′(q) =
∑
(5)
with wa being the Gaussian weight
wa = exp( −2(a /p)2 )
∑ ∑ i = 1 ,.., M
k
nk b k /M Vlayer
(7)
where nk is the number of coarse-grained beads of type k found in layer j, bk is the coherent scattering length for the coarsegrained bead type (see Table 2), and Vlayer is the volume of the layer defined by the xy surface of the simulation box times the bin-size in the z-direction (4 Å). In turn, the reflectivity curves corresponding to the MD simulation are calculated using eqs 1−6 for the found average sld profile. Before the final comparison of the calculated reflectivity curves with the experimental data, two small corrections have to be applied. The first is related to the imperfect coverage of the substrate by the lipid membrane so that patches of solvent (∼5% in the present case) have to be taken into account.23 The second correction is related to the bulk sld value of the solvent that needs to be slightly adjusted due to the apparent density of the coarse-grained water model that may differ from the density of D2O. Coarse-Grained Molecular Dynamics. All simulations were performed using GROMACS v5.1.1 simulation suite41−44 and the coarse-grained (CG) model MARTINI v. 2.0 proposed
from which finally the reflectivity is calculated as
R(q) = |M11/M 21|2
2.07 3.47 1.73 −0.37 6.22 6.32 −0.55
j
(2)
# of layers
∏
Si SiO2 DPPC heads DPPC tails D2O buffer @ 333 K D2O buffer @ 293 K H2O buffer @ 333 K
ρlayer =
and the system’s matrix M is given by M=
sld ( × 10−6 Å−2)
introduced by applying Névot-Croce correction terms to the Fresnel coefficients (eq 1).39 Additional “hydration” layers at the interface between bilayer and substrate are added according to the needs of the fit. The reported fitted parameter uncertainties were estimated through calculation of each parameter’s standard deviation after executing multiple simulated annealing minimizations. In order to reduce the ambiguity of the fits we have tried to reduce the number of free parameters through simple arguments based on the molecular structure of lipids. Since lipid headgroups and tailgroups are chemically connected, the area (A) per head and tail should be the same. Given the molecular volumes40 Vh = 326.3 Å3 and Vt = 889.2 Å3 occupied by the two lipid species, we can relate the head and tail region sld values21 so that we may reduce the required parameters by one. For the comparison between simulation and experiment the average sld profile perpendicular to the supporting surface had to be calculated from the MD trajectories. For this purpose several snapshots (i) (during the MD simulation) of the system are discretized (binned) along the z-axis in (j) layers of equal thickness (d = 4 Å) and the sld of each one is calculated as
where kn = [q /4−4π(ρn − ρ0)] is the neutron wavevector in layer n. For each layer a characteristic matrix is defined as 2
layer
(6) 11476
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B
controlled during simulation using a V-rescale thermostat with a temperature scaling factor of 1.5 ps. The cut off for nonbonded Coulombic interaction was rc = 1.2 nm and the force was smoothly shifted to zero at cutoff. The neighbor list for nonbonded interactions was updated every 10 steps. Periodic boundary conditions in the xy plane and a time step of 30 fs for simulations with the standard water model and of 20 fs for the polarizable water model was used. For all systems 3 to 5, 1 μs-MD simulations were performed, which are equivalent to about 4 μs of effective time.35 For the calculation of reflectivity curves only the last 500 ns of the MD trajectory were considered by averaging 500 system snapshots (one snapshot every 1 ns). For the same system parameters, the calculated reflectivity curves show no appreciable variation. Finally a set of complementary simulations (in the NPT ensemble) of free-standing nonsupported membranes was performed for comparison purposes.
Table 2. Coherent Scattering Length (bk) for Each MARTINI CG Bead Typea MARTINI CG bead type
atoms
bk (fm)
Q0(choline) Qa(phosphate) Na(glycerol linkage 1) Na(glycerol linkage 2) C1(first three hydrocarbon beads) C1(last hydrocarbon bead) P4 or PW(4D2O molecules)
C5NH13 PO4 C3O2H3 C2O2H2 C4H8 C3H7 D8O4
−6.01 28.33 31.55 17.42 −3.32 −6.23 76.56
a
Note that for the two glycerol linkages, two different values for the same bead type have to be considered due to their different stoichiometry after coarse-graining. The same applies for the last hydrocarbon bead type that contains one carbon and one hydrogen atom less than the rest.
■
by Marrink.35 Phospholipid (DPPC) molecules are represented by 12 CG interaction sites (beads), with 2 charged (type Q0 and Qa) beads representing the phosphate-choline headgroup, 2 (type Na) hydrophilic beads representing the glycerol linkages and 8 (type C1) hydrophobic beads, 4 per each hydrocarbon tail. Water is modeled either by single hydrophilic beads (type P4) representing four real water molecules (standard model) or by three interaction sites (neutral, negatively and positively charged) representing again four water molecules, where the orientational polarizability of water is reproduced (polarizable model).45 With the exception of negative and positive interaction sites of the polarizable water model, all other CG beads interact pairwise via a Lennard-Jones potential with the interaction strengths defined in the MARTINI model.35 Initial configurations for each simulation where constructed by first preparing a square D × D (D = 156 Å) solid support consisting of completely fixed hydrophilic beads (type Nda) at the bottom of the simulation box which had a height of 200 Å. For a limited set of system parameters, tests were also conducted with D = 78 Å and D = 312 Å system sizes in order to confirm the absence of size effects in the production runs with D = 156 Å. Two types of supports where considered (a) a perfectly smooth surface where beads were placed on a grid with a spacing of 3 Å and (b) a “rough” surface where bead positions were those of a thin (20 Å) layer extracted from an independent simulation of bulk water CG beads at T = 310 K.46 In more detail, a separate simulation of a CG-water filled periodic 156 × 156 × 200 Å3 box was performed for several ns. From the final system state the positions of water beads within a 156 × 156 × 20 Å3 slice are assigned to hydrophilic (Nda) beads that are made immobile and are placed at the bottom of the simulation box, serving as a “rough” substrate. Using the insane script47 we have prepared initial DPPC membrane configurations containing N lipids, where the area per lipid was equal to 2D2/N. These nonequilibrated membranes were placed at different distances from the hydrophilic surface in order to test various water layer thicknesses between the membrane and the substrate. Then the system was solvated and relaxed by performing energy minimization. During molecular dynamics simulations, system area was kept constant along the xy plane (NAPzT ensemble) and a constant pressure was kept fluctuating around 1 bar along the zaxis, controlled using a Berendsen barostat and by setting the compressibility equal to 5 × 106 bar−1. Temperature was
RESULTS Model Fitting of Bilayer Formation and Phase Transition Data. After cleaning the silicon block in order to render the surface hydrophilic, we first characterize the native SiO2 thin layer on the Si surface by acquiring a neutron reflectivity curve of the Si/D2O interface. The fit of the curve gave a layer sld ρSiO2 = 3.48 × 10−6 Å−2 and thickness dSiO2 = 9.7 Å. These values were held fixed during the subsequent fits. Following this initial step, we set the liquid cell temperature at T = 60 °C and inject the preheated DPPC SUV dispersion, bringing final lipid concentration to ≈1 mg/mL in the cell. The choice of this elevated temperature above the lipid transition temperature (Tm) ensures that the vast majority of the initially adsorbed vesicles on SiO2 will quickly rupture and form a smooth supported bilayer.48 The vesicles were allowed to fuse on the hydrophilic surface for 60 min (a time that is much larger than the typical fusion kinetics5,48) and then the cell was flushed with preheated D2O buffer in order to remove any unfused vesicles. Acquisition of reflectivity curves before and after supported DPPC bilayer formation (Figure 1) clearly evidence the passage from a quasi-Fresnel R(q) ∝ q−4 behavior to the appearance of
Figure 1. Neutron reflectivity of the silicon/D2O interface before and after the formation of supported DPPC bilayer by vesicle fusion. The temperature of the solution in contact with the substrate was held at T = 60 °C in both measurements. Full lines represent model fits of the experimental data using the methods described in the Data Analysis section. 11477
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B
and liquid crystalline phase are compatible with previously reported data.49 Liquid Membrane Phase Molecular Dynamics. As described in the previous section, from the model fitting analysis of the reflectivity curves corresponding to the liquid state of the DPPC membrane it was found that the area per lipid is equal to 62 Å2. So for the fluid-state MD simulations, membranes composed of 784 lipids were prepared at various distances from the substrate. In this way different amounts of trapped water (in the range of 1.5−5 water CG beads per lower leaflet lipid) between the substrate and the membrane were tested. Initial trials began with MD runs of a perfectly smooth substrate and the nonpolarizable water model (“smooth” nonPW) at the 300−330 K temperature range which is above the DPPC transition temperature of 295 ± 5 K that was found previously using the MARTINI force-field.50 Not surprisingly we observed that the substrate acted as a nucleation center for the water beads,35 leading to quick freezing (crystallization) of the whole system. In order avoid this problem we have tested two alternatives (a) the use of the polar water model (“smooth” PW) and (b) the weakening of water−water nonbonded interaction to 76% of its original value30,31,51 (“smooth” non-PW 76%). For these two systems, the CG water beads do not freeze in the bulk, however in the trapped water layer near the substrate water layering and depletion is clearly observed. Additionally with the weak-water model the membrane expands to a new equilibrium area per lipid31 equal to 76 Å2. This overall change of membrane properties leads to a much higher disagreement36 with the experimental reflectivity curve, comparing to the results with the polarized-water model (Figure 3)
characteristic Kiessig fringes related to the presence of a bilayer membrane and especially to its high contrast hydrocarbon core relative to D2O. The position of the fringe is quite sensitive to the overall membrane thickness and is expected to displace to lower-q for larger membrane thicknesses. Model fitting with the three-slab model and an additional thin hydration layer of about 3−4 Å between the membrane and the substrate gave an overall bilayer thickness of 47.7 Å, with layer roughness of ≈4 Å. This thickness value is fully compatible with results reported in the literature5,10 for DPPC layers in the liquid state. The combination of fitted sld values and thicknesses of head and tail regions indicate that almost full surface coverage is achieved (≈ 95%), while the area occupied by each lipid is equal to 62 Å2. Following the characterization of the DPPC layer after its deposition at T = 60 °C, we focus on the effects of temperature on the structure of the supported membrane. It is well established that at a given temperature a lipid bilayer can exist either in the gel (solid) or liquid (liquid-crystalline) state. This phase behavior is dominated by van der Waals attractive interactions between lipid tails, and is affected by the length and degree of hydrocarbon tail saturation. The two phase states are associated with large differences in the overall fluidity and molecular packing. Reflectivity curves were acquired at 7 different temperatures (T = 20, 30, 35, 40, 45, 50, 60 °C), placing the DPPC transition temperature Tm = 40.9 °C at the center of the experimental window. All curves together with their corresponding fits are plotted in Figure 2. (See Supporting Information for reflectivity
Figure 2. Reflectivity curves and model fitting results (full lines) corresponding to the supported DPPC bilayer at different temperatures above and below Tm. The same set of data in Rq4 vs q representation can be found in Figure S2. In the inset the obtained total and tail region thickness of the bilayer against temperature are plotted and the full lines serve as guides to the eye.
Figure 3. Experimental neutron reflectivity results (in Rq4 vs q representation where curve oscillations due to the presence of the membrane at the interface are more easily visible) for a supported DPPC membrane in the liquid phase (T = 60°C), compared to the theoretical curves obtained from coarse-grained MD trajectories. Results from three different simulation system setups are displayed. (a) Rough substrate with polarizable water model (“rough” PW); (b) a perfectly smooth surface with the polarizable water model (“smooth” PW); and (c) a perfectly smooth surface with the nonpolarizable water model and water−water nonbonded interactions reduced to 76% or their original strength (“smooth” non-PW 76%). In all cases the “trapped” water layer is the one corresponding to 3.6 CG water beads per lower leaflet lipid. In the inset the sld profile of the different simulations are depicted. Note that when a smooth surface is used, sld largely fluctuates near the surface due to the layering and depletion of water.
curves in Rq4 vs q representation.) The most characteristic change as we decrease the temperature, is the shift of the frindge to lower-q, a clear sign of layer-swelling. From the curve fits we get an increase of the total bilayer thickness from 47.7 to 57.8 Å and of the hydrocarbon core thickness from 28.7 to 35.6 Å (see Figure S1 for the sld profile evolution and also the inset of Figure 2). We also note that in the gel phase, the fits gave a reduction of the area per lipid A = 49 Å2 which is explained by the more compact and ordered arrangement of lipid tails. We also note that both derived values of molecular area in the gel 11478
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B Since a smooth surface seems to induce artifacts in the simulation, leading also to relatively poor agreement with experimental reflectivity results (see Figure S3 and S4), we prepared a nonsmooth (rough) substrate in the form of a thin rectangular layer (156 Å × 156 Å × 20 Å) where the position of each bead was taken from a snapshot of bulk CG water simulation (see Figure 4). As shown in a previous work,46 such
Figure 5. Mass density profile of lipid heads, lipid tails, and water for a supported membrane system (top) and for a free-standing membrane (bottom). The free-standing membrane profile was shifted so that the hydrocarbon core centers of the two profiles are at the same position. The shaded area in the top figure represents the substrate.
0 corresponds to random orientation, and the value P2 = −0.5 corresponds to anti-alignment.) The dynamics of the lipid molecules were examined by calculating the lateral diffusion coefficient Dxy of each separate leaflet, since it is expected that due to the interaction of the lower leaflet of the supported membranes with the substrate, its dynamics should be slower comparing to the upper leaflet that faces the bulk solution.52 We have found that in the freestanding system, Dxy−free = (2.6 ± 0.1) × 10−7cm2/s which is quite close to reported experimental values of about (1.4 ± 0.05) × 10−7cm2/s.52 In the supported bilayer, both leaflets are slowed down with the lower leaflet being more affected, Dxy−lower = (1.2 ± 0.05) × 10−7cm2/s and Dxy−upper = (1.9 ± 0.05) × 10−7cm2/s (see Figure S8 and S9). Gel Membrane Phase Molecular Dynamics. For simulations of the DPPC membranes in the gel state, we followed the same protocol as in the case of the liquid state simulations described above and by running the MD simulations at a lower temperature range 270−290 K (below the MARTINI transition temperature for DPPC of 295 ± 5 K50). Taking into account the model-fitting estimation of the area per lipid in the gel state (49 Å2), we prepared membranes composed of 992 DPPC lipids at varying distances from the 156 Å × 156 Å × 20 Å “rough” substrate which is composed of fixed hydrophilic Nda MARTINI particles. Test runs with both the polarizable and nonpolarizable water model were performed. Initial observations showed that the use of nonpolarizable water leads to complete crystallization (freezing) of water in the system, so all subsequent attempts were performed with the polarizable water model. Note that decreasing the temperature down to 270 K led to no observable layering effects with the polarizable water model. It is interesting that the simulation run that presented the best agreement with the experimental data (Figure 6) is the one having the same number of CG water beads trapped beneath the membrane as in the liquid-state simulations that also provided the best agreement with experiment. Except from the evident swelling of the membrane layer in the gel-state (see Figure S5), further analysis of the MD trajectories showed that the average lipid order parameter reached the value P2 ≈ 0.79 for both leaflets, indicating the alignment of the majority of the hydrocarbon tails along the axis normal to the substrate. Note
Figure 4. Typical MD snapshot after 1 μs simulation time of a DPPC membrane in the liquid state near a “rough” substrate (using nonpolarizable water model).
setups prevent water CG bead freezing near the surface. Using these “rough” substrates in the MD simulations with both polar and nonpolar water models (“rough” PW and non-PW, see Figure S5), we were able to obtain quite good agreement with the experimental data when the trapped water layer was the one corresponding to about 3.6 CG water beads per lower leaflet lipid (see Figure 3 and Figure S6). Agreement of comparable quality was also observed with reflectivity data in H2O for the same system that were acquired after solvent exchange in the cell (see Figure S7). It is interesting to note that by varying the hydrophilicity of the substrate (by changing the substrate MARTINI bead type to Na or N0) and by keeping all other simulation parameters fixed, no appreciable change in the calculated reflectivity curves was observed. For this specific system that reproduces the measured reflectivity and shows no signs of artificial CG bead layering in the trapped-water layer, further analysis was performed in relation to its structure and dynamics. As can be seen in Figure 5 the density distribution of various species (lipid heads, tails, water) within the system is not significantly altered from the one found from a CG MD simulation of a free-standing DPPC membrane having the same area per lipid. However, density distributions in the free-standing case appear broader and tend to overlap more with each other, something that has mainly to do with the suppression of characteristic membrane undulations in the supported membrane. Note that lipid heads populate the entire “trapped”-water layer and interact with the attractive hydrophilic substrate. The calculated average lipid order parameter P2 for the free-standing membrane was found equal to 0.425 and increased slightly to the value 0.44 for the supported membrane, where the upper leaflet shows only marginally higher ordering comparing to the lower leaflet. 1 (Note: P2 = 2 (3⟨cos(θ )2 ⟩ − 1) where θ is the angle between the bonds of two adjacent hydrophobic beads and the bilayer normal, the value P2 = 1 corresponds to perfect alignment, P2 = 11479
DOI: 10.1021/acs.jpcb.6b05433 J. Phys. Chem. B 2016, 120, 11474−11483
Article
The Journal of Physical Chemistry B
Consequently for the setup of our CG MD simulations that were finally compared against the acquired experimental curves we have used the standard MARTINI force-field parameters and all initial membrane configurations where prepared near specially prepared rough hydrophilic substrates. By carefully tuning the trapped water layer composition (3.6 CG waters or equally ≈14.5 real water molecules per lower leaflet lipid in the liquid phase), we were able to obtain agreement between simulation and experiment as far as the molecular distribution within the membrane layer is concerned, while the obtained mass density profiles (Figure S10) are in agreement with those reported in previous simulation53 and experimental54,55 studies. In these simulations the thickness of the trapped water layer can be estimated by finding the average distance of choline groups from the substrate which is equal to about 6.5 Å. In simulations at temperatures above the DPPC transition (300− 330 K) where the membranes are in the liquid phase, both the nonpolar and polar water model gave identical results concerning the structure of the lipid membrane and the behavior or trapped water between membrane and substrate. However, this is not the case for simulations at lower temperatures (270−290 K) below the liquid-gel transition at 295 ± 5 K in the MARTINI model.50 In this range of temperatures only the use of the polarizable water model leads to MD trajectories that can be compared with experiment, although somewhat larger discrepancies are observed comparing to the liquid-phase results. The observed semiquantitative agreement between experiment and simulation for DPPC membranes in the gel state probably has also to do with the absence of tilt in the simulated lipid tails with respect to to the bilayer normal. It is well-known that DPPC forms a tilted gel phase (of about 30°)56 due to the imbalance between the headgroup size and the optimal packing distance of the lipid tails. The inability to observe the tilted phase in our simulations is related to the fixed size of lipid tail beads in the MARTINI model.50 In comparison to free-standing membranes, the presence of the supporting surface beneath the membrane appears to affect mildly the structure and dynamics of the whole bilayer in the fluid phase, with the bottom leaflet appearing about 40% less fluid due to its interaction with the attractive substrate (Figure S8), a trend that is observed for all simulated membranesubstrate separation distances (Figure S9). In contrast, by analyzing MD trajectories from membranes near a smooth supporting surface, we found that both the structure and dynamics of the lower leaflet are severely altered. In these simulations, the hydrocarbon beads follow a similar layering trend as the solvent (see Figure S11) while the lateral diffusion coefficient of the lower leaflet is about an order of magnitude smaller comparing to the diffusion coefficient of the upper leaflet (Figure S12 and S13). Presumably apart from the direct interaction of lipid heads with the substrate, these effects are related to the lateral diffusion ability of water in the bulk and at the interface as suggested in a previous simulation work.34 It has been shown experimentally25 that the qualitative behavior of the relative lipid mobility in the two leaflets of supported membranes varies significantly as a function of membrane preparation method, chemical nature, and roughness of the supporting substrate and salt content of the solvent. In this respect the observed mild slowing down of the lower leaflet in the simulations that agree with our neutron reflectivity results, may be a specific feature of zwitterionic head lipid membranes supported on silicon substrates. Nevertheless, the
Figure 6. Experimental neutron reflectivity results (in Rq4 vs q representation) for a supported DPPC membrane in the gel phase (T = 20°C), compared to the theoretical curves obtained from coarsegrained MD trajectories. Results from simulation system setup with a rough substrate and the polarizable water model (“rough” PW). In the inset the corresponding sld profile is depicted.
that the ordering of chains in the free-standing membrane was found to be slightly higher (P2 = 0.82). Adop ting χ as a figure of merit (defined as N
χ=
∑i =points 1
(R exp(q) − RMD(q))2 2 σexp
Npoints − 1
, where Npoints is the total number
of acquired experimental points and σexp the standard deviation of each data point), we noticed that the agreement between simulation derived reflectivity and experimental results is slightly lower comparing to the liquid-state case (χliquid ≈ 3.5, χgel ≈ 6), as it is also evident from a visual inspection of Figures 3 and 6. Finally, the lipid diffusion within the membrane is also affected in the gel state due to the presence of the substrate. Both membrane leaflets are slowed down comparing to the liquid-phase, reaching a lateral diffusion value Dxy−upper ≈ Dxy−lower ≈ (0.65 ± 0.05) × 10−9cm2/s, which is around half of the value found for the free-standing membrane in the gel state Dxy−free ≈ (1.2 ± 0.05) × 10−9cm2/s.
■
DISCUSSION Having thoroughly examined the conditions under which the coarse grained simulations can reproduce the experimental neutron reflectivity data from supported DPPC membranes we may summarize the following points: (a) the presence of a perfectly smooth substrate in relation to the coarse-grained nature of water, leads to layering and depletion effects near the substrate and in the case of the standard water-model in the MARTINI force field gives rise to apparent water freezing even at relatively high temperatures, (b) these effects near the substrate produce sld profiles that show clear deviations from the experimental behavior of the system, (c) a nonsmooth specially prepared substrate resolves these issues especially at temperatures above the liquid-gel lipid transition, (d) at lower temperatures (
