Divergent Diffusion Coefficients in Simulations of Fluids and Lipid

Article pubs.acs.org/JPCB

Divergent Diffusion Coefficients in Simulations of Fluids and Lipid Membranes Martin Vögele and Gerhard Hummer* Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue-Straße 3, 60438 Frankfurt am Main, Germany S Supporting Information *

ABSTRACT: We investigate the dependence of single-particle diffusion coefficients on the size and shape of the simulation box in molecular dynamics simulations of fluids and lipid membranes. We find that the diffusion coefficients of lipids and a carbon nanotube embedded in a lipid membrane diverge with the logarithm of the box width. For a neat Lennard-Jones fluid in flat rectangular boxes, diffusion becomes anisotropic, diverging logarithmically in all three directions with increasing box width. In elongated boxes, the diffusion coefficients normal to the long axis diverge linearly with the height-to-width ratio. For both lipid membranes and neat fluids, this behavior is predicted quantitatively by hydrodynamic theory. Mean-square displacements in the neat fluid exhibit intermediate regimes of anomalous diffusion, with t ln t and t 3/2 components in flat and elongated boxes, respectively. For membranes, the large finite-size effects, and the apparent inability to determine a welldefined lipid diffusion coefficient from simulation, rationalize difficulties in comparing simulation results to each other and to those from experiments. thus a concern.16 Simulations of membranes have been performed with edge lengths of tens of nanometers17 and have recently reached the mesoscale of more than 100 nm,18 giving important insights into the organization and dynamics of lipids and membrane proteins. These length scales are within reach of experimental techniques such as stimulated emission depletion19 and are therefore at the bridge between theory and experiment. A theoretical prediction for the special case of lipid membrane diffusion with increasing edge lengths has recently been presented20 but not yet extensively tested. Preliminary results from simulations of protein and lipid diffusion supported this prediction of strong finite-size effects.20 The aim of this work is to investigate the dependence of diffusion coefficients on box size and geometry, first in a general manner and then with special focus on the application to lipid membrane simulations. We start out by reviewing the hydrodynamic theory of system-size corrections to diffusion coefficients formulated for neat fluids on the basis of the periodic form of the Oseen tensor.12 Applied to noncubic boxes, we recover corrections in accordance with previous findings by Kikugawa et al.21,22 For membrane systems, we adopt the hydrodynamic theory by Camley et al.,20 which is built on the basis of the Saffman−Delbrück model.23 We test

1. INTRODUCTION Diffusion is one of the most fundamental dynamic processes in microscopic and mesoscopic systems, from the translational and rotational diffusion of individual molecules to the Brownian motion of colloidal particles. The accurate treatment of diffusion plays a central role in computational biophysics.1,2 In reflection of the central role of diffusion in the description of dynamic processes,3 the celebrated Stokes−Einstein relation, D = kBT/6πηR, approximates the diffusion coefficient, D, in terms of particle radius R, thermal energy kBT, and fluid viscosity η, thus connecting microscopic dynamics with macroscopic hydrodynamics and thermodynamics. The calculation of diffusion coefficients is a common objective of molecular dynamics (MD) simulations. Most MD simulations are performed under periodic boundary conditions (PBC) to mimic an infinite system without surface.4 In their analysis, it is tacitly assumed that if the simulation boxes are large enough, the artifacts of periodicity are negligible. However, this is not always the case, especially for long-range Coulomb5−8 and hydrodynamic interactions.9−11 For diffusion in liquids,12−15 hydrodynamic self-interactions have been shown to result in large deviations from the infinite-system limit. For a wide range of systems, the correction was found to decay slowly, as 1/L, with increasing length L of cubic simulation boxes.12 In many cases, a cubic simulation box is not the geometry of choice. Membranes in particular are commonly simulated in flat boxes to minimize computational cost. Finite-size effects are © XXXX American Chemical Society

Special Issue: J. Andrew McCammon Festschrift Received: May 20, 2016 Revised: July 1, 2016

A

DOI: 10.1021/acs.jpcb.6b05102 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B the theoretical predictions for flattened and elongated boxes in simulations of a Lennard-Jones (LJ) model of fluid argon. This simple system allows us to study the essence of the finite-size effects without any electrostatic, intramolecular, or surfacerelated interactions influencing the diffusion process. From this general picture, we then turn to the more complex system of a fixed layer of atoms surrounded by a fluid that is effectively confined between the layer and its periodic image. The fixed layer is a simple model for a membrane, and the behavior of the fluid is similar to that in the case of nanoscaled slits, for example, in porous media. As membranes are the most obvious application of asymmetric box shapes, we study the diffusion of lipids and a carbon nanotube (CNT) in a 1-palmitoyl-2-oleoylsn-glycero-3-phosphocholine (POPC) lipid membrane up to near-macroscopic box sizes. Membrane-spanning CNTs have recently been realized in experiment.24 We conclude by discussing open challenges, in particular, concerning the apparent inability to determine a well-defined lipid diffusion coefficient from MD simulations.

DPBC = D0 +

+

k(≠ 0)

r→0

DPBC(L) = D0 −

∑ k(≠ 0)

e−ik·r ⎛ kk ⎞ ⎜I − ⎟ 2 ⎝ k ηV k2 ⎠

(1)

ΔT = (2)

1 ⎛ rr ⎞ ⎜I + ⎟ ⎝ 8πηr r2 ⎠

(5)

kBTξ 6πηL

(6)

1 ⎛ 6α ⎞ ⎜S1 + Sn + Sk − 1/2 I⎟ ⎝ ⎠ 6πη π

(7)

where 2 2 ⎡ 3 erfc(αn) ⎛ 9α ⎟⎞ e−α n ⎤⎥ 3 2 ⎢ ⎜ S1 = ∑ I + 3α n − ⎝ 2 ⎠ π 1/2 ⎥⎦ 4n ⎢ n(≠ 0) ⎣

(3)

Sn =

and for the infinite system T0 =

2

with ξ ≈ 2.837297. For anisotropic simulation boxes, the hydrodynamic formulation in eqs 1−4 suggests that the correction becomes dependent on direction even for an isotropic D0 = D0I. The tensorial correction ΔT depends on both the size and shape of the simulation box. For right-square-prism simulation boxes, Lx = Ly, the anisotropy correction depends only on the ratio of the widths, c = Lz/Lx, and ΔT is diagonal, with elements of the form ΔTii = f ii(c)/6πηLx.22 To evaluate the correction eq 1 in general for the full diffusion tensor, we need a rapidly converging form of the full Oseen tensor under PBC. For this, we use Beenakker’s expression for the Rotne−Prager tensor under PBC,11 leaving out all terms proportional to the squared particle radius, a2. The Oseen tensor correction then becomes

