Brownian Translational Dynamics on a Flexible Surface: Nuclear Spin Relaxation of Fluid Membrane Phases.

A general model for nuclear magnetic resonance (NMR) relaxation studies of fluid bilayer systems is introduced, combining a mesoscopic Brownian dynamics description of the bilayer with atomistic molecular dynamics (MD) simulations. An example is given for dipalmitoylphosphatidylcholine in 2H2O solvent and compared with the experiment. Experimental agreement is within a factor of 2 in the water relaxation rates, based on a postulated model with fixed parameters, which are largely available from the MD simulation. Relaxation rates are particularly sensitive to the translational diffusion of water perturbed by the interface dynamics and structure. Simulation results suggest that a notable deviation in the relaxation rates may follow from the commonly used small-angle approximation of bilayer undulation. The method has the potential to overcome the temporal and spatial limitations in computing NMR relaxation with atomistic MD, as well as the shortcomings of continuum models enabling a consistent description of experiments performed on a solvent lipid and added spin probes. This work opens for possibilities to understand relaxation processes involving systems such as micelles, multilamellar vesicles, red blood cells, and so forth at biologically relevant timescales in great detail.

Introduction

Lyotropic phases display complex dynamics involving, in addition to the intramolecular degrees of freedom (e.g., molecular rotation and vibrations), thermally excited interface modes that feature collective molecular motion of lipids and solvents.[1] Understanding these processes provides insight into how fundamental electrostatic and van der Waals interactions lead to colloidal interactions, such as the repulsive entropic force between bilayers.[1] This is of fundamental importance for the understanding of biological systems.[2] Likewise, in refining practical studies, for example, antibiotic and model membranes,[3] a microscopic understanding of the system is required.

In this context, nuclear magnetic resonance (NMR) is a supreme method that provides information via its spectral line shape, relaxation, and diffusion in a noninvasive manner. Particularly, sensitivity to slow timescale phenomena is seen in low-field magnetization decay (T1–1) and transverse (spin–spin) relaxation (T2–1) rates.[4−6] Experimental NMR data inform, for instance, on the molecular order via the quadrupole splitting in the spectra, as well as on molecular dynamics (MD) via the relaxation rates and line width.[4,7] However, because experimental NMR data, on their own, only provide qualitative information, efforts have been made during the past four decades to also mathematically model bilayers and to understand NMR relaxation experiments performed on them. The continuum theory of Helfrich[8] has inspired the description of the dynamics included in NMR relaxation models. One of the first assumptions therein is that the interfacial small-amplitude undulation is the only mechanism that modulates the director system of the spin moiety and, thus, is the cause of T–1 (i = 1 and 2) relaxation. This model provides a dependence of T1–1 proportional to ω0–1 (where ω0 is the resonance frequency).[9,10] Further development has addressed the overall spherical geometry, aiming at an NMR relaxation model for micelles.[11]

The continuum theory has been taken further by considering stacked membranes with translationally diffusing spin moiety, which provides an additional mechanism in modeling the relaxation rates.[12] Because of the approximations required by the continuum theory, the spectral density function is formulated to second[12] and fourth order[13] in the director fluctuations of the spin moiety.

An attractive complement to the continuum theory is provided by atomistic MD simulations of bilayer systems, with an explicit water solvent. As an example, the electron paramagnetic resonance (EPR) line shape was, for the first time, computed directly in the time domain without any perturbative approximations, based on a 100 ns long MD trajectory.[14] This provided a model for the spin probe dynamics of the dipalmitoylphosphatidylcholine/water (DPPC) carbon chain. 2H2O relaxation was also investigated based on the same simulation.[15] Common to these works is that static properties, such as the order parameters of the lipid chain and water, can be reproduced in a reasonable agreement with the experiment. This provides support for detailed models of the water structure in the vicinity of the bilayer, as obtained from the MD simulation. However, other important conclusions that ultimately motivate the present work included the fact that the dynamics of the EPR spin probe seem to be too fast by at least a factor of 10, and that of water by a factor of 100, to reproduce the results of line shape and relaxation experiments, respectively. A smaller deviation is found when the high-field magnetization of a 2H-labeled hydrocarbon chain is modeled with MD.[16] The latter observation can be ascribed to the dominant dependence on the fast (subnanosecond) MD that determines the magnetization relaxation at high field, and the fact that MD simulations describe the relevant motional modes well.

At an intermediate model size between the realms of atomistic and continuum theory, Brownian dynamics (BD) simulation has been developed to model diffusion on curved interfaces, with one application area in the study of nonlamellar phases.[17] Within the usual assumption that translational motion on the curved surface is the dominating dynamical model for quadrupole relaxation, this simulation technique has provided valuable additional insight on the curvature of the bicontinuous cubic phase. Simulations have complemented what can be concluded from X-ray studies, as the latter provide only limited information of powder samples.[18]

The continuum modeling approach has limitations that hamper further progress and which may be summarized as follows: (i) studying the physics of close membrane–membrane interactions is not possible because the molecular details (for instance, the steric hindrance) of the lipid head group and the water structure are crucial;[1] (ii) small-angle approximations of the director field are required in the case where the exact range of the validity of the computed spectral density is not known; (iii) the processes of (a) the translational diffusion of the spin moiety and (b) interface dynamics occur on the same timescale, which generally rules out models where the two are statistically separated;[12] and (iv) finally, the presence of exceptionally slow dynamics typically requires a generalization to a nonperturbative relaxation calculation, and the notion that realistic modeling of this motion demands the application of a trajectory-based simulation method[14,19,20] rather than a continuum approach. Although continuum modeling has been performed for liposome systems,[11,21] the interplay of curved, undulating membrane, correlated translational motion along the interface and its effect on relaxation observables, still presents a challenge. These topics can, in part, be addressed with the methodology presented in the present paper and, in part, through forthcoming developments.

