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.


Contents
List of Figures

ESI-1. Electric field gradient
A summary is given on the electric field gradient (EFG) at the deuterium nucleus of heavy water molecules [1], and how it enters the NMR observables, i.e., the static (or residual) EFG anisotropy and its time correlation function (TCF).
Residual electric field gradient anisotropy. The second-rank spherical component of the EFG tensor is expressed using the various coordinate frames in Figure 1 of the article, considering the successive transformations L → D → n → d → M, as V L p (t) = m,k,l,q (−1) q D (2) pm (Ω LD )D (2) mk (Ω Dn t )D (2) kl (Ω nd t )D (2) where D (2) mk (Ω XY ) denotes the second-rank Wigner rotation matrix component [2] for frames X and Y, with the indices {m, k} in the range −2, −1, 0, 1, 2. The transformations follow from the closure and symmetry rules [2]. First we need to distinguish the modes of molecular dynamics that cause quadrupole relaxation (as well as broadening of the lineshape) from the processes that are static in the NMR timescale, and which result in quadrupole splitting of the spectra of the locally ordered systems [1]. This implies that an assumption of statistically independent processes is made [1,3]. Referring to Figure 1 in the article, the processes relevant for the reorientation of the water molecular frame are (i) the local water dynamics Ω dM t ; (ii) the exchange between different water sites Ω nd t modulating the orientation of the water site director system d relative to the (local) bilayer normal n and (iii) the large-scale motion of bilayer normal, [see Figure 1 (c)]. This allows the calculation of the static EFG as with the second-rank order parameter given by S XY ij = D (2) ij (Ω XY ) XY computed as ensemble averages denoted as · XY for orientational distribution related to reference systems XY and A MP q A P 0 defining the average in the molecular frame, where V P 0 = ( √ 6/2)V P zz is the dominating principal component of the EFG. The molecular frame average, which is typically one of the adjustable parameters in experiments is, however, also open to computation by combining atomistic MD with quantum chemistry property surface (QCPS) [4], the latter to reduce the cost by eliminating a large number of QC calculations. To identify "bound" water, responsible for both the quadrupole splitting and the large T −1 2 rate seen in these systems, the water region with non-zero second-rank order parameter in an atomistic MD simulation [5], is selected. It should be noted that the order parameter is defined relative to a static bilayer centre in Ref. [5]. The eq S2 may be simplified further by assuming that there exists at least a three-fold symmetry in the partial ensemble averages (the is supported by the MD simulation subject to finite sampling [5]), giving S dM lq = δ l0 δ q0 S dM 0 (δ l0 = 1 if l = 0, otherwise zero) and similarly also for the Dn and nd reference systems, where the former follows from the cylindrical symmetry of the membrane interface. Hence, with these symmetry properties eq S2 may be written as V L,static In this work, the quantity A MP q V P 0 (with q=0) is assumed also for TCF calculation.