is the hydrodynamic self-interaction, approximated as the difference between the Oseen tensors for PBC TPBC(r) =

n(≠ 0)

erfc(αn) n

The first sum extends over the real-space lattice vectors (n; excluding the origin) defined by the PBC and the second sum over the reciprocal lattice vectors (k; excluding the origin). α is an arbitrary convergence factor (α > 0) chosen to balance the contributions of real-space and Fourier-space sums. For a cubic simulation box of length L, the lattice vectors are n = L(lx, ly, lz) and k = 2πL−1(lx, ly, lz), respectively, where lx, ly, and lz are integers. Evaluation of the lattice sum results in12

where the correction term ΔT = lim[TPBC(r) − T0(r)]

∑

∑

⎤ 4π e−k /4α π ⎥ − k 2V Vα 2 ⎥⎦ 2

2. THEORY 2.1. Hydrodynamic Corrections for the Dependence of the Diffusion Tensor on Box Size and Shape. 2.1.1. General Considerations. Hydrodynamic self-interactions under PBC have been found to account near-quantitatively for the large system-size dependence of self-diffusion coefficients observed in MD simulations of a range of systems.12 In this hydrodynamic formulation, diffusion tensor DPBC (obtained from the slopes of the mean-squared displacements (MSDs) with respect to time in the usual way4) deviates from the ideal infinite-system diffusion tensor, D0, as DPBC = D0 + kBT ΔT

⎡ kBT ⎢ 2α − + 6πη ⎢⎣ π 1/2

Sk =

(4)

2 2 ⎡ 3 erfc(αn) ⎛ 3α ⎟⎞ e−α n ⎤⎥ nn ⎢ − ⎜3α 3n2 − ⎝ 4n 2 ⎠ π 1/2 ⎥⎦ n2 ⎢ n(≠ 0) ⎣

∑

∑ k(≠ 0)

at the origin, r → 0, with kB being the Boltzmann constant and T the absolute temperature. In the Fourier representation of TPBC for the periodic system, eq 3, the sum extends over the vectors, k, of the reciprocal lattice corresponding to the PBC applied, excluding the origin; kk is the matrix with elements kikj of k-vector components; k = |k|; I is the identity matrix; V is the box volume; and η is the bulk shear viscosity of the fluid. The trace of the diffusion tensor defines the scalar diffusion coefficient, DPBC = Tr(DPBC)/3, whose correction is formally analogous to the Coulomb self-interaction of a point charge in a Wigner lattice25 (i.e., a lattice of point charges embedded in a uniform neutralizing background), which can be evaluated by rewriting the formal definition of the trace of TPBC in eq 3 in terms of a rapidly converging Ewald sum12,14

(8)

(9)

2 2 6π e−k /4α ⎛ kk ⎞ k2 k 4 ⎞⎛ ⎜I − ⎟ + + 1 ⎜ 4 ⎟⎝ 2 2 8α ⎠ ⎝ 4α kV k2 ⎠

(10)

are rapidly converging lattice sums in real (n) and Fourier space (k), whose relative contributions are balanced by the coefficient α > 0. 2.1.2. Diffusion Correction for Elongated Simulation Boxes. The dominant asymptotic behavior of the correction term can be calculated analytically in the limit of rectangular elongated simulation boxes of widths Lx = Ly ≪ Lz (i.e., for a matchstick-shaped square cuboid). In the limit of Lz → ∞, a divergence arises in the xx and yy diagonal components of Sk because the reciprocal lattice points k = (0, 0, 2πlzL−1 z ) with lx = ly = 0 approach the origin, where the summand is singular. To estimate the resulting asymptotic corrections to the diagonal components, Dxx = Dyy, of the diffusion tensor, D0, we separate B


Article

The Journal of Physical Chemistry B out the singular part in Sk/6πη and approximate the remainder by an integral −kz2 /4α 2 ⎛

kz2 2