Addressing NMR relaxation modeling purely with atomistic MD simulation is hampered by the typical system dimension accessible, circa 4–20 nm and the time span of tens to hundreds of nanoseconds, ranges that are too small to capture the key dynamic modes of relevance to low-field T1 or T2 observables. Furthermore, it is noted that in the majority of MD simulations, a single bilayer is modeled with periodic images, which, therefore, provides an incomplete account of the entropy of the system. Expanding the simulated system size, for example, to consider a stack of bilayers, vastly increases the simulation time. This suggests that MD simulation of bilayers on its own is far too costly to fulfill the task of a generally applicable magnetic resonance relaxation model. It is noted that the so-called coarse-grained MD simulations[22] constitute an alternative route when atomistic MD is insufficient. However, even such methods drastically increase the computational cost relative to the BD simulation approach presented here and, additionally, are expected to be hampered by statistical noise in the computed observables.

In this work, we report a proof-of-principle study of the 2H2O solvent relaxation and line shape, in which the shortcomings of atomistic MD[15] are largely overcome by combining it with a mesoscopic BD model. A spatially discretized bilayer sheet is used to simulate a time-dependent interface,[23] which allows for interface-correlated translational dynamics.[17] With a timescale separation in the relaxation model, all of the fast molecular processes are accounted for in atomistic detail. This provides a close connection to both MD simulations and continuum models, avoiding, however, the approximations of the latter.[12] Furthermore, with the present trajectory-based method, generalization to the nonperturbative Liouville formalism is readily achievable, as has been illustrated with a water relaxation model with a trajectory length in the range of milliseconds.[20]

Atomistic MD plays an important role in drastically reducing the number of adjustable parameters in relaxation models. In future work, quantum-chemically (QC) computed, time-dependent interaction tensors appropriate to individual MD snapshots open the possibility for achieving increased realism.[24,25]

The simulation of large sheets of Brownian interfaces has been of interest[23,26] also concerning translational diffusion on the membrane interface.[27,28] However, in these works, the focus has been on the minor impact on diffusion observable at typical experimental conditions. In the present work, the NMR relaxation properties are addressed and it is shown that the type of simulation introduced here can have a real experimental impact even in the case of a stiff membrane.

The paper is organized as follows. After a review and an adaptation of the existing perturbative water relaxation model[29,30] in section , section explains the simulation of the undulating surface implemented in the Fourier space (FS). Brownian motion on a time-dependent, undulating surface is discussed in section . This is followed by an elaboration of the results in section . First, a lipid or other spin-labeled moiety is modeled in a general way independent of the solvent, followed by a discussion of the relaxation of heavy water. Finally, the conclusions in section summarize the new insights gained on water dynamics.

NMR Relaxation Model

We will discuss the line shape, as well as the T1–1 and T2–1 parameters, that is, the magnetization and transverse relaxation rates, respectively. These observables are compactly available from the Bloch–Wangsness–Redfield (BWR) relaxation theory.[7,31] The key equations for an isotropic powder sample of a bilayer system are given, following the existing theory.[29,32] The system is introduced in Figure , where panel (a) displays a stacked bilayer system with the coordinate system fixed to the director (D) at an angle βLD relative to the laboratory system (L). The interface dynamics enters the Euler angles ΩDn relating D with the normal of the bilayer patch n, depicted in panel (b). After the residence time τm, the water molecule sitting at site d moves to the bulk region or to a new water site d′ with a large number of different sites present, for instance, depending on where on the lipid the water approaches. The relatively fast processes include change of sites (d → d′) that provides the time dependence of Ωnd as well as the local motion of the water molecular axis system (M), described by the Euler angles ΩdM. Finally, panel (c) illustrates models of the time-dependent n(t), incorporating both translational diffusion and interface undulations.

Figure 1

Panel (a): Model system of stacked bilayers separated by the distance(s). The time-dependent Euler angles ΩDn relate the director (D) with the local bilayer normal frame (n). Powder samples are considered, which implies averaging over the laboratory-director orientations (βLD). Panel (b): Angles Ωnd relate the local bilayer normal to the water site director system (d). Local dynamics is described by ΩdM specifying the water molecular axis system. Panel (c): Modes for local bilayer normal frame (n) dynamics.

The equation of motion for the deuterium spin (nuclear spin quantum number I = 1), expressed with the reduced density matrix ρ(t), is given bywhere L0 and R are the coherent and relaxation matrices, respectively. The Zeeman Hamiltonian contributing to L0 has the three energy states labeled H = 1̅, 0, 1. The quadrupole Hamiltonian[7]HQ(t) provides the relaxation mechanism in R. The relaxation matrix is nine-dimensional (including coupling between spin populations and coherencies). However, at high magnetic fields,[7] the 2 × 2 subset is sufficient. The matrix components arewhere ω0 is the Zeeman resonance frequency and ωQLD is the coherent part of the quadrupole interaction. By solving the Fourier–Laplace form of eq , a spectral line shape with angular (βLD) dependence is obtained as[33]from which the line shapeis obtained by numerical integration over a spherical distribution of angles.

Model Dynamics