S2
Time correlation function. The TCF for bound water [1,6] (denoted by the subscript B) is given for the EFG with zero mean [∆V L where · LM is the ensemble average over orientational dynamics and in the following we are interested in the TCF for an isotropic powder sample for which the forms V L p , p = 0, 1, 2 are equal. The TCF in eq S5 is not directly accessible to trajectory-based methods due to the six orders-of-magnitude difference between the timescales of local water dynamics and the slow collective membrane dynamics. In the following, an explicit model is derived, whereby the slow membrane dynamics is introduced, and which enables a direct comparison with the previously established relaxation models [1,5,6,7]. The TCF is formed by combining eqs. S5 and S1 and only in the initial discussion of TCFs introduce the approximate factor 1/5, [1] resulting from a powder average over an isotrope Ω LD -distribution: Note that the TCFs for different static Euler angles Ω LD have in general different decay as become apparent from the Ω LD -dependent weights of the TCF-components (see eq 2.3 page 692 in Ref. [9]). Thus, the approximate powder average is subsequently compared with the exact counterpart (obtained by explicitly sampling the Ω LD -distribution, [c.f. sec. ESI-2]).
To progress from eq S6, attention is given to the various timescales, in particular the slow processes that can cause significant quadrupole relaxation. It is first noted from the MD simulation [5], that there is a short (picosecond) residence time of the bound water sites. This motivates a TCF model where a water molecule is either bound or undergoes a site exchange process, but where the explicit site-residence time is too short to have impact on the NMR relaxation rates (other than in the configurational averaging). The physical model for the membrane-perturbed water molecule is based on considering the following two processes, (1) water undergoes partial reorientation (Ω dM t ) at the local site and is modulated by (Ω Dn t ), i.e., with fixed (Ω nd ), (2) water molecule undergoes site exchange on an averaged interface. Dividing eq S6 into two ensemble averages, one for both (1) and (2), gives, without "double counting" any of the processes, first Here and in the following, use has been made of the previously mentioned statistical independence of the three processes Ω Dn t , Ω dM t and Ω nd t (due to the timescale separation [1,3]). Also at least three-fold symmetry in the angular distributions has been considered [1], together with the orthogonality of the Wigner rotation matrix elements [2]. The set of Wigner elements for Ω nd in the first term of eq S7 equals unity by definition [2], S3 and applying the criteria leads to Introducing normalised TCF defined by and entering these into eq S8 gives Finally, expanding the bracket in eq S10, withG . The MD TCFs are adapted from Ref. [5] with the correlation times τ dM = 1.6 ps, τ nd = 110 ps, and the last product g Dn (τ )g dM (τ ) ≈ g dM,eff (τ ) is approximated using Dn , with these computed in the single-exponential format. The TCF for perturbed water is combined with "free water" TCF (0.74 ps correlation time) in computing the spectral densities and relaxation rates, in the main article.

ESI-2. Analytical properties of FSBD
The TCF of a fluctuating free membrane in the Fourier space is given by [8] where h k are the Fourier components of the interface amplitude and the real-space TCF follows from the inverse Fourier transform of [see eq 18 in the article] The correlation time is and it should be noted that a tensionless (σ = 0) membrane is considered in this work.
Involved are the viscosity of the surrounding medium µ, k= k , the Boltzmann constant k B , bending modulus k c , the side length of the square-shaped periodic membrane simulation cell L and temperature T . The relevant TCF, g Dn , for an NMR relaxation study concerning membrane fluctuations (Ω Dn t ) is given by eq S9. The TCF may, equivalently, be written in terms of a second-order Legendre polynomial, P 2 [as follow for a cylindrical distribution from eq 2.25, page 28 in [2]]: where P 2 [x] = 1/2 (3x 2 − 1) and the local surface normal is given by The corresponding TCF g Ln (τ ) for a static director orientation Ω LD is readily available by rotating the local surface normal [n(τ )] Ln = A(Ω LD )[n(τ )] Dn , where A is the Euler rotation matrix. The eq S15 cannot be connected with the analytical solution (eq S13) for the general case of undulating interface. This, motivates developing the present simulation methodology. The need for a simulation approach is emphasised when considering the generally non-separable translational diffusion process on an undulating surface. However, to connect simulation results to analytical theory, an approximate form in eq S15 is discussed in next section.
Approximate analytical TCF. As introduced in Refs. [9,10] the case of small fluctuations is considered where the local membrane normal may be approximated as In addition, as an approximation to second order in the fluctuations of the normal, the TCF takes the form [9]g whereñ ⊥ is the component of the normal in the base (x (1) , x (2) ) plane. In the long-time limit of eq S18, it is the slowest-decaying Fourier component h k , denoted h max , that remains as the relevant one. Furthermore, the surface normal in eq S17 is computed by multiplying the Fourier coefficients h k by the factor −ik a (a = x (1) , x (2) ) followed by the inverse Fourier transform. Hence, the components of eq S18 may be written as the product of a time-independent factor and h * max (0)h max (τ ) . The slope of the TCF tail, following from eq S12, becomes 1 From eq S19 and the definition of k it follows that the tail of the TCF has a correlation time that is equal to two times the largest component of eq S14:

ESI-3. Numerical FSBD implementation
To compute the discrete inverse Fourier transform [eq 18], as well as the surface derivatives and the Hessian of main the paper, particular care is needed considering the specific definition of the ifft routine within the programming language used. In Matlab [12] and Octave [11] used here, the function ifft2 is defined (for a square 2D matrix) as: , it is noted that the derivatives, and not the surface itself, are the relevant quantities. The overall scale factor of the surface amplitude may be chosen according to the dimensionality of the simulation cell. The critical part is the numerical accuracy of the derivatives and that the gradient and Hessian are correctly scaled relative to each other. In comparison with grid-based derivatives discussed below, the scale factor M = (N/L) 2 is introduced for the interface in the real space: The spatial derivatives, considering eq S20 and the scale factor M D = N/L, are given by ifft2 of the following variables: Note that the structure of the Fourier space is retained [see eq S22], thus the procedure does not require a shift of the output, unlike what is normally needed for numerical ifft2.
Grid-based gradient. In particular, the surface normal computed with the above steps and eq S16 is found to correspond to the grid-based implementation in the Matlab function [ nx , ny , nz ]= surfnorm (X, Y, h ) ; However, the latter contains a spatial discretisation error. Since this error is not present in the presently adopted Fourier-transform approach, it is more cost-effective and reduces the level of complexity for the user. . This may in principle be done by initiating a flat surface [h k (0) = 0] that is subsequently equilibrated [13]. In this approach, a significant time (of the order of τ max ) is spent in reaching the equilibrium, however. Instead, direct generation of complex Gaussian random variables for h k is performed, with zero mean and variance given by eq S12. It should be noted that the element k = 0 is a singularity in general. However, in the present application the center of mass S6 is fixed [13] and is not a dynamical quantity and, therefore, is set to Λ −1 k=0 = 0 (whereas Λ k = 1/4µk for the other grid points). (3) h k,p is propagated according to the first line in eq 22 of the article, in which ∆V k,p are complex Gaussian random numbers with zero mean and variance (4) The local interface normal [see eq S16] and second derivatives are obtained following from eqs S22-S27. Only the derivatives are needed for the NMR observables. Note that the explicit powder average is obtained at this point by sampling g Ld with Ω LD = [0, β LD , 0], drawn from an isotropic angular distribution and rotating each trajectory in eq S16 (see appendix A in [14]). (5) For the derivatives required for propagating the translational Brownian motion [see eq 21 of the main article] the above ifft2 results are interpolated to the current spatial location of the diffusing particle using the function qinterp2 [15] (vide infra).
Discretisation in simulation. In FSBD, two kinds of discretisation need to be considered: spatial discretisation of the surface and the time-step discretisation concerning the surface undulations and translational diffusion. The time-step discretisation of the surface undulation simulation and static position of spin system (see the case [h(t),D = 0] in Figure 4 of the main article) was tested by comparison with the analytical formulas in eqs S12, S13. In particular, min(τ k ) in eq S14 dictates how small time step is required in the surface-FSBD simulation. A cost-effective route to determine the relevant time step for the Fourier-space simulation is with a focus on this, fastest-decaying mode. For the DPPC membrane considered, min(τ k )=0.31 ns and ∆t = 0.004 ns [see eq 22 in the article] is found to be a sufficiently small time step. Both the simulated and analytical real-space TCFs of the interface (real space) amplitude are illustrated in Figure s1. This TCF is found to have significantly larger statistical variability than the corresponding TCF relevant for NMR relaxation [see eq S15]. Thus, this increases the statistical error of the TCF. However, within the statistical error, a good agreement between the simulated and analytical results is found.
Concerning the translational Brownian motion on a static surface, it was concluded in previous work [16] that it is the spatial variability, for instance minimal Gaussian curvature or, alternatively, minimal wavelength that dictates the maximal step size of time discretisation [see eq 22 in the article]. Here, no analytical solution is available in the general case. However, in time step regime where the discretisation error is dominating over the statistical error, the theoretically known, (O(∆t 1 )) truncation error of the numerical scheme has been verified considering |S Dn 0 (∆t) − S Dn 0 (∆t min )| for both frozen and undulating interfaces [16]. Furthermore, the TCF in eq S15, together with T −1 2 were inspected as a function of the time step. Focusing on a small surface model (N = 4) in exploring the time step of translational BD, reduces the computational cost and emphasises the larger curvatures and faster fluctuations that represent a bigger challenge for the choice of the time step. Typically ∆t = 0.004 ns is found sufficient for D 0 = 1.6 · 10 −11 m 2 s −1 and was also used for the case D 0 /4, i.e., in all the production runs.
For the translational process, the surface normal and Hessian are needed at an arbitrary spatial location. Hence, interpolation methods were used to obtain these, featuring refined spatial discretisation. Similarly to time discretisation, the spatial discretisation is expected to need refinement when the translational displacement (estimated as √ 4D 0 ∆t = 0.02 nm − 0.05 nm) is close to the minimal wavelength occurring in the simulation (ca. a = 7 nm). The precise meaning of "close" remains to be explored in the simulation however. The components of the spatial Hessian are obtained with the 2D-interpolation function qinterp2, a faster version of the Matlab interp2 developed for evenly-spaced spatial discretisation [15] and in this work also found to work also in the Octave platform. As a route to improve the accounted spatial variation, a linear interpolation was explored. In Figure s2, linear interpolation is used in the computed TCF (blue circles) and compared with the qinterp2 flag "nearest" (finding the nearest grid point), without significant difference to within the statistical error. Hence, the somewhat less costly nearest interpolation method is used in this work, with the discretisation of L/N = 7 nm.
The production runs were carried out with N = 32 and included both an undulating surface and translational Brownian dynamics. The trajectories were 8 µs long and take 1 h 50 min of cpu time each in MATLAB [12] and about four times longer in Octave (clock speed 2133 MHz, requiring only 430MB of memory). The Matlab-based FSBD script computes single-processor batch jobs for each surface trajectory. Typically 400  Figure s2. Comparison of the TCFs of g Dn (τ ) [see eq S15] computed both with the qinterp2-interpolation "nearest", that finds the nearest grid point (black line) and "linear" interpolation (blue circles). The "nearest" case is computed with 400 trajectories, each 10µs long, whereas the "linear" case is computed from one hundred 5 µs long trajectories.
independent surface trajectories where simulated, each computing four mutually noninteracting, translational Brownian processes on each surface.