z = 0) Sk(k, zzz = 0) ≈ 2 Sk(k, xx

⎡ 2 2 −k 2 /4α 2 ⎛ k 2 ⎞⎤ 3π 2 ⎢ (6α + k 0 ) e 0 ⎜0, 02 ⎟⎥ 2 ≈ + Γ V δk 2 ⎢⎣ α2 ⎝ 4α ⎠⎦⎥

kz4 ⎞ ⎟ 4

1 (kx = ky = 0) e ⎜1 + = ∑ + Sk , xx 2 6πη 8α ⎠ 4α kz(≠ 0) k z ηV ⎝ ∞ Lz ≈ ∑ 2 2 2 + Jk lz = 1 2π ηLx lz

(18)

For small k0, they diverge as

(11)

z = 0) ≈ const. + Sk(k, zzz = 0) ≈ 2 Sk(k, xx

The integral starts at δk/2, where δk = 2π/Lz is the spacing of k-vectors in the kz-direction. It can be evaluated as ∞

2

ΔTzz ≈ 2ΔTxx ≈ const. +

(12)

which approaches 0 in the limit of δk → 0 (Lz → ∞). The lz −2 sum can be evaluated using ∑∞ = π2/6, which results in n=1n Lz 12ηLx2

(13)

asymptotically for Lz ≫ Lx = Ly. Thus, the hydrodynamic theory predicts that for elongated rectangular boxes the diffusion coefficients, Dxx = Dyy, in the xy plane normal to the long axis diverge linearly with increasing box asymmetry Lz/ Lx. By contrast, the diffusion correction in the z direction of the long axis has no diverging component ΔTzz ≈ const.

Dxx = Dyy ≈ D0 +

Dzz ≈ D0 −

because the zz component of I−kk/k2 in Sk with lx = ly = 0 is zero. 2.1.3. Diffusion Correction for Flat Simulation Boxes. For flat rectangular simulation boxes, Lx = Ly ≫ Lz, both the inplane (Dxx = Dyy) and the out-of-plane diffusion coefficients (Dzz) are predicted to diverge with increasing asymmetry Lx/Lz, as shown in the following. We hold Lz fixed and let Lx = Ly increase to infinity, now replacing the sums in Sk for kz = 0 by integrals in two-dimensional (2D) polar coordinates

∫k

12π 2 V δk 2

∫k

Sk(k, zzz = 0) ≈

∞

∫0

dk 0

2π

dϕ Q (k)(1 − cos2 ϕ)

d k Q (k )

Dzz ≈ D0 +

(21)

6πηLx

(22)

⎤ kBT ⎡ 3 ⎛ Lx ⎞ ⎢ ln⎜ ⎟ − ξxx ⎥ ⎥⎦ 6πηLz ⎢⎣ 2 ⎝ Lz ⎠

⎤ kBT ⎡ ⎛ Lx ⎞ ⎢3 ln⎜ ⎟ − ξzz ⎥ 6πηLz ⎢⎣ ⎝ Lz ⎠ ⎦⎥

(23)

(24)

with ξxx ≈ 2.8897 and ξzz ≈ 2.77939 from numerical fits. At the extreme of a cubic box, L = Lx = Ly = Lz, the asymptotic approximations in eqs 21−24 are in good agreement with those in eq 5. Thus, asymptotic approximations of the full lattice sum are remarkably accurate and cover the entire range of Lx/Lz ratios for right-square-prism boxes. We note that the linear divergence of DPBC in elongated boxes, the logarithmic divergence for flat boxes, and the 1/L convergence in cubic boxes all reflect the close analogy to electrostatics. In one-dimensional (1D), 2D, and threedimensional (3D) systems, Coulomb interactions depend linearly, logarithmically, and inversely on distance, respectively. We note further that in the elongated and flat systems, we can think of the periodic replicates of particles as forming infinite sheets and lines that move “collectively”. 2.2. Connection to the Saffman−Delbrü ck Model. Camley et al.20 have recently modeled the diffusive dynamics in MD simulations of lipid membranes under PBC using an Oseen tensor developed on the basis of the Saffman−Delbrück model23

(15)

(16)

with 2 2 e−k /4α ⎛ k2 k4 ⎞ Q (k ) = + ⎜1 + ⎟ k ⎝ 8α 4 ⎠ 4α 2

⎞ kBT ⎛ π Lz − χxx ⎟ ⎜ 6πηLx ⎝ 2 Lx ⎠

kBTχzz

Dxx = Dyy ≈ D0 +

∞ 0

(20)

with coefficients χxx ≈ 4.3878 and χzz ≈ 2.9252 determined by numerical fits. By contrast, for flat boxes, Lx = Ly ≫ Lz, we have the following asymptotic approximations

(14)

6π z = 0) ≈ Sk(k, xx V δk 2

ln(Lx /Lz) 2πηLz

to the diagonal components of the diffusion tensor in flat rectangular boxes, Lx = Ly ≫ Lz, with different constants for xx and zz. Thus, correction terms diverge logarithmically with the width-to-height ratio irrespective of the precise choice of δk. 2.1.4. Summary of Asymptotic Behavior. The hydrodynamic theory predicts that Dzz for elongated boxes and all diagonal components for flat boxes diverge in the limits of Lz → ∞ and Lx → ∞, respectively. For elongated boxes, Lx = Ly ≪ Lz, we arrive at asymptotic approximations

2

= [(4 + δk 2/4α 2) e−δk /16α − 4]/δk 2ηV

ΔTxx = ΔTyy ≈ const. +

(19)

Substitution of δk = 2π/Lx and V = LxLyLz results in corrections

⎡ e−kz2 /4α 2 ⎛ kz2 kz4 ⎞ 2 1 ⎤⎥ ⎜ ⎟ + + − Jk = dk z ⎢ 1 δkηV δk /2 ⎢⎣ kz2 ⎝ 8α 4 ⎠ kz2 ⎥⎦ 4α 2

∫

α2 6π 2 ln V δk 2 k 02

(17)

where k2 = k2x + k2y and δk = 2π/Lx. As here we did not break out the diverging summands, the magnitude of the approximations depends on the exact location k0 of the integral boundary (e.g., k0 = δk/2 or δk). However, as will be shown in the following, the scaling behavior does not. Both integrals can be evaluated analytically in terms of incomplete Γ functions C


Article

The Journal of Physical Chemistry B TPBC(r) =

⎛ 1 e−ik·r kk ⎞ ⎜I − ⎟ ∑ 2 A k(≠ 0) ηmk + 2ηf k tanh(kH ) ⎝ k2 ⎠ (25)

ηm is the surface viscosity of the membrane, ηf is the viscosity of the surrounding fluid, and H is half of the height of the water layer between two periodic images of the membrane. The sum extends over the nonzero vectors of the reciprocal space in the membrane plane of area A. For a square box of width L = Lx, we have A =L2 and (kx,ky) = 2πL−1(lx,ly) with lx,ly integers. In addition to the changed dimensionality, the Fourier coefficients differ from those in the neat system, eq 3, by the tanh term, which accounts for the influence of the surrounding solvent on the diffusion inside the membrane. The correction to the 2D diffusion coefficient in the membrane plane is calculated in the same way as in eq 1 for 3D diffusion. We have DPBC ≈ D0 + kBT ΔT, where ΔT is half of the trace of the difference between the Oseen tensor, eq 25, and its limit for L,H → ∞. This definition of the correction and its evaluation deviate slightly from that used by Camley et al.,20 where the Oseen tensor was modified by a Gaussian function to account for the finite particle radius, which is scaled by a numerical parameter to fit the Saffman−Delbrück−Hughes− Pailthorpe−White hydrodynamic theory of diffusion in membranes. However, the numerical results for the correction agree qualitatively and near-quantitatively. For the case of L ≫ Lz corresponding to typical simulations of lipid membrane systems, Camley et al.20 predicted an increase in the diffusion coefficient with ln(L). This logarithmic divergence follows again from the approach of the innermost k vectors to the singular point at k = 0 as L → ∞. If one approximates the k sum by an integral, one finds DPBC ≈ D0 +

kBT ln(L*) − ξ(H *) 4πηm (1 + H *)

Figure 1. Theoretical prediction for the difference between the observed and actual diffusion coefficients of lipids in membranes simulated under PBC. Colored lines show DPBC − D0 in units of kBT/ ηm from eq 28 as a function of the reduced box width L* = L/LSD for 16 logarithmically equidistant values of H* = H/LSD from 0.1 to 100 (rainbow colored curves top to bottom). Black dotted lines show the approximation eq 26 for L ≫ H in the range used for the linear fit to obtain the correction term, ξ(H*), in eq 26. The inset shows the resulting ξ(H*) as symbols, and the approximation eq 27 as line. The vertical dashed lines indicate the Saffman−Delbrück length.

2ΔT ≡ −

⎡ ⎛ 1 1 = lim ⎢ 2 ∑ ⎜⎜ 2 σ →∞⎢ L ⎣ k(≠ 0) ⎝ ηmk + 2ηf k tanh(kH ) 2 2 2 ⎤ e−u [π erfi(u) − E i(u 2)] ⎥ 1 − e−k /2σ ⎞⎟ − − ⎟ ⎥⎦ 4πηm ηmk 2 + 2ηf k ⎠

(28)

where u = (√2)ηf/ηmσ, erfi(x) is the imaginary error function and Ei(x) is the exponential integral. To extrapolate from numerical results for finite σ of about 10π/L to infinity, we take advantage of the approximately linear convergence in 1/σ2. As shown in Figure 1, the hydrodynamic correction is well approximated by eq 26. Overall, the hydrodynamic theory predicts that DPBC approaches the limiting value D0 only for boxes that are orders of magnitude larger than the Saffman− Delbrück length in all three dimensions, L ≈ H ≫ LSD. Formally, the above results have been developed for large membrane-spanning objects, not for lipids confined to leaflets. However, for small heights of the simulation box, the correction for lipids, and other objects spanning only one leaflet of the membrane, is similar to the one for transmembrane objects. Adopting a model by Camley et al.26 with two lipid leaflets, the term 1/[ηmk2 + 2ηfk tanh(kH)] in the Fourier expansion of the Oseen tensor in eq 25 is replaced by C(k)/[C2(k) − B2(k)] with C(k) = ηmk2/2 + ηfk coth(2Hk) + b and B(k) = b + ηf k/ sinh(2Hk), where b is the intermonolayer friction coefficient and ηm/2 is the monolayer surface viscosity.20 For small k and H, we find C(k)/[C2(k) − B2(k)] ≈ 1/[ηmk2 + 2 ηfk tanh(kH)], which agrees with the Fourier coefficients in eq 25. As a consequence, the box size correction for lipids and membranespanning molecules is predicted to be nearly identical for small H.

(26)

for L ≫ H, where H* = H/LSD and L* = L/LSD are the system dimensions scaled by the Saffmann−Dellbrück length, LSD = ηm/2ηf. This approximation predicts a logarithmic increase of DPBC asymptotically for large widths L. The H-dependent correction term ξ(H*), obtained numerically from eq 28 below, is well approximated by ⎛ π ⎞ ξ(H *) = ln⎜1 + H *⎟ + e − 1 ⎝ 2 ⎠

∫

1 1 ∑ 2 2 L k(≠ 0) ηmk + 2ηf k tanh(kH ) d2k 1 (2π )2 ηmk 2 + 2ηf k

(27)

as shown in the inset of Figure 1. A plot of the expected absolute error in the diffusion coefficient, DPBC−D0, highlights the problem of determining a lipid diffusion coefficient from MD simulations (Figure 1). To evaluate ΔT, we subtracted the long-wavelength part from the lattice sum and instead approximate it by an integral, which results in D


Article

The Journal of Physical Chemistry B

Figure 2. Effects of simulation box size and shape on single-particle diffusion in a 3D argon fluid. Left: the box height Lz = 4.75 nm is kept constant and the width Lx,y ≡ Lx = Ly is varied (logarithmic scale). Right: the box width Lx = Ly = 6.75 nm is kept constant and the height Lz is varied (linear scale). Red circles show the simulation results for the lateral diffusion coefficient (Dxx + Dyy)/2, averaged over the symmetry-equivalent x and y directions. Blue squares show the vertical component, Dzz. Continuous lines in the respective colors show the theoretical predictions using the viscosity calculated from a set of cubic-box simulations (depicted in the inset) with the respective uncertainty as a shaded region around the lines. The black dotted lines show the asymptotic approximation eqs 21−24. The horizontal gray dotted line shows the limiting diffusion coefficient, D0, obtained by extrapolation to infinite system size from simulations in cubic boxes. Filled symbols show diffusion coefficients calculated from slopes of linear fits to the MSD from 2 to 4.5 ns. Small and large open symbols show the results of fits to earlier time windows of 0.3−0.5 and 0.7−0.9 ns, respectively.

type with a total particle density of 10.97 nm−3 (0.43 in LJ units) and mass density of 728 kg m−3. For the 2D and 2D + 1D systems, we varied the width Lx = Ly, keeping Lz = 4.5 nm fixed in the latter. Size and shape were varied in the same way as for the flat 3D systems. All other simulation parameters of these systems were the same as for the 3D systems (further details are provided in the Supporting Information). 3.2. Lipid Membrane. We simulated POPC lipid membranes in boxes of varying size and shape using the MARTINI coarse-graining scheme.30 In these simulations, a planar lipid bilayer spans the box in the xy plane, with water above and below. MD simulations were performed using Gromacs 4.5.628 with protocols and parameters typical for the MARTINI model.30 Membrane systems in boxes of different sizes and shapes were set up using the script insane.py available at the MARTINI homepage (http://www.cgmartini. nl). Above and below the lipid bilayer, water layers of desired thickness were added. Some of the water beads are of the MARTINI bead type WF, which prevents the bulk from freezing. After 5000 steps of steepest-descent energy minimization, each system was equilibrated for 20 ns by MD, followed by data production for 2 μs. The time step was 0.02 ps. A temperature of 300 K was maintained with a Berendsen thermostat (τT = 1 ps), and a pressure of 1 bar with a Parrinello−Rahman barostat (τp = 12 ps). Note that the effective time scale of coarse-grained models varies for different dynamic observables.31 In the first set of simulations (set A), the initial box height was fixed at Lz = 15 nm, and the width varied from L = Lx = Ly = 10 to 45 nm. In the second set (B), the initial box width was fixed at L = 10 nm, and the height varied from Lz = 10 to 115 nm. In the third set (C), the simulation box was cubic with widths of L = 10−40 nm. In the fourth set (D), we pushed into the near-macroscopic regime, with box widths, L, from 10 nm to 0.4 μm, at initial box heights of Lz = 9 nm. This corresponds to systems with 336 up to 540 800 lipids. In these runs, otherwise strong undulations of the lipid bilayer were

3. MODEL SYSTEMS 3.1. Argon Fluid. 3.1.1. 3D Systems. To test the predictions of the hydrodynamic theory, we simulated a neat LJ fluid in boxes of different sizes and shapes. As a test system, we used argon at T = 240 K and with a mass density of 728 kg m−3, equivalent to a particle density of 10.97 nm−3. For argon OPLS-AA force-field parameters,27 this state point corresponds to a reduced temperature of 2 and a reduced particle density of 0.43 in LJ units. All simulations were performed using Gromacs 4.5.628 in an NVT ensemble using a velocity rescaling thermostat29 with a time constant of 0.4 ps. We used an integration time step of 2 fs and a cut-off radius for the LJ interactions of 1.2 nm. These parameters and the thermostat were chosen to be consistent with typical simulations. However, the absolute values of the diffusion coefficients will reflect in part the choice of thermostat. After an initial equilibration of 1 ns, data were collected in production runs of 100 ns for the different box sizes. We varied the width from Lx = Ly = 3.6 to 180 nm at a fixed height of Lz = 4.5 nm (flat systems), resulting in boxes that contained from 640 to 1 600 000 particles. We also varied the height from Lz = 6.75 to 67.5 nm at a fixed width of Lx = Ly = 6.75 nm (elongated systems), resulting in boxes that contained from 3375 to 33 750 particles. For reference, we performed simulations in cubic boxes of length L = 3.15−18 nm, containing from 343 to 64 000 particles (further details are provided in the Supporting Information). 3.1.2. 2D and 2D + 1D Systems. To explore the effect of dimensional restriction, we also simulated 2D Ar systems and 2D + 1D systems, in which the motion of a subset of otherwise identical Ar particles was restricted to the xy plane. In the purely 2D systems, all Ar atoms are confined to a plane. The area density of the 2D-confined atoms was 4.94 nm−2 (0.57 in LJ units). We varied the edge length from Lx = Ly = 3.6 to 180 nm, resulting in boxes with 64−160 000 particles. In the 2D + 1D system, we embedded the 2D system of confined atoms in a box of 10 nm height, filled with unrestrained atoms of the same E


Article


Figure 3. Derivatives of the MSD with respect to time for flat 3D argon systems. Curves for increasing Lx/Lz are colored from purple to red. Top: MSD′(t) vs t. Bottom: MSD′(t) vs ln t. Left: (MSD′xx + MSD′yy)/2. Right: MSD′zz. For the largest systems, a logarithmic fit to the first 1.5 ns is shown as a black dashed line. The gray background shows the region used for the calculation of the diffusion coefficients.

Figure 4. Derivatives of the MSD with respect to time for elongated 3D argon systems. Curves for increasing Lz/Lx are colored from purple to red. Top: MSD′(t) vs t. Bottom: MSD′(t) vs t1/2. Left: (MSD′xx + MSD′yy)/2. Right: MSD′zz. For the largest systems, a square-root fit to the first 1.5 ns is shown as a black dashed line. The gray background shows the region used for the calculation of the diffusion coefficients.

suppressed by a weak harmonic restraint acting on the z coordinate of the center of mass of a quarter of the lipids. This approach has already been tested and used for mimicking the effect of a cytoskeleton network.17 Further details are provided in the Supporting Information.

As an example for a membrane-spanning object, we used a CNT whose rims were functionalized with polar groups. For each setup in sets A and B, otherwise identical systems were simulated with and without a CNT embedded in the membrane. The topology of CNT is adapted from an existing F


Article


Figure 5. Derivatives of the MSD with respect to time for 2D argon systems. Curves for increasing L are colored from purple to red. Top: linear time scale. Bottom: logarithmic time scale. Left: purely 2D system. Right: 2D system embedded in a 3D fluid. For the largest systems, a logarithmic fit to the first 1.5 ns is shown as a black dashed line. The gray background shows the region used for the calculation of the diffusion coefficients.

model for capped nanotubes.32 We removed the caps and instead added a ring of weakly polar beads (SNda) at the two ends, which keep the CNTs upright in the membrane. In all simulations, the same nanotube configuration was used: 12 rings of 8 CG beads. This results in a length of 4.48 nm and a diameter of 1.23 nm (distances from the positions of the CG beads). See the Supporting Information for further details. 3.3. Diffusion Analysis. The diagonal elements of the diffusion tensor were determined from least-squares fits of linear functions a + 2Dt to the MSD(t) in a time window of t0 to t1 in the i direction (i ∈ {x,y,z}).

increase is possibly due to thermostatting. Remarkably, the predictions of the hydrodynamic theory appear to hold over the entire range of box-shape asymmetries, both for the flat and elongated box shapes. Apparent deviations for highly distorted box shapes disappear if we move the window of the MSD fit further out in time (small open symbols for early time windows and solid symbols for late windows). An intermediate regime of anomalous self-diffusion is the origin of this dependence of Dii on time in highly distorted boxes. For flat systems, the derivatives of the MSDs with respect to time exhibit a logarithmic phase, MSD′(t) ∝ ln(t), before the plateau corresponding to normal diffusion is reached (Figure 3). Thus, in this regime, the MSD contains a t ln(t) component expected for ideal 2D systems.35 For elongated systems, we find MSD′(t) ∝ t1/2 for the dynamics in the xy plane at intermediate times (Figure 4). By contrast, the derivative of the MSD in the z direction is constant for elongated boxes. Importantly, the anomalous regime invariably transitions to normal diffusion after some time. However, this time increases as the box geometry becomes more extreme. Reliable estimates of the diffusion coefficients in distorted boxes thus require extensive simulations, both to establish normal diffusion and to suppress statistical noise. 4.1.2. 2D and 2D + 1D Systems. To explore the effects of extreme dimensional restriction, we also simulated purely 2D systems, and systems in which some particles were restricted to 2D motion, but surrounded by a fluid without restrictions. In the purely 2D systems, the MSD exhibits a pronounced t ln(t) dependence (Figure 5 left). However, this hallmark of 2D diffusion35 disappears after some time. This time to establish normal diffusion increases with the width, Lx, of the 2D square system. In the 2D + 1D system, the t ln(t) phase also develops, albeit in a less-pronounced manner (Figure 5 right). For the related case of a colloid confined in 2D and surrounded by a 3D

N

MSDii (t ) ≡ ⟨(ri(t + τ ) − ri(τ ))2 ⟩ ≈

∫0

T−t

(riα(t + τ ) − riα(τ ))2 dτ

1 ∑ 1 N α=1 T − t (29)

where the sum extends over the N identical particles. The MSD was calculated using an efficient Fourier-based algorithm.33 Statistical errors were estimated by block averaging, dividing the simulation runs into 20 blocks that were assumed to be statistically independent. Because of the lateral symmetry of the systems, we average (Dxx + Dyy)/2.

4. RESULTS AND DISCUSSION 4.1. Argon Fluid. 4.1.1. 3D Systems. Figure 2 shows the self-diffusion coefficients of the argon fluid simulated in boxes of different sizes and shapes together with the prediction from hydrodynamic theory. The effective viscosity coefficient η entering the correction formulas for Dxx, Dyy, and Dzz was obtained by a least-squares fit of eq 5 to a series of simulations of cubic systems (inset in Figure 2). The value of η = 7.64 × 10−5 Pa s obtained in this way agrees well with literature values interpolated to our state point (7.5 × 10−5 Pa s).34 The small G


Article

The Journal of Physical Chemistry B fluid, a divergence of the collective diffusion has been reported.36,37 We note that the free fluid in the 2D + 1D simulations is effectively a thin fluid film confined between two plates. We show in the following paragraph that in this case the same logarithmic divergence of the in-plane diffusion coefficient of unrestrained particles is expected as for the flat boxes without an embedded layer of fixed particles. Thus, the effects seen in the MSDs of the 3D system and attributed to dimensional restriction are indeed present in an actual 2D system and a 2D + 1D system. 4.1.3. Confined Films. The hydrodynamic theory predicts the same logarithmic increase for flat fluid films confined between two planes as it does for the flat periodic systems. The confinement in the vertical direction only affects the summation over kz. The arguments for the case kz = 0 that lead to eq 15 are the same regardless of the summation in the z direction and the boundary conditions (stick, slip, etc.) applied. Thus, the theory predicts a logarithmic increase for thin films as well. Only the numerical constant ξ is expected to differ. The qualitative conformity of the two types of systems can be seen in the simulations (Figure 6). The freely moving argon

the hydrodynamic theory and demonstrate its applicability also to dynamics in confined fluids. 4.2. Lipid Membrane. The MSD of lipids diffusing in the membrane shows an initial rapid increase followed by a transition to normal diffusion at longer times (Figure 7). After about 3 ns, the MSDs are linear for all box sizes. The lipid diffusion coefficients obtained from linear fits to the MSD in a time window from 4 to 9 ns are shown in Figure 8. Also shown are the much slower diffusion coefficients of CNTs embedded in the membranes (right-hand axis). In all cases, the diffusion coefficients increase proportionally to the logarithm of the box width Lx = Ly. Up to the largest box sizes, the diffusion coefficients of lipids and CNTs exhibit the same logarithmic dependence on the box size. The identical slopes in Figure 8 may seem surprising because lipids are contained in leaflets and CNTs span the entire membrane. Moreover, the Saffman−Delbrück model assumes the diffusing objects to be much larger than the constituent molecules. However, the observed independence of the PBC correction from the size of the object and its detailed location within the membrane is consistent with expectations from a two-leaflet hydrodynamic theory sketched above using the formulation by Camley et al.20,26 Note that transient lipid− CNT interactions explain the slightly lower diffusion coefficient of the lipids in the smallest systems with a CNT (Figure 8, blue triangles) compared to that in the systems without a CNT (gray circles). The dominant hydrodynamic corrections to the apparent diffusion are absolute, not relative. As a consequence, slowly diffusing molecules suffer from finite-size effects in a morepronounced manner. Indeed, the diffusion coefficient of the membrane-embedded CNT almost doubles over the range of box sizes considered (Figure 8, green squares). Even more problematic is the fact that this logarithmic growth does not settle for boxes as large as 100 nm (for CNTs) and 400 nm (for lipids). The logarithmic increase holds true for the whole range over which lipid bilayer simulations have been performed. These results make it clear that from these data (or, for that matter, any other simulation of a membrane system using a similar setup and resolution level) one cannot safely extract a reliable diffusion coefficient that could be compared to that from experiment. We also examined the dependence of lipid diffusion on box height at a constant width of 10 nm. As shown in Figure 9, increasing the height, Lz, to 115 nm (and thus increasing H) has no effect. This behavior is expected from the hydrodynamic model. Note that eq 26 is not applicable here because L is small. Instead, we can approximate tanh(kH) ≈ 1 in eq 25. The denominator and thus the correction become independent of H, consistent with the data in Figure 9. For the narrow box, diffusion inside the membrane is only negligibly influenced by the height of the simulation box. Indeed, the hydrodynamic theory (Figure 1) predicts a significant influence of the height only for widths far beyond the Saffman−Delbrück length. Typical values for MARTINI lipid simulations are around 10 nm (e.g., LSD = 8.6 nm for DPPC lipids20,38). A noticeable deviation from the logarithmic divergence is thus only expected for L,H > 50−100 nm. Our results confirm the predicted logarithmic divergence of DPBC with the box width in membrane simulations and the independence of H for small L. However, to test the hydrodynamic theory of Camley et al.20 and eqs 25−28 more thoroughly, simulations of systems large in all three dimensions

Figure 6. Comparison of the lateral diffusion coefficients (Dxx + Dyy)/ 2 of unrestrained particles in the flat systems, in which argon particles were allowed to diffuse freely in all directions (red), and in the 2D + 1D systems containing an added 2D layer of particles with fixed z coordinates that confines the intervening fluid (blue). The red line is the prediction for a periodic system. The blue line is shifted by addition of a constant, to match the confined 2D + 1D dynamics. Red dotted and blue dashed lines are the asymptotic estimates. The value D∞ of an infinite cubic system is depicted as horizontal gray dotted line.

atoms of the 2D + 1D system can be interpreted as a thin film of fluid confined between two surfaces. Figure 6 compares the lateral diffusion coefficient of this system to the one of the freely moving fluid in a flat box. Despite small differences in the density, because a fraction of 10% of the particles is fixed in a plane of the 2D + 1D system, the logarithmic increase is very similar for both systems. The absolute value of the diffusion coefficient is shifted to slightly higher values for the confined film. From a fit of the simulation results to the functional form predicted by the theory (eq 23), the numerical constant, ξ, can be estimated. Excluding the first four data points (where the dependence is not yet logarithmic and the approximation is not applicable) and the last two (for which the MSD has not quite reached the regime of linear t dependence), a value of ξ ≈ 2.19 is obtained. These simulation results confirm the predictions of H


Article


Figure 7. Diffusion in lipid membrane. Left: MSD curves of diffusing POPC lipids in lipid membrane simulations for different box widths between 10 and 400 nm at a constant box height of 9 nm. Right: derivatives of the respective MSD curves. The gray background shows the region used for the calculation of the diffusion coefficients.

Figure 9. Effects of simulation box height Lz on the lateral diffusion coefficient in lipid membrane simulations with fixed L = Lx = Ly = 10 nm. The results of two sets of POPC lipid membrane simulations are shown together with linear fits (lines). Blue triangles: Dlip from unrestrained membrane simulations including a CNT. Green squares: DCNT of the CNT. Red circles: Dlip from unrestrained membrane simulations without CNT.

Figure 8. Effects of simulation box size on the lateral diffusion coefficient in lipid membrane simulations. The results of four sets of POPC lipid membrane simulations are shown as a function of the box width, L = Lx = Ly, together with fits to a logarithmic function. Blue triangles: lipid diffusion coefficient Dlip from simulations of an unrestrained membrane in boxes of height Lz = 15 nm including a CNT. Green squares: diffusion coefficient DCNT of the CNT as a model transmembrane object (right-hand scale). Gray circles: Dlip from simulations with unrestrained membranes in cubic simulation boxes. Black dots: Dlip for a membrane in boxes of constant width Lx and varying heights Lz. Red stars: Dlip for a membrane whose undulations in z were gently restrained, with Lz = 9 nm. The inset shows side views of the lipid head groups in a freely undulating membrane (bottom, blue) and a restrained membrane (top, red).

preclude the possibility of determining meaningful diffusion coefficients from such simulations. However, not all is lost. Simple hydrodynamic theories account quantitatively for the observed behavior, akin to the success of continuum electrostatic corrections for Coulomb periodicity artifacts.5−8 For the neat fluid, hydrodynamic selfinteractions arise as a result of imposing periodicity on the solutions of hydrodynamic transport theory.12,13 The magnitude of these self-interactions can be estimated from the difference in the Oseen tensors of the periodic and free system. In our simulations of fluid argon, the resulting correction works well even for the most extreme box asymmetries. One caveat is that for the most strongly distorted systems one has to wait a considerable time until normal diffusion is established. Before, the MSD as a function of time exhibits the hallmarks of dimensional restriction, containing significant t3/2 and t ln t components for elongated and flat boxes, respectively. For lipid membranes, the situation is more problematic. Consistent with the predictions of the hydrodynamic theory based on the Saffman−Delbrück model,20 the logarithmic divergence persists even for near-macroscopic systems of 0.4 μm width. As a result, one cannot immediately extract the proper diffusion coefficients of a lipid or of a membrane-

would have to be performed, L ≈ H ≫ LSD, which is beyond our current possibilities.

5. CONCLUSIONS In MD simulations of a simple fluid and of lipid membranes, we found that single-particle diffusion coefficients do not converge to the proper infinite-system limit even for boxes approaching the infinite volume limit. This pathological behavior arises when the simulation box is not scaled uniformly. Such anisotropic scaling is common in simulations of membrane systems, molecular rods, or filamentous structures. In highly asymmetric boxes, the apparent diffusion coefficient depends strongly on the direction and size. For elongated boxes, diffusion normal to the long axis speeds up linearly with increasing box height, without apparent bound. For flat boxes, the diffusion in all directions increases logarithmically with the ratio of width to height. In essence, these divergences appear to I



■

embedded molecule. Assuming that hydrodynamic theory holds, one would have to simulate systems that are large not only in the xy plane but also in the z direction, 1−2 orders of magnitude larger than the Saffman−Delbrück length LSD. By varying L and H, one could then use eq 26 to fit D0 and ηm, with ηf, H, and L known. Currently, however, such simulations are prohibitively expensive, even for the coarse-grained MARTINI model. As an alternative, one can determine ηm from nonequilibrium simulations38 and use eq 28 or the approximation eq 26 to estimate the correction term. Our analysis suggests that a logarithmic divergence arises also in simulations of confined fluids,39 beyond bulk fluids and membrane systems studied here. Diffusion in fluid-filled slit pores, for instance, serves as a model of transport in porous media. Considering the practical relevance of this and related problems, it will be interesting to explore the impact of PBCs in simulations also of dimensionally confined systems. The divergence of the diffusion coefficient with increasing box width is a serious problem when it comes to comparing dynamic properties obtained from different simulations and from the experiment. Unlike 3D diffusion in cubic boxes, one cannot easily extrapolate to infinite system size. A particularly insidious aspect is that the logarithmic dependence appears to be small when probed over a typical range of box lengths, say, from 10 to 20 nm for a lipid membrane, and possibly smaller than the statistical errors. However, it is unbounded. The correction term (or alternatively the size of the system at which the infinite-size value is reproduced) can be computed theoretically only if the exact values for the viscosity of the fluid and of the membrane are known. Thus, obtaining an accurate estimate of the correction is not straightforward for complex membranes consisting of many different components. The theory suggests that by increasing the box height, H, together with the width, L, to very large values, it would be possible to get close to the infinite-system value. In such a calculation, enormous computational effort would be “wasted” on calculating solvent−solvent interactions, to model the hydrodynamics of the bulk. Clearly, this is not desirable and asks for novel theoretical and algorithmic formulations. In the meantime, great care should be exercised when comparing diffusion coefficients of lipid membranes from experiment and simulation.

■

REFERENCES

(1) Ermak, D. L.; McCammon, J. A. Brownian Dynamics with Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352−1360. (2) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; et al. Electrostatics and Diffusion of Molecules in Solution: Simulations with the University of Houston Brownian Dynamics Program. Comput. Phys. Commun. 1995, 91, 57−95. (3) Zwanzig, R. Nonequilibrium Statistical Mechanics; Oxford University Press: New York, 2001. (4) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (5) Neumann, M. Computer Simulation and the Dielectric Constant at Finite Wavelength. Mol. Phys. 1986, 57, 97−121. (6) Hummer, G.; Pratt, L. R.; García, A. E. Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100, 1206−1215. (7) Hünenberger, P. H.; McCammon, J. A. Ewald Artifacts in Computer Simulations of Ionic Solvation and Ion-Ion Interaction: A Continuum Electrostatics Study. J. Chem. Phys. 1999, 110, 1856−1872. (8) Lin, Y. L.; Aleksandrov, A.; Simonson, T.; Roux, B. An Overview of Electrostatic Free Energy Computations for Solutions and Proteins. J. Chem. Theory Comput. 2014, 10, 2690−2709. (9) Oseen, C. W. Neuere Methoden und Ergebnisse in der Hydrodynamik; Mathematik und ihre Anwendungen in Monographien und Lehrbüchern; Akademische Verlagsgesellschaft: Leipzig, 1927. (10) Yamakawa, H. Modern Theory of Polymer Solutions; Harper & Row: New York, 1971. (11) Beenakker, C. W. J. Ewald Sum of the Rotne-Prager Tensor. J. Chem. Phys. 1986, 85, 1581. (12) Yeh, I. C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108, 15873−15879. (13) Dünweg, B.; Kremer, K. Microscopic Verification of Dynamic Scaling in Dilute Polymer Solutions: A Molecular-Dynamics Simulation. Phys. Rev. Lett. 1991, 66, 2996−2999. (14) Dünweg, B.; Kremer, K. Molecular Dynamics Simulation of a Polymer Chain in Solution. J. Chem. Phys. 1993, 99, 6983−6997. (15) Fushiki, M. System Size Dependence of the Diffusion Coefficient in a Simple Liquid. Phys. Rev. E 2003, 68, No. 021203. (16) Klauda, J. B.; Brooks, B. R.; Pastor, R. W. Dynamical Motions of Lipids and a Finite Size Effect in Simulations of Bilayers. J. Chem. Phys. 2006, 125, No. 144710. (17) Ingólfsson, H. I.; Melo, M. N.; Eerden, F. J. V.; Arnarez, C.; López, C. A.; Wassenaar, T. A.; Periole, X.; Vries, A. H. D.; Tieleman, D. P.; Marrink, S. J. Lipid Organization of the Plasma Membrane Lipid Organization of the Plasma Membrane. J. Am. Chem. Soc. 2014, 136, 14554−14559. (18) Koldsø, H.; Sansom, M. S. P. Organization and Dynamics of Receptor Proteins in a Plasma Membrane. J. Am. Chem. Soc. 2015, 137, 14694−14704. (19) Eggeling, C.; Ringemann, C.; Medda, R.; Schwarzmann, G.; Sandhoff, K.; Polyakova, S.; Belov, V. N.; Hein, B.; von Middendorff, C.; Schönle, A.; Hell, S. W. Direct Observation of the Nanoscale Dynamics of Membrane Lipids in a Living Cell. Nature 2009, 457, 1159−1162. (20) Camley, B. A.; Lerner, M. G.; Pastor, R. W.; Brown, F. L. H. Strong Influence of Periodic Boundary Conditions on Lateral Diffusion in Lipid Bilayer Membranes. J. Chem. Phys. 2015, 143, No. 243113. (21) Kikugawa, G.; Ando, S.; Suzuki, J.; Naruke, Y.; Nakano, T.; Ohara, T. Effect of the Computational Domain Size and Shape on the Self-Diffusion Coefficient in a Lennard-Jones Liquid. J. Chem. Phys. 2015, 142, 024503. (22) Kikugawa, G.; Nakano, T.; Ohara, T. Hydrodynamic Consideration of the Finite Size Effect on the Self-Diffusion Coefficient in a Periodic Rectangular Parallelepiped System. J. Chem. Phys. 2015, 143, No. 024507.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b05102. System parameters and results of all simulation runs (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +49 69 63032501. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Drs. Richard Pastor, Jürgen Köfinger, Lukas Stelzl, and Max Linke for helpful discussions. This work was supported by the Max Planck Society. J


Article

The Journal of Physical Chemistry B (23) Saffman, P. G.; Delbrück, M. Brownian Motion in Biological Membranes. Proc. Natl. Acad. Sci. U.S.A. 1975, 72, 3111−3113. (24) Geng, J.; Kim, K.; Zhang, J.; Escalada, A.; Tunuguntla, R.; Comolli, L. R.; Allen, F. I.; Shnyrova, A. V.; Cho, K. R.; Munoz, D.; et al. Stochastic Transport through Carbon Nanotubes in Lipid Bilayers and Live Cell Membranes. Nature 2014, 514, 612−615. (25) Nijboer, B. R. A.; Ruijgrok, T. W. On the Energy per Particle in Three- and Two-Dimensional Wigner Lattices. J. Stat. Phys. 1988, 53, 361−382. (26) Camley, B. A.; Brown, F. L. H. Diffusion of Complex Objects Embedded in Free and Supported Lipid Bilayer Membranes: Role of Shape Anisotropy and Leaflet Structure. Soft Matter 2013, 9, 4767− 4779. (27) Jorgensen, W. L.; Tirado-Rives, J. The OPLS [Optimized Potentials for Liquid Simulations] Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (28) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; Van Der Spoel, D.; et al. GROMACS 4.5: a High- Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (29) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, No. 014101. (30) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; Vries, A. H. D. The MARTINI Force Field : Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (31) Marrink, S. J.; Tieleman, D. P. Perspective on the Martini Model. Chem. Soc. Rev. 2013, 42, 6801−6822. (32) Baoukina, S.; Monticelli, L.; Tieleman, D. P. Interaction of Pristine and Functionalized Carbon Nanotubes with Lipid Membranes. J. Phys. Chem. B 2013, 117, 12113−12123. (33) Calandrini, V.; Pellegrini, E.; Calligari, P.; Hinsen, K.; Kneller, G. R. nMoldyn-Interfacing Spectroscopic Experiments, Molecular Dynamics Simulations and Models for Time Correlation Functions. Collect. SFN 2011, 12, 201−232. (34) Rowley, R. L.; Painter, M. M. Diffusion and Viscosity Equations of State for a Lennard-Jones Fluid Obtained from Molecular Dynamics Simulations. Int. J. Thermophys. 1997, 18, 1109−1121. (35) Keyes, T.; Oppenheim, I. Bilinear Hydrodynamics and the Stokes-Einstein Law. Phys. Rev. A 1973, 8, 937−949. (36) Bleibel, J.; Domínguez, A.; Oettel, M. Hydrodynamic Interactions Induce Anomalous Diffusion under Partial confinement. J. Phys.: Condens. Matter 2015, 27, No. 194113. (37) Lin, B.; Cui, B.; Xu, X.; Zangi, R.; Diamant, H.; Rice, S. A. Divergence of the Long-Wavelength Collective Diffusion Coefficient in Quasi-One- and Quasi-Two-Dimensional Colloidal Suspensions. Phys. Rev. E 2014, 89, No. 022303. (38) den Otter, W. K.; Shkulipa, S. A. Intermonolayer Friction and Surface Shear Viscosity of Lipid Bilayer Membranes. Biophys. J. 2007, 93, 423−433. (39) Mittal, J.; Truskett, T. M.; Errington, J. R.; Hummer, G. Layering and Position-Dependent Diffusive Dynamics of Confined Fluids. Phys. Rev. Lett. 2008, 100, No. 145901.

K


Divergent Diffusion Coefficients in Simulations of Fluids and Lipid

Recommend Documents