For further characterization of the experimental observables, a physical model of the water order and dynamics is summarized in this section. The focus is on the dominating mechanism, that is, the quadrupole relaxation of 2H2O. A formulation is sought that incorporates information from both MD and BD simulations into the NMR relaxation and line shape model developed in the context of BWR theory, eq . We follow the extensive body of theoretical and experimental work done on water relaxation in bilayer systems.[29,30,34,35] First, we need to distinguish the modes of MD that cause quadrupole relaxation (as well as broadening of the line shape) from the processes that are static on the NMR timescale, which result in quadrupole splitting of the spectra of locally ordered systems.[29] To address this, the ESI-1 provides detailed derivations of both static, denoted as, residual electric field gradient (EFG) anisotropy(29) and its time correlation function (TCF).

On the basis of the assumption of statistically independent processes often assumed for the system with the local order,[29,34] the static EFG may be represented in terms of separate ensemble averages. Referring to Figure , the processes relevant for water frame reorientation are (i) the local water dynamics ΩdM; (ii) the exchange between different water sites, which modulates the orientation Ωnd of the multiple water site director systems d (and d′) relative to the bilayer normal n; and (iii) the large-scale angular dependence of the bilayer normal, Figure c, ΩDn. This allows the calculation of the static second-rank EFG (see eqs S1–S4 in the ESI-1) aswhere the second-rank order parameters are given by SXY = ⟨D(2)(ΩXY)⟩XY computed as ensemble averages (denoted by the angular brackets) of the Wigner rotation matrix component Dmk(2) for orientational distribution related to the two reference systems XY. The parameter A0MP represents the tensor value averaged in the molecular frame (see eq S2), and V0P is the principal frame component.

The static tensor is directly proportional to the quadrupole resonance frequency, in eq ,and the quadrupole splitting observable in a spectra of an ordered system (isotropic powder average),[29] with the additional dependence on the fraction of perturbed waters PB

Both eqs and 7 include the quadrupole coupling constant, given by χ = |2eQA0MPV0P|/(6h), including the nuclear quadrupole moment (eQ) and Planck’s constant (h).

Time Correlation Function

A model for the TCF of perturbed water is derived introducing the slower membrane dynamics in a setting that makes it directly comparable to several previous models[15,29,30,35] (with details in ESI-1). In particular, it enables using the structural and dynamics information drawn from atomistic MD.[15] This model assumes a statistical independence of the three processes ΩDn, ΩdM, and Ωnd. Furthermore, two separate relaxation mechanisms are considered: (1) water undergoes partial reorientation (ΩdM) at the local site and is modulated by (ΩDn), that is, with the associated reorientation (Ωnd) fixed, and (2) water molecule undergoes site exchange on an averaged membrane interface. This leads to total TCFwhere gXY(τ) is the TCF normalized to [0, 1] [see eq S9 in ESI-1]. In eq , three types of membrane-perturbed dynamic modes of the water molecules are represented (see Figure ). The first term is decaying slowly, in a sub-microsecond timescale, because of the surface undulations and diffusion along the interface which modulate the environment experienced by the water molecule. The second term, on the other hand, is rapidly decaying (in the picosecond timescale) because of the local site reorientation (ΩdM). The third term encompasses water exchanging to various sites (d, d′) in a sub-nanosecond process suggested previously as based on a MD simulation.[15] The format in eq is for the analytical (approximate) powder average, obtained after averaging over the ΩLD Euler angles. This is complemented with the exact powder average by sampling the set of gLn(τ) TCF over the isotropic ΩLD distribution, instead of gDn(τ). Finally, the fourth term brings only a minor additional contribution that, similarly to the second term, decays rapidly at a picosecond timescale [see eq S11 in ESI-1].

The spectral density components are computed from the Fourier–Laplace transformation of the TCFswhere q represents the quantum transition order, and from eq , the total spectral density, including the free (bulk) water contribution (superscript F), is given bywhere PB and PF = 1 – PB are the fractions of bound and free water, respectively.

The magnetization and transverse relaxation rates, as well as the relaxation matrix components, are expressed as

Concerning the underlying approximations, the BWR theory requires that , that is, the characteristic correlation time of the slow dynamics (τc) is shorter than the inverse of the corresponding spin Hamiltonian amplitude. For the systems of this study, this limit corresponds to less than circa 10 μs. This requirement may be eliminated by employing a nonperturbative computational approach as has earlier already been exemplified for a similar 2H2O-surfactant system.[20] Given the computational effort in the presented BD simulations, this is expected to be a readily feasible approach. Furthermore, the time dependence due to exchange between water sites, as well as the complete QC tensor, may, in principle, be computed directly from atomistic MD snapshots,[25] thus avoiding the approximation of using a single experimental quadrupole coupling constant (χ) as done here. Such more advanced calculations can be expected to further enrich the information that can be drawn from NMR relaxation experiments. However, the approximations used in the present work are tolerable for the first characterization of the new BD/MD simulation method presented in this paper.

Flexible Membrane Dynamics

To introduce the undulating interface degrees of freedom [see Figure c] in a membrane, which is on the average flat and has no overhangs (i.e., 1–1 map between the position and the interface amplitude), we used the Fourier space Brownian dynamics (FSBD) method developed by Lin and Brown.[23,26] The work of Helfrich[36] for an incompressible fluid sheet was followed, relating the curvature of the membrane to the free energy. A linearized model for the free energy iswhere L is the projected side length of the square-shaped interface, L2 is the area, kc is the bending modulus, and r is the xy position in the projected plane. The Laplacian and gradient are taken from the function h(r), the height above the projected reference interface. Finally, Hint represents arbitrary external interactions, for instance, membrane–membrane interactions or the presence of membrane proteins.[23] A tensionless membrane was considered in this work (σ = 0).

The time evolution of the system may be formulated via a nonlocal Langevin equation for a two-dimensional (2D) periodic interface[23] and is conveniently solved for in the FS aswhere the 2D reciprocal space vector = (l, m)2π/L, with the integers l and m adapting values in the range −N/2, ..., N/2. Λ = 1/(4μk), with the wavenumber (k = ||), is the mobility (Oseen) tensor.[26] This is in the diagonal form relevant for small undulations.[37] The viscosity of the surrounding fluid is given by μ, and f is the force per area (drift term).[26]V(t) is the Gaussian white noise satisfying the fluctuation–dissipation theorem[37]and kB and T are the Boltzmann constant and temperature, respectively. We considered a free membrane where the Fourier modes decouple and with the drift term given by f = −kck4h,[26] where h is the Fourier amplitude of the height h. This choice (with Hint = 0) is done for consistency with the MD simulation[15] that does not explicitly include multiple membranes. The linearized curvature model[23] of eq is widely used since its first introduction in 1973 (ref (36)) and may appear as a restrictive assumption. However, it should be noted that the Hint term may also be used as a correction that accounts for the nonlinear curvature. This possibility will be explored in the forthcoming work.

The key steps of the numerical propagation of the FS (eq )[23,26] are provided in the following. The FS is discretized and truncated to the range given by −N/2 ≤ {l, m} ≤ N/2 – 1, thus giving the shortest real-space length scale at the interface equal to a = L/N. This is already a characteristic of molecular length scales, sufficient in the present case as explicit molecules are not included in the FSBD simulation. The desired simulation output is the real-space interface height function h(,t), and thus, FS is explicitly filled according to the conjugate symmetry (h* = h–). The key points of the simulation are (i) the initial FS structure of h(0) is given by a complex Gaussian random variable with zero mean and variance ⟨|h|2⟩, for which the analytical form is known for a free membrane;[37] (ii) the Fourier components are propagated according to the Euler numerical integration scheme of the first order in the time step[38] (vide infra); and (iii) the real-space height function is extracted with the inverse 2D Fourier transform

In the ESI-3, we provide details of the implementations (Octave[39] and MATLAB[40]), together with a verification that the simulation reproduces the analytically known TCF of h(,t).

Itô Diffusion on the Time-Dependent Interface

Static or Fluctuating Interface

Reference (17) introduced an efficient method to perform BD simulation on a curved surface, for instance, to model a spin-bearing molecule propagating along a lipid/solvent interface. In the previous work, NMR relaxation applications assumed a static interface[17,18,20] in modeling highly curved structures, in which undulations can be reasonably neglected. The BD of the translational process resides in the time dependence of in , with the model surface defined by the implicit representation ϕ() = 0. By definition, the dynamics takes place at the 2D interface. In the present work, the specific formis used, with the dependence on the height function h emphasized with the superscript. Hence, in general terms, the problem is a stochastic process on an s-dimensional interface (s-submanifold) embedded in s + 1.[17]

The mathematical form is a Langevin equation of dimension s + 1where the s + 1-dimensional vector is the Brownian motion in satisfying the fluctuation–dissipation theorem,[17]() is an (s + 1) × (s + 1)-dimensional matrix function, and () is again an (s + 1)-vector function.

Considering a Brownian process defined by the differential dϕ() = 0,[17] that is, where a particle initially at the surface will remain at the surface through all future times, leads to the following and in eq Here, ∇∇Tϕ is the Hessian matrix of ϕ (with transpose denoted by superscript T) and D0 denotes the translational diffusion constant at the interface. The gradient of the interface model is ∇Tϕ = [−∂h(), −∂h(), 1] (with ∂ ≡ ∂/∂x(), and || || denotes the Euclidean norm. Hence, projects and scales the three-dimensional (3D) Brownian process () to the interface tangent of h, and accounts for the surface curvature by including the Hessian matrix and projects parallel to the normal of the interface. A numerical scheme for eq is available[17] and has a sound mathematical foundation in the calculus of stochastic processes.

It should be noted that a 2D alternative to the presently employed, embedded method is possible. In this case, a parameterized surface provides the Langevin equation expressed with the metric tensor and local curvilinear coordinates.[41] The two eqs and 21 of the present method are advantageous in the sense that no curvilinear coordinates are required, making the present approach easier to work with. Furthermore, an implicitly defined surface, that is, ϕ, generalizes to more complex model surfaces.[17] The numerical scheme that follows from eq is more cost-effective in comparison with the empirical approach[17] that is suggested by Hołyst et al.[42]

Concerning the numerical implementation, we note that the gradient and Hessian in real space are conveniently available at the discretized surface by differentiating the inverse Fourier transform representation [see eq ]. The values at an arbitrary location can be obtained by interpolation. The numerical details for these steps are given in ESI-3.

Time Propagation Algorithm

A temporal discretization scheme was used for the coupled stochastic differential eqs and 20, considering either a static or time-dependent interface h, as

Hence, the interface model h is propagated in the FS, and the height function h(,t) at t = (p + 1)Δt, that is, h, is extracted with a 2D inverse Fourier transform [see eq ]. The interface curvature information incorporated in and allows the location of the spin-bearing molecule to be propagated in time, that is, from the time point p to p + 1. The key to a relevant relaxation study is that surface fluctuations feed into the translational process. The opposite, that is, translational process influencing the surface dynamics, is not considered relevant because of the large difference in mass. Time steps Δt of equal length are used with Δ = – . The integration algorithm for is of the predictor-corrector form, accurate to the first order in the time step.[38]

Simulation Setup and Parameters

A small-size bilayer MD simulation captures the molecular-level motion and small-scale collective behavior, as well as sets a natural upper limit on the parameter a = L/N. Typical box sizes in MD simulations are in the interval 4–12 nm,[43] here a was set to 7 nm, same as the discretization in previous BD simulations.[23] In this work, the spectral densities from ref (15) are incorporated [see jF, jnd, and jdM in eq ]. The MD simulation involves a 100 ns long trajectory of 64 DPPC lipids at the water/lipid ratio of 23, performed at a constant pressure and temperature of 1 atm and 50 °C, respectively.[16] From the trajectory, perturbed waters are defined as those with nonzero order parameter and the remainder denoted as bulk water. To determine bound water order parameter, as well as TCFs for various regions, the MD-TCFs are interpreted in terms of a model consisting of a sum of exponentials providing S0dM = 0.184 and the remaining perturbed waters have an order parameter S0nd = 0.027.[15] The spectral densities [jF, jnd, and jdM] follow from integrating the exponential TCF model.[15] The residence time of water in the two regions is estimated to be on the picosecond timescale.[15] Because of this fast timescale, the exchange-related relaxation mechanisms, caused by changing interaction tensors, are not included beyond the effect of reorientation caused by exchange. Concerning D0, the bound water is assumed to follow the DPPC lipid, with the experimental translational diffusion constant at 1.6 × 10–11 m2 s–1,[44] in accordance with MD predictions.[16] Note that the faster diffusing waters in the vicinity of the bilayer[15] have a relaxation contribution contained in the jnd spectral density within this model. The experimental heavy water quadrupole coupling constant in bilayer, χ = 220 kHz, is used.[29] The bending modulus of the bilayer may either be experimentally determined or extracted from MD simulations; here, we adopt kc = 4 × 10–20 J based on MD.[43] The temperature in the lamellar-phase experiment and MD simulation is T = 323 K.[35] The viscosity of the 2H2O medium that surrounds the membrane is η(T) = 0.651 mPa·s.[45]

An important question for periodic models in general is how large a system (here parameterized by the simulation cell side length L) is required to represent the relevant physics. With the FSBD simulation defining the system size, the atomistic MD simulation may be kept small because the underlying idea of the model is to account only for the fast processes with the MD simulation and the slower ones with BD. In contrast, a model containing only a MD simulation would feature no computational tool to predict how large a system size would be needed to capture the relevant physical phenomena. Such a model cannot presently be extended to be large enough.[43]

The situation is drastically different when conclusion from the TCF computed with BD simulations assist in determining the largest relevant system size, as is done in section . Hence, a significant progress is made in this work, in determining the system size based on model criteria. However, the choice of the system size also fixes the interface undulation amplitude with consequences for the quadrupole splitting [see eq ] and the amplitude of surface BD contribution to the NMR relaxation model [eq ]. Hence, also, the experimental quadrupole splitting in eq is considered, when selecting L. Here, L = 224 nm is used in production runs with S0Dn = 0.922.

In taking the modeling one step further, by explicitly simulating a stack of interacting membranes, it is noted that at large enough system size, the persistence length[12,46] provides a natural truncation of the longest-wavelength transverse modes without involving experimental data.

The production runs were done with the time step of 4 ps. Typically, 400 independent 8 μs long interface trajectories were used. The simulation is trivially parallel with one 430 MB trajectory produced by each processor, which requires about 1 h 50 min of cpu time in the MATLAB[40] implementation and, with the same implementation, about 8 h in Octave.[39] Octave is a freely available platform designed to provide the same functionality as MATLAB. Additional, but not dominating, computational cost is needed for powder angle sampling,[20] where each trajectory is rotated about 50–100 times. The computational details are discussed in ESI-3.

Results and Discussion

This section is arranged as follows. From the computed mean-square displacements (MSDs) of the NMR nuclei, in eq , it is illustrative to first compare the presented algorithm with previous simulations and analytical theory. Subsequently, the general aspects of the TCF gDn [see eq 8 and ESI-2] obtained from FSBD are discussed. These features are independent of whether the aim is to build a model for water, lipid, or other spin probes. Finally, the specific insights gained for 2H2O relaxation are discussed.

Obstructed Translational Diffusion

The effective translational diffusion constant, measurable with field-gradient NMR experiments,[44] may be related to the MSD aswhere the dimensionality of the system is d = 2. In the particular case of translational diffusion projected to the membrane plane, eq relates in a nontrivial way to the true curvilinear diffusion constant D0, where the latter quantity is of interest when interpreting experiments. For instance, in the case of diffusion and relaxation studies of systems with a complex geometry, the obstruction factor D/D0 plays an important role in understanding the results.[18] For a bilayer system at the limit where the membrane undulations are much faster than the translational diffusion process (this situation denoted the annealed limit), the obstruction factor for a free membrane is theoretically known as[12,47]where n⊥ denotes the xy component normal to the membrane plane. Although the annealed limit does typically not occur in a fluid bilayer system, eq has still been shown to coincide with simulations even at the other extreme, in which the surface dynamics is quenched.[12,27,28] This occurs provided that harmonic modes and a relatively stiff membrane are considered. Thus, a comparison shown in Figure of the MSDs computed with the algorithm of eq with eq illustrates consistency with previous observations.[12,27] The cases of either employing the DPPC parameters or those of a more flexible surface are illustrated in panels (a) and (b) of Figure , respectively. The obstruction factor is small (3%) for diffusion along the stiff DPPC surface, suggesting that from the modeling perspective, the use of an explicit BD simulation to compute the obstruction factor is of minor importance. From the experimental point of view, the key observation is that 2D rather than a 3D diffusion process takes place and atomistic MD simulation is the tool to understand the translational diffusion of the lipid and waters.[15,43] Furthermore, the result also highlights the complication mentioned in section , that is, the system size (L) must be set prior to interpretation of the model [see eq ]. This question is ignored in previous works on undulating membrane simulations[23,27,28] that, therefore, do not directly connect with experiments.

Figure 2

MSD for translational diffusion of water on an undulating lipid surface (black circles), DPPC Lα lipid diffusion constant D0 = 1.6 × 10–11 m2 s–1 ref (44), theoretical MSD for a free particle (dashed blue line) and the slope predicted by the analytical model (red line). (a) DPPC parameters used with L = 224 nm, with the analytical model giving D/D0 = 0.97. (b) With kc reduced to one-quarter and L = 112 nm, with the analytical model predicting D/D0 = 0.90.

The information required to compute the TCF gDn(τ) [see eq S9 in ESI-1] is the orientation ΩτDn at position τ and time τ, given the orientation Ω0Dn at position 0 and time 0, as an ensemble average over all possible configurations. This is the information accessible from the Langevin equation that maps the surface trajectories ([]) into the second-rank TCF. Equivalent information[38,41] resides in the conditional probability P(ΩτDn,τ|Ω0Dn,0) (lower case is used for deterministic functions). However, to obtain the latter, a high-dimensional Fokker–Planck equation needs to be solved, which is a difficult task and not readily extendable to the case of interacting surfaces. This makes the Langevin approach the more suitable approach. Using a continuum description, a model for free as well as interacting membranes has been constructed.[12,13] This involves approximations, however, making the presently introduced simulation method a suitable tool to test some of them.

The continuum description[12,13] builds on the obstruction factor [see eq ] and, thus, approximates the TCF calculation with a translational process independent of the time-dependent surface orientation: P(ΩτDn,τ|Ω0Dn,0) → Θ(ΩτDn|Ω0Dn)Ψ(τ|0), with conditional probabilities Θ and Ψ. In addition, an approximate local normal of the bilayer is required to compute the TCF. This means that from the exact unit vector along the direction of the normal, n(t)/|n(t)|, the components ñ⊥(t) ≈ −(∂xh(t),∂xh(t)), valid for |ñ⊥(t)| ≪ 1, are used. In addition, the approximation ⟨ñ⊥(0)ñ⊥(τ)⟩/⟨ñ⊥2⟩ is employed for gDn(τ). This type of separation has been in the past exploited for both a planar[12] and micellar surface.[21] It has also been indirectly assumed when considering the director fluctuations only, in the relaxation model.[9−11,48] Within the trajectory-based simulation approach, this approximation may be tested for free 2D Brownian motion with the translational diffusion constant D given by eq and independent realizations of the surface, h(t), denoted [h(t), D], to mimic the Brownian motion of the spin probe along the surface.

In Figure , the TCFs at director angles 0° [panels (a,c)] and 90° [panels (b,d)] are given for the exact process [], with the separated form [h(t), D] (dashed red line) and with ⟨ñ⊥(0)ñ⊥(τ)⟩/⟨ñ⊥2⟩ (dashed green line). Note, first, that there is a significant difference in the decay rate depending on the static director angle [compare panel (a,b)], hence, illustrating that the anisotropy of the TCF cannot be neglected. Thus, to model relaxation rates in an isotropic powder sample, explicit average over the director angles is necessary. Second, both approximate forms of TCFs provide an accurate representation at a director angle 0°, whereas at 90°, the deviation is greater. The agreement seen in panels (a) and (c) may be explained by the parameter regime of the system, that is, a stiff interface: given the harmonic surface model, both the separation of translational diffusion and surface undulations, as well as the approximated normal, work well.[12,21] However, we require TCFs for all possible director angles, and panel (b) illustrates that the first required approximation, separation of translational diffusion and surface undulations, introduces an error that cannot be ignored. This deviation in panel (b) corresponds to a 13% overestimation of T2 at high field (for this particular orientation).

Figure 3

Comparison of the second-rank TCF [gLn(τ)] for the exact process (blue line) with the result of approximate TCF [h(t), D] (dashed red line) and with ⟨ñ⊥(0)ñ⊥(τ)⟩/⟨ñ⊥2⟩ (dashed green line) for ΩLD = [0, 0, 0] [panel (a,c)] and ΩLD = [0, π/2, 0] [panel (b,d)], respectively. The inset shows the difference TCF(approx)–TCF(complete). The statistical error is estimated to be 1 × 10–3.

Typically, when the TCF is approximated to the second order in the surface normal, an assumption that the shape of the approximate TCF is independent of the director angle (ΩLD) is used.[9,10,48] Hence, the full anisotropy of the TCF is not accounted for in this approximation,[12] as it is in the simulation approach. A potential remedy for this is to go to a more complex TCF model.[13] However, given the interest in this type of model,[9−11,21,48] the same TCF is used in panels (c) and (d) of Figure , to explore this assumption. At the director angle 90°, this leads to an overestimation of T2 by a factor of 3. Thus, Figure highlights the need for simulation methods to make progress in NMR relaxation modeling for D2O, as well as for bilayer systems in general.

Figure shows the outcome of trajectory-based simulations, illustrating the TCFs at various system sizes for an undulating surface [h(t), D = 0] (i.e., using a static position of the spin system), for a translational spin system on a static surface [], and for a dynamic spin carrier on a fluctuating surface []. The TCF for the last case mentioned [] is also computed both with the DPPC lipid diffusion constant and by a factor of four slower translational Brownian motion. An approximate analytical model for the tail of the TCFs, τmax = 8μL3/[(2π)3kc], derived in ESI-2, is included in the [h(t), D = 0] case. Several observations can be made. The case with explicit surface diffusion [panels (b,c)] provides TCFs approximately with the same decay rate for large enough system size (L). The decay process is more efficient with the combined surface undulation and translational diffusion, in panel (c) as compared to diffusion on a static surface (b), and in the former case, the system sizes from 112 nm and larger only show TCFs approximately decayed after 1 μs independent of L. In panel (a), the opposite and significant simulation cell size dependence is present (note the logarithmic scale) and the TCF has not at all decayed after 10 μs for the larger systems. The tail model, τmax, reproduces the simulated results in this case with a static spin carrier on an undulating surface, which shows an unbounded (τmax ∝ L3) dependence on the system size. Panel (d) illustrates that the slope of the TCF tail is proportional to the translational diffusion constant (and not to the τmax surface undulations), for the process .

Figure 4

TCFs gDn(τ) for (a) surface undulation model [h(t), D = 0], (b) surface translational diffusion on the static surface [] and (c) translational diffusion on the flexible surface [], with the system size L in the rage 28–448 nm. The tail in panel (a) is the derived τmax [see eq S19 in ESI-2]. Panel (d) displays TCF with the DPPC diffusion constant (D0 = 1.6 × 10–11 m2 s–1, blue line) and with D0/4 (dashed red), L = 224 nm.

2H2O Relaxation

Here, we focus on the three main relaxation model assumptions that are compared with high-field experimental data for the hydrated DPPC bilayer powder sample (22 waters per lipid).[35] The resulting line shapes and relaxation–dispersion curves are given in Figure and Table (relaxation rates). For all three models, we include the spectral densities simulated by atomistic MD. The three models are: (a) pure atomistic MD simulation [h = 0, D = 0], (b) flexible surface, [h(t), D = 0], and (c) translational diffusion on a flexible surface, []. The theoretical line shapes and the field dependence of the relaxation rates are given in the left and right columns of Figure , respectively.

Figure 5

Line shapes and magnetic field dependence of relaxation rates (T1–1, T2–1) for three relaxation models, all based on the spectral density from MD: (a and a′) pure atomistic MD,[15]h = 0, D = 0; (b and b′) static spin carrier location, flexible (FSBD) surface [h(t), D = 0]; (c and c′) translational diffusion on the flexible surface, . The dashed black line illustrates the ω0–1 dependence for comparison. High-field (9.4 T) experimental relaxation rates are represented in (a′), (b′), and (c′) as a triangle (T2–1) and a circle (T1–1).

Table 1

Experimental 2H2O NMR Relaxation Rates in the DPPC/Water System (22 Water Molecules per Lipid, Field 9.4 T),[35] as well as Rates Simulated with BD and MD Spectral Densities

	T1–1 (s–1)	T2–1 (s–1)	T2–1/T1–1
(1) experiment	3.57	250	70
(2) h = 0, D = 0a	2.59 ± 0.06	2.59 ± 0.06	1.0
(3) h(t), D = 0	2.75 ± 0.04	540 ± 13	196
(4) Xth(t)	2.6 ± 0.1	137 ± 6	53

Pure MD relaxation rates.

First, it is noted that all three approaches predict the high-field T1–1 reasonably well. This independence of the employed surface model is expected because this rate is mainly determined by the fast (with timescale 1/ω0) molecular processes at high field. The exceptional challenge to the model is to interpret the T1–1 relaxation dispersion and the transverse relaxation rates. The second row of Table shows a transverse relaxation rate that is 2 orders of magnitude too slow, when pure MD is considered.[15] This is also noted from the line shape[15] in Figure a, which is by far too narrow to reproduce the typical experimental line width for heavy water.[34,35]

Under the assumption that the flexible interface is the sole phenomenon that provides the slow dynamics, illustrated in Figure b′, it is seen that the linear field dependence[9,10,48] expected for this dynamical model is produced by the simulation. Within this model, the absolute agreement with the experiment is also improved. However, the spin–spin relaxation rate is found to be more than a factor of two too large (the third row in Table ), which is also visible from the too broad line shape.[34,35] To approach experimental results, a low-frequency cutoff of surface undulations would be needed. Such a truncation may serve as an additional adjustable parameter and requires a theoretical motivation.[49] All this renders the undulation model [h(t), D = 0] a poor description. With the complete process, , a slight deviation from the linear dependence is seen at high frequency (panel c′). The majority of the model contribution to T2–1 comes from jDn (averaged over βLD powder angles), with 2% from the exchange (jnd) and local site reorientation (jdM) combined. Translational Brownian motion [Figure c′] significantly slows down the NMR relaxation rates and is a key to provide experimental relevance of the model. Within the complete process , there is no apparent need to introduce a low-frequency cutoff of surface undulations. Comparing the case (fourth row in Table ) with experiment again reveals a significant improvement as compared to purely atomistic MD interpretation, going from by a factor of 100 to a factor 2 in underestimated T2–1. The width of the line shape can be visually verified to be in the regime of experiments.[34,35] Furthermore, the line shapes (a) display a small dip at zero frequency (βLD = 57.4°); a feature present when significant MD occurs in the timescale regime of 1/ω0.[30] Because this feature is not observed for spectra with a fully hydrated bilayer[34,35] and is not observed in the line shape (c), it gives support of the model.

To summarize, three types of waters are considered [see eq ]: (i) waters bound to the membrane experiencing a fast local reorientation[15] (jdM) and are modulated by translational diffusion (at the rate of lipid diffusion) as well as membrane undulations (jDn); (ii) MD simulations dictate that there are waters diffusing much faster than lipid diffusion,[15,51] and these are included in the jnd spectral density computed from the MD simulation; and (iii) finally, bulk waters (jF) are modeled in the MD simulation.

A refined water relaxation model may be considered in view of the derived TCF in ESI-1. For instance, on the NMR timescale (milliseconds), waters may move several hundred nanometers, suggesting that a model where also diffusing waters in jnd are coupled to the slower modes of the membrane would be a useful route to explore. This would complement postulated model,[30] or for instance, “reorientations mediated by translational displacements”.[52] In addition, packing defects of the lipid bilayer may come into play by influencing T2–1. A simulation model may in the future explore this with a concentration-dependent relaxation sink. Similar approach may be used to model lipid relaxation.[49] Given the additional experimental information provided for oriented bilayers[4,49,50] and field-dependent relaxation studies,[4] such experimental studies will clarify which additional mechanisms are needed.

Conclusions

For more than four decades, the dynamical information from NMR relaxation studies on lipid bilayer systems, relevant on the biological timescale, has presented a challenging modeling topic.[9,10,12,29,34,48] However, during this period, atomistic MD simulation has only played a minor contribution in confirming the static properties for the high-field magnetization relaxation, whereas the slow timescale dynamics relevant for the transverse relaxation rate has not been studied at a quantitative level of accuracy, in general.[14,15]

Proposed here is a mesoscopic simulation approach incorporating atomistic MD results to model NMR relaxation in fluid bilayer systems. An example is given for solvent relaxation in a model of DPPC in 2H2O. A division is made between the fast MD of the spin-carrying moiety (water in the example case) and the, assumed statistically independent, slow processes (surface undulation and translational diffusion) modeled with the Brownian motion. In particular, the molecular order and fast motion of the spin moiety, be it solvent, lipid, or some other spin probe, are accessible from atomistic detail available from MD simulations.[14,15,19,25] The slower processes are modeled with BD simulation, and they consist, in the example case, of reorientation by translational dynamics. In this, the translational motion is intertwined with interface undulations simulated with Fourier space surface dynamics. This allows for a close connection to both analytical models[9,10,12,29,34,48,52] and pure MD modeling.[14−16,53] The key novelty of the presented simulation methodology is that the surface structure and dynamics explicitly feed into the translational motion, thus avoiding the assumption of statistical independence of the two processes.[12,21]

In the relaxation test case at hand, the performance of the combined MD and BD model shows significant progress as compared to pure atomistic MD, explaining within a factor of two in the transverse (T2–1) relaxation rate. This suggests that the combination of interface undulation and translational diffusion provides a plausible explanation for the apparent shortcomings in the pure MD simulation. It is found that translational diffusion plays a significant role in the relaxation, which cannot be ignored. In contrast, for unilamellar micellar systems, considering the dipolar relaxation mechanism of lipid molecules, it is argued that surface undulation alone is the dominating mechanism, whereas translational diffusion can be ignored.[21] The present results for quadrupole 2H relaxation motivate us to apply the combined simulation approach to a number of different soft matter systems.

Similar to the previous works,[27,28] translational diffusion is seen to agree well with an analytical model for a free membrane (see Figure ). Hence, it can be questioned where a BD simulation for translation along the surface is really needed. However, we show that the same simplifying assumptions introduce deviation in the relaxation model (compare Figures and 3), suggesting that NMR relaxation with the simulation method is an important and sensitive tool to explore membrane dynamics.

Comparing the BD/MD simulation methodology with conventional analytical models,[9,12,48] we see a great benefit in the spirit of Occam’s razor: a drastic reduction of the number of, or no adjustable parameters at all. However, this is not first principles modeling, either. Thus, the hypotheses made can serve as parameters. A plausible explanation for the remaining discrepancy from the experiment may be that even though a fully hydrated membrane at large (28 Å) interlayer separation is considered, entropic, long-range undulation forces are present.[1] Thus, extending the simulation methodology to a stack of explicitly interacting membranes, both at BD and MD levels and comparing with relaxation experiments at different hydration levels, will be of interest in the future. In contrast to the continuum approach, the BD simulation can be extended in a relatively straightforward manner to interacting membranes, even at low hydration levels, and/or to corrections for linearized free energy [see eq ]. Such refinements of the simulation methodology would allow transcending beyond the approximation of from harmonic surface undulations. Furthermore, concerning the spin dynamics, the approximation of a single static quadrupole interaction tensor can readily be refined with the corresponding time-dependent tensor available from the appropriate quantum chemistry property surface.[25]

The presented and prospected developments are needed in expanding the exploration of soft matter via spin dynamics observables, to address at a great level of detail many biologically relevant problems.[2,3] Further applications of the model, together with what alternative spin systems, such as long-lived states,[25,54] are envisaged for future work.