Anomalous and anisotropic nanoscale diffusion of hydration water molecules in fluid lipid membranes

We have studied nanoscale diffusion of membrane hydration water in fluid-phase lipid bilayers made of 1,2-dimyristoyl-3-phosphocholine (DMPC) using incoherent quasi-elastic neutron scattering. Dynamics were fit directly in the energy domain using the Fourier transform of a stretched exponential. By using large, 2-dimensional detectors, lateral motions of water molecules and motions perpendicular to the membranes could be studied simultaneously, resulting in 2-dimensional maps of relaxation time, t, and stretching exponent, b. We present experimental evidence for anomalous (sub-diffusive) and anisotropic diffusion of membrane hydration water molecules over nanometer distances. By combining molecular dynamics and Brownian dynamics simulations, the potential microscopic origins for the anomaly and anisotropy of hydration water were investigated. Bulk water was found to show intrinsic sub-diffusive motion at time scales of several picoseconds, likely related to caging effects. In membrane hydration water, however, the anisotropy of confinement and local dynamical environments leads to an anisotropy of relaxation times and stretched exponents, indicative of anomalous dynamics.


Introduction
Diffusion is the fundamental mechanism of motion for lipids, proteins and water molecules throughout a membrane. While diffusion describes a stochastic, intrinsically non-directed motion, many biological processes are related to highly directed motions. The incoherent interaction between active processes in the cytoplasm was recently shown to lead to random fluctuation forces, 1 which in return may change the character of a directed motion into a thermally-induced diffusion process.
In biology, diffusion often occurs in crowded media, which for instance, have been found to lead to anomalous diffusion of proteins. [2][3][4] The dynamics of the atoms in protein or lipid molecules are more complex than those in ordinary liquids. In a simple liquid, atoms move ballistically at short times, followed by a crossover to Brownian diffusion. In dense fluids a ''caging'' effect is observed, where the atoms can be trapped momentarily by their neighbours. The motions of protein and lipid atoms are further complicated by their polymeric structure, characterized by connectivity and flexibility.
Sub-diffusion of lipid molecules in membranes has been observed in computer simulations. [8][9][10] The experimental observation of sub-diffusive lipid or hydration water motion, however, has turned out to be more difficult and a conclusive observation has yet to be reported. 11,12 While numerous studies addressed nanoscale lipid dynamics, see, e.g., [5][6][7][13][14][15][16][17][18][19][20][21] significantly less work has been done to study dynamics of membrane hydration water, most likely due to the corresponding computational and in particular experimental challenges. However, hydration water was shown to have distinct properties as compared to bulk water, such as a slower diffusion, 22 reduced freezing temperature 23,24 and glassy dynamics. 25 Properties of hydration water play an important role for the dynamics of biological systems up to the point that their corresponding dynamics have been suggested to be ''slaved'' by the solvent's fluctuations. 26 Protein dynamics and folding were reported to strongly depend to hydration water dynamics. [27][28][29] Lipid structural dynamics have a complex relationship with hydration water 30 and coupling of membrane and hydration water dynamics has been reported. 24,31 From computer simulations of lipid bilayers, sub-diffusive and ballistic regimes of hydration water motions have been predicted on short time scales and correspondingly small distances. 32,33 Water molecules are expected to move anisotropically when diffusing along the membranes or approaching and penetrating the bilayers, as reported from nuclear magnetic resonance imaging 34 and are shown to slow further from the introduction of proteins by sum frequency generation spectroscopy. 35 Time-resolved vibrational spectroscopy has shown that water in DOPC membranes is distributed anisotropically and forms nano-clusters 36 for a range of hydration levels, supported by molecular dynamics simulations. 37 While water molecules can penetrate lipid bilayers, 38 the perpendicular diffusion, D > , is typically at least an order of magnitude lower than lateral diffusion along the bilayers. Water motion through a bilayer is not a simple diffusion process as free energy and diffusion rate strongly depend on the position in the bilayer. 39 In this work we combined neutron diffraction, incoherent quasi-elastic neutron scattering and computer modelling to study nanoscale diffusion of hydration water molecules throughout lipid bilayers on length scales of B3-100 Å. By preparing oriented membranes, the experimental results give unambiguous access to the anisotropic local molecular confinement geometry to which the diffusive processes are restricted. Large, 2-dimensional detector arrays further enabled us to differentiate between motions of water molecules along and across the membranes separately.
Hydration water diffusion was also modelled using Brownian dynamics in a coarse-grained simulation of the trajectory of water molecules represented by center of mass coordinate for 10 ms (in time steps of 0.1 ps). Diffusion constants and potentials for water molecules are specific to the region, such as hydrophobic core and hydration water layer, with friction included as an effective term. From the computer model, we find that the anisotropy of confinement and local dynamical environments of membrane hydration water lead to an anisotropy of relaxation times and stretched exponents, indicative of anomalous diffusion.

Materials & methods
2.1 Experimental sample preparation 2.1.1 Bulk water sample. The bulk water sample consisted of a cylindrical aluminum annulus filled with ultra pure H 2 O with a resistivity of 18.2 MO cm. The walls of the annulus were separated by 0.1 mm and the can was held at room temperature (T = 21 1C = 294 K) during the measurement. Vanadium foil was bent to the same shape as the annulus and used to measure the energy resolution at the different incident neutron energies.
2.1.2 Membrane samples. Highly oriented, multi-lamellar stacks of 1,2-dimyristoyl-3-phosphocholine (DMPC) and chain deuterated DMPC (DMPC-d54) were prepared on 2 00 singleside polished Si wafers of thickness 300 mm. A solution of 16.7 mg mL À1 DMPC in 1 : 1 chloroform and 2,2,2-trifluoroethanol (TFE) was prepared. The Si wafers were cleaned by alternate 12 minute sonications in ultra pure water and methanol at 40 1C (313 K). This process was repeated twice. The cleaned wafers were placed on a heated sample preparation surface, which was kept at 35 1C (308 K). This temperature is well above the main phase transition for DMPC that occurs at 23.4 1C (296.4 K), thus the heated substrates ensured that the lipids were in the fluid phase during deposition and the self-assembly of the lipids. An aliquot of 1.2 mL of the lipid solution was deposited on each Si wafer and allowed to dry. The wafers were kept under vacuum overnight to remove all traces of solvent. Following this protocol, B3000 highly-oriented, stacked membranes were supported on each wafer with a total thickness of B10 mm.
Deuterated lipids, DMPC-d54, hydrated with H 2 O, were used for the incoherent experiment to emphasize the incoherent scattering of the hydration water molecules. In order to study water dynamics from the quasi-elastic neutron scattering experiments, eighteen Si wafers were stacked with two 0.3 mm aluminum spacers placed in between each wafer to allow for the membranes to be properly hydrated. Wafers and spacers were stacked into an aluminum sample can, as shown in Fig. 1(a)-(c). Aluminum foil between the wafers and sample can lid provided tension to keep the wafers in place. Samples were then hydrated with ultra pure H 2 O and incubated at 30 1C (303 K) for 24 hours before the measurements. Energy resolution was determined by cooling a sample with a large incoherent scattering contribution consisting of protonated lipids hydrated with H 2 O down to 15 K to freeze all molecular dynamics, such that the observed scattering was considered to be purely elastic.
In this experiment, the bilayers were hydrated from water vapor to avoid high absorption and background signals from bulk solvent hydration. 40 Solvent-saturated filter paper was placed above loosely packed balls of aluminum foil and cadmium shields were mounted carefully to avoid scattering of the filter paper and bulk water. The can was mounted inside of a temperature controlled cryostat with a temperature stability of AE0.01 K. Data were taken while the membranes were fully hydrated and in their fluid phase, at a temperature of T = 30 1C (303 K). Due to the typically low scattering of the sample (B10%), we did not observe evidence of multiple scattering events.
The membrane structure was determined via neutron diffraction in which a single Si wafer was mounted in an aluminum can, as shown in Fig. 1(d) and aligned in a temperature controlled chamber on a neutron reflectometer. Temperature was set to T = 30 1C (303 K) and controlled by a circulating water bath to better than AE0.1 K. Samples were hydrated with ultra pure H 2 O or heavy water (D 2 O), and incubated at 30 1C (303 K) for 24 hours before the measurements.
2.2 Elastic and quasi-elastic neutron scattering 2.2.1 Neutron diffraction on OFFSPEC. The structure of the bilayers was studied through neutron reflectivity at the OFFSPEC reflectometer at ISIS, the pulsed neutron and muon source at the Rutherford Appleton Laboratory in Oxfordshire, UK. 41,42 OFFSPEC is a vertical, flexible, low background, reflectometer and was used in its unpolarized mode.
We used OFFSPEC to probe the structure of the lamellar bilayers perpendicular to the plane of the bilayers. The thickness of the hydration water layer and water distribution throughout the bilayers, bilayer head group and tail region thickness were determined from these measurements. Neutron scattering has the advantage of scattering differently from hydrogen and deuterium atoms due to their vastly different scattering cross sections. This difference allows measurements to highlight the structure of particular atoms and functional groups through deuterium labelling. The corresponding reflectivity profiles are shown in Fig. 2: (red) chain-deuterated DMPC hydrated with H 2 O to emphasise the scattering contribution of the hydrophobic membrane core, (blue) DMPC hydrated with heavy water (D 2 O) to highlight contribution of the hydration water, and (green) DMPC and H 2 O to emphasise the head group in the bilayers.
The equally spaced Bragg peaks are indicative of a well organized lamellar membrane structure. While up to 12 diffraction orders can be observed in partially-hydrated DMPC membranes, [43][44][45][46] thermal fluctuations in fully hydrated fluid lipid membranes usually strongly suppress higher diffraction orders. 47 Four to five diffraction orders were observed in each diffractrogram in Fig. 2 and used for Fourier synthesis of the scattering length density.
The OFFSPEC measurements covered a q z range from 0.04 to 1.10 Å À1 , however Bragg peaks could be observed only to a maximum of q z B 0.50 Å À1 . The out-of-plane and in-plane component of the total scattering vector, Q, are denoted q z and q xy , respectively, in what follows. Fitting of Lorentzian peak profiles was used to determine the location of each reflection. The peak areas were integrated by manual determination of the corresponding peak bounds. An example of such an area is marked in Fig. 2 in the third reflection of the DMPC-d54/H 2 O sample (red).
The neutron scattering length density (NSLD), r(z), can be approximated by a 1-dimensional Fourier analysis, where the different Fourier components are observed in the experiment as the integrated intensities of the reflectivity Bragg peaks: ffiffiffiffiffiffiffiffi ffi I n q n p n n cos 2pnz d z

;
(1) where N is the highest order of the Bragg peaks observed in the experiment. This method is equivalent to determining electron densities from X-ray scattering. 24,45,48 The square root of the integrated peak intensities I n , multiplied by q n give the form factors, F(q z,n ). 49,50 The bilayer component form factor F(q z ), which is in general a complex quantity, is real-valued in the case of centro-symmetry. The phase problem of crystallography, therefore, simplifies to the sign problem F(q z,n ) = n n |F(q z,n )|, where the phases, n n , can only take values AE1. These phases are used to reconstruct the neutron scattering length density profiles from the diffraction pattern following eqn (1). When the form factor F(q z ) is measured at several q z values, T(q z ), which is proportional to F(q z ), can be fit: In order to determine the phases quantitatively, the form factor has to be observed at different q z -values by measuring at different contrast conditions or using the so-called swelling technique. 46 Here, fitting the experimental peak intensities and comparing them to the analytical expression for T(q z ) in eqn (2) allowed us to assess the phases, n n . Up to 5 evenly spaced diffraction peaks were observed for these curves and were used to determine the neutron scattering length density profiles for each sample. Each peak in the diffraction pattern has a corresponding phase of AE1. Phase arrays for each sample were as follows: The lamellar spacings, d z , were determined from the peak spacing in Fig. 2  As the lamellar spacing includes the thickness of the bilayer plus the thickness of the water layer between the stacked bilayer, it is an established measure of the level of hydration reached. Values of d z as a function of hydration are well documented for DMPC by the Nagle group and d z for DMPC was reported to show a power-law dependence on membrane hydration. 51 d z values in excess of 63 Å denote a significant increase in the thickness of the water layer. The different membranes prepared for this study can, therefore, be considered to be fully hydrated with more than 25 water molecules per lipid molecule such that the results of each system can directly be compared.
2.2.2 Incoherent quasi-elastic neutron scattering on Let. The incoherent quasi-elastic neutron scattering experiments were carried out at the Low-energy-transfer (Let) spectrometer 52,53 at ISIS, the pulsed neutron and muon source at the Rutherford Appleton Laboratory in Oxfordshire, UK.
In order to cover as large a range of momentum transfer, Q, and dynamical range as possible different incident neutron energies, E i , were combined: 0.64, 0.96, 1.62, 3.28, and 9.8 meV corresponding to neutron wavelengths of l = 11.3, 9.23, 7.11, 4.99 and 2.89 Å, respectively. Position-and time-sensitive 3 He detectors collected scattered neutrons with an angular range of À101 to 701 (up to 1401 for bulk water) horizontally and AE301 vertically.
A 3-dimensional spectrum (E, Q, S(Q, E)) for each incident energy was acquired in .nxspe (Nexus event-mode) format. A Matlab (Mathworks Inc.) script tailored to stacked membranes used a cylindrical geometry where the components of scattering vector, Q, are the membrane plane, q xy , and the cylindrical axis, q z .
Bulk water data were sliced using a powder geometry in DAVE (Data Analysis and Visualization Environment) software from NIST Center for Neutron Research. 54 All data were cut into ten equal slices in q z and up to 22 equal slices in q xy to make a grid of 220 squares in reciprocal space. For E i s of 0.64, 0.96, 1.62, 3.28, 9.80 meV the widths in q z and q xy were 0.027, 0.033, 0.043, 0.062, and 0.106 Å, respectively.
Each square results in a spectrum of S(q (xy,z) , E) vs. E and was fit with functions that will be described in detail in the next section (Section 2.2.3). Only data from detector positions of greater than B101 (the sixth slice and greater) in the lateral direction were used in the analysis due to contributions of the super mirror guides. The super mirror guides cause a beam divergence, which in turn affects the lower scattering angles.
Resolution measurements for the membrane system were taken using a B15 K DMPC/H 2 O sample and bulk H 2 O resolution measurements using a vanadium sample. Resolutions and quasi-elastic spectra were cut to the same dimensions in reciprocal space. The fits are based on the Matlab nonlinear least-squares minimization algorithm ''lsqnonlin''. Fits were a numerical convolution of the resolution (Gaussian function(s)) with the d-shaped function and the model for the QENS, with the addition of a background. The resolution function of Let depends on the incident energy and scattering angle and may deviate from a simple Gaussian function. Therefore, the resolution spectra were fit with one Gaussian function (E i = 0.64, 0.96, 1.62 meV) or two Gaussian functions (E i = 3.28, 9.80 meV). A Toeplitz matrix was constructed from each resolution spectrum, forming a square matrix dependent on the size of the scattering function spectrum. Convolution was carried out by multiplication of the Toeplitz matrix with the sum of the d-shaped function and broadening model.

Diffusion models.
There are different models used to describe diffusive motions. 55 For a particle diffusing via thermally-induced Brownian motion, the 1-dimensional mean squared displacement (MSD) of the particle is a function linear in time, where D is the translational diffusion coefficient. With this characteristic length scale, one can define a self time-dependent pair-correlation function for incoherent scattering, which is a solution of Fick's law, The corresponding intermediate scattering function decays exponentially, and can be Fourier transformed to give a Lorentzian-shaped incoherent scattering function Thus, quasi-elastic scattering, which exhibits itself as Lorentzian broadening with a full-width at half-maximum (FWHM) that has a Q 2 dependence, is indicative of a continuous diffusion process with E FWHM = 2 hDQ 2 . Water dynamics are more complex than that of ordinary liquids. In a simple liquid 56 atoms move ballistically at short times. The ballistic MSD h[dr(t)] 2 i p h(vt) 2 i p t 2 of the single molecule is eventually followed by a crossover to Brownian diffusion, characterized by h[dr(t)] 2 i p t, for long times. In dense or strongly interacting fluids, the caging effect leads to a plateau in h[dr(t)] 2 i. The motion in this regime is sub-diffusive with h[dr(t)] 2 i p t b with an exponent b o 1. This anomalous regime of hydration water dynamics was observed in computer simulations of lipid bilayers. 32,33 The models used to fit the incoherent neutron data incorporate a constant background, a resolution function, R(t) (in the time domain) or R(E À E 0 ) (in energy domain), which is numerically convolved with the sum of a d-like function and a function that describes the quasi-elastic broadening.
In the first model, motion of the water molecules is modelled by a single relaxation time, t: , and I N (Q) is the so-called elastic incoherent structure factor, that represents effects of confinement, i.e., correlations in the long-time limit 57,58 and S R (Q, E) is the resolution-broadened incoherent scattering function. Eqn (8) describes a single Brownian diffusion process of the water molecules with a Lorentzian quasi-elastic broadening: where B is the amplitude of the d-like function and is equal to I N (Q). A and C bkgnd are constants and G is the FWHM of the Lorentzian function. In the second model, the observed dynamics are described by two Brownian diffusion processes, resulting in two Lorentzian functions: Here, A 1 and A 2 are constants. These functions typically describe a slow motion (narrow width, G 1 ) and a fast process of the molecules (broad width, G 2 ). In the third model, we used the Kohlrausch-Williams-Watts (KWW) function to describe motion of the water molecules. 59 This function describes a distribution of relaxation processes characterized by relaxation time t and stretching exponent b. Values of b less than unity signify a sub-diffusive regime, where a stretching exponent of unity is indicative of Brownian diffusion.
Using this model one obtains resulting in The choice of which of the above functions to use to model incoherent quasi-elastic neutron scattering data is not a purely mathematical problem. The use of a single Lorentzian peak implies a relaxation process with a single relaxation time t. Two Lorentzian peaks imply two well defined dynamical processes. Stretched exponential correlation functions are often more appropriate to describe dynamics in heterogeneous systems to take into account the heterogeneous local structure in which water molecules move. b exponents of less than 1 describe a distribution of relaxation times. The dynamics are then ''glassy'' and the corresponding diffusion dynamics slower than a Brownian diffusion. This regime is referred to as ''sub-diffusive''. The difference in the distribution of relaxation times between a single and stretched exponentials, is pictured in Fig. 3 for different values of b.
A clear distinction between the different models often turns out to be difficult due to low statistics in experiments. Examples of different models fit to quasi-elastic scattering are shown in Fig. 4. While the one-Lorentzian fit shown in parts (a)-(c) fits the data reasonably well, fitting with two Lorentzians (parts (d-f)) or the Fourier transform of the KWW function (in parts (g-i)) provides a significantly better fit.
In order to quantitatively compare the different models, the reduced chi-squared values (w 2 ) of the different models were calculated using 61 where N is the number of data points summed over and n is the number of constraints in the fits. The reduced w 2 value gives a statistically valid measure of the goodness of fit: in the ideal case ofw 2 = 1, the fit is within the error bars of the spectrum. Forw 2 4 1, the spectrum is under-fitted, i.e., there is more information in the spectrum that has not been reproduced in the fit. Forw 2 o 1, the spectrum is over-fitted.
Each spectrum was fit with all three models and the corresponding reduced w 2 calculated. The results are listed in Table 1. A detailed list is given for the q x -averagedw 2 for an incident energy of 3.28 meV to show that the fit quality also depends on q z . Meanw 2 s were calculated and are given for all energies. It is observed that the fit of a single Lorentzian provides a significantly poorer fit while two Lorentzian peaks and a KWW function seem to adequately describe the quasi-elastic data with almost identicalw 2 s.
We note that this problem is not unique to incoherent quasielastic neutron scattering experiments. Lipid motion determined by fluorescence correlation spectroscopy (FCS), for instance, was reported to be equally well described by anomalous diffusion or two Brownian diffusion processes. 11 Following eqn (9), (11) and (13), the corresponding fitting functions involve 2, 4 and 3 free parameters. While a larger number of free parameters typically leads to a better fit, the fundamental question is if motion of molecules in anisotropic and heterogeneous environments is more reasonably described by a series of Brownian diffusion processes or rather by a distribution of relaxation times, as described by the KWW function. We will present evidence that the KWW provides the correct description of hydration water dynamics based on results from the neutron experiment and computer simulations further below. In the remaining sections, the quasi-elastic broadening is characterized by the Fourier transform of a KWW function to describe anomalous dynamics, and t and b parameters are determined.

2.2.4
The effect of instrumental resolution. If there is an instrumental energy resolution, R(E), then the resolution broadened dynamical structure factor, S R (Q, E) is given by By introducing the Fourier transform, R(t) of R(E), 62 S R (Q, E) can be written as, dt e iEt= h IðQ; tÞRðtÞ: When assuming a Gaussian resolution function, common in quasi-elastic incoherent neutron scattering experiments, R(E) can be written as where s is the Gaussian width of the energy resolution. This resolution has the Fourier transform When where s and E FWHM are in energy units. Molecular motions slower than the experimental observation time can still be accessed in experiments, however potential additional statistical noise in a spectrum can lead to a systematic error, which tends to overestimate the actual relaxation time, t, as will be discussed below.

Simulations
2.3.1 Atomic-detail molecular dynamics simulations of bulk water diffusion. In order to study the intrinsic character of water diffusion, the motion of water in bulk was studied using molecular dynamics at full atomic detail. We employed GROMACS to simulate 895 water molecules (SPCE model) in a NPT ensemble with the built-in Nose-Hoover thermostat (T = 300 K, t = 1 ps) and the built-in Parrinello-Rahman barostat ( p = 1 bar, compressibility 4.5 Â 10 À5 /bar, t = 1 ps).  (12) and (13). 60 b values range from 0.1 to B1 to illustrate the range of relaxation distributions for sub-diffusive systems as opposed to the d-like distribution for Brownian diffusion (approximated by b = 0.9999).
We used a time step of 1 fs and simulated for 50 ns. To ensure proper equilibration, only the trajectory beyond 1 ns was further analyzed. The start configuration is constructed with standard GROMACS routines for SPC water.
2.3.2 Coarse-grained Brownian dynamics simulations of water diffusion in a membrane system. In order to focus on the dynamical features introduced by the membrane due to different local dynamical environments and geometrical confinement, profiles for diffusivity D(z) and free energy V(z) were chosen to mimic those obtained in atomic-detail molecular dynamics simulations 32 and neutron reflectivity.
In particular, we choose: and with the bulk diffusion coefficient D bulk = 2.7 nm 2 ns À1 , the inner-membrane diffusion coefficient D in = 0.12 nm 2 ns À1 , and the parameters a = 4/nm, b = 3.6, c = 0.43/nm 2.75 . V(z) is given in units of k B T. The potential V(z) is shown in Fig. 5(b) and was modelled based on the experimentally determined distribution of water molecules in the bilayers. Note that since on sub-nanosecond time scales the crossing of water through membranes is negligible, these profiles describing only one water layer between two membranes is appropriate. The displacement Dr of the center-of-mass of one water molecule within a time step Dt is calculated via where the Wiener increment DW is drawn from a normal distribution with hDWi = 0 and hDW 2 i = 1. Note that the second summand is the apparent drift term from the Ito interpretation of the Langevin equation 63 which ensures relaxation to equilibrium at long times. The trajectory for one water molecule is simulated with Dt = 0.1 ps for a time of 10 ms. As starting condition we choose the center of the water layer, in agreement with the experimental data in Fig. 5. Intermediate scattering function and fitting procedure. From the simulation trajectories, the intermediate scattering function I(Q, t) was calculated directly for each combination of q x and q z , respectively. For the case of bulk water, I(q x , q z , t) was calculated for the center-of-mass of each molecule for an exponential time range from 0.05 ps to 30 ns, respectively, as shown in Fig. 6, and then averaged over all molecules. For the case of membrane water, an exponential time range from 0.1 ps to 300 ns was used.
The intermediate scattering functions were fitted with a stretched exponential plus a constant (non-zero only for q x = 0), consistent with the experimental analysis as laid out in eqn (13). To render the fitting more robust, we employed the following fit procedure: first, for the case of hydration water and q x = 0, I N (q x , q z ) is fixed to the mean of the last 3 points of the correlation function. For all other cases, I N (q x , q z ) is set to zero. Second, we fit a stretched exponential to the re-scaled data.
In order to investigate conditions closer to the experiment, we include a hypothetical Gaussian experimental resolution via fitting weights R(t) = exp(À 1 2 t 2 /t R 2 ) with t = 500 ps. Fig. 6 shows a set of I(Q, t) profiles with corresponding fits and artificial resolution function. Except for very small correlation times, the model of a stretched exponential provides good fits over several decades in times.

2.3.4
Mean-squared displacement and time-dependent stretch exponent. The MSDs parallel, hdx 2 i + hdy 2 i, and perpendicular, hdz 2 i, to the membrane were calculated from the center-of-mass coordinate separately. In the case of bulk water, the MSDs are averaged over all molecules. According to the MSD for an anomalous diffusion, hdr 2 i = 2(Dt) b(t) , we calculate the timedependent stretch exponent as an alternative measure of motional anomaly.

Membrane structure
In order to confirm fully hydrated bilayers neutron reflectivity measurements were done to determine the bilayer structure. Fig. 5 shows the scattering length density profiles, as determined from the reflectivity curves in Fig. 2.
The DMPC/H 2 O sample allowed us to determine the extent of the head group, the water layer and lipid tail regions. A combination of X-ray scattering, neutron scattering and simulations have been used to determine bilayer structures of a variety of phospholipids (PC). [64][65][66] Since these lipids have the same PC head group, these publications were instructive in finding the extent of our head group region. We define the extent of the head  (9)), two Lorentzians (eqn (11)) and the KWW function (eqn (13)) of each spectrum. Exemplary fits are shown in Fig. 4 Energy (meV) q z (Å À1 )  group to be the center of the distributions of its most interior/ exterior molecules: carbonyl/choline groups. The probability distributions for the head group components 65,66 show the head group/water boundary at approximately the centre of the choline distributions. The head group/water boundary should be approximately at the inflection point of the head group peak bordering the water at z = 24 Å, as shown in Fig. 5.
We assumed the maximum of our NSLD peak is centered between the phosphate and carbonyl groups. As the distribution of the carbonyl 64-66 has a half-width at half-maximum of B3 Å, we approximate the head group/tail boundary at 15.5 Å.
Half the average lamellar repeat distance minus the head group/water boundary (34.9-24.0 Å) results in half a water layer of 10.9 Å. Based on this analysis, the membrane core was found to have a thickness of 31 Å (15.5 Å per leaflet), the head group region was 8.5 Å thick and the hydration water layer had a total thickness of 21.8 Å (10.9 Å on each side of the density profile). These values were used in the simulations and confirm that the phospholipid bilayers were in their fully hydrated fluid phase. We note that the values for the bilayer dimensions are in excellent agreement with reference values for DMPC by the Nagle group 51 of D HH = 35.3 Å, D W = 19.2 Å and 9 Å for the head group thickness.
The raw NSLDs for each sample were scaled in the following way: the profile edges were scaled to the density of bulk water (r H 2 O = À5.6 Â 10 À7 Å À2 or r D 2 O = +6.39 Â 10 À6 Å À2 ) and bilayer centers were scaled to a CH 3 (or CD 3 ) group (r CH 3 = À1.83 Â 10 À6 Å À2 ) or (r CD 3 = +8.89 Â 10 À6 Å À2 ) plus the NSLD of 0.10 water molecule. The center of the bilayer in the NSLD profiles are typically scaled to a CH 3 group, however, differences between tail regions of DMPC/H 2 O and DMPC/D 2 O profiles point to non-zero water content at the center of the hydrophobic core.
We determine the water content at the center of the hydrophobic core to be 0.10 solvent molecules (H 2 O or D 2 O) per lipid such that the NSLD profiles of the two DMPC samples overlap from B6 Å to 10 Å.
Though this technique cannot determine the absolute number of water molecules in the entire tail group region, our experiments provide evidence for B1 water molecule for every 10 lipids. The presence of water in the center of lipid bilayers is supported by IR spectroscopy experiments, 36 as well as simulations. 37,39 While the density of water molecules in the hydrophobic core is low compared to the hydration water layer, the center of the bilayer is a region of low lipid density, with a larger free volume fraction 39 and consequently lower free energy and higher probability for water molecules to be present.
The NSLD profiles in Fig. 5 are in very good agreement with results from recent computer simulations in the same system by Hansen et al. 33 and Yang, Calero and Marti. 67 Both our experiment and simulations agree that water molecules penetrate the bilayer into head group and tail region, however the simulations do not provide evidence for an increased water density in the center of the bilayers. We note that the simulations used a potential to implicitly model the membrane where water was not explicitly accounted for within the membrane core.

Water dynamics from experiments
3.2.1 Bulk water. The diffusion of water molecules in two systems are discussed in this section: (1) the isotropic diffusion of bulk water (H 2 O) contained in a cylindrical annulus at room temperature, and (2) the anisotropic diffusion of hydration water (H 2 O) in fluid, oriented, deuterated lipid bilayers at T = 30 1C (303 K). In order to describe motion of water molecules in crowded or strongly interacting environments, which may induce a distribution of relaxation times, the dynamics were fit directly in the energy domain using the Fourier transform of a KWW function, where b is the exponent of the stretched exponential and t the corresponding relaxation time. Fig. 7(a) and (b) show the results for bulk water. A sketch of the scattering geometry is shown in Fig. 7(c). Stretching exponent ( Fig. 7(a)) and relaxation time (part (b)) are displayed in 2-dimensional maps as a function of parallel and perpendicular momentum transfer, q x and q z . t is plotted in units of picoseconds and on a logarithm scale due to its multiple orders of magnitude range. The map covers motions of water molecules from B3 Å to B100 Å in the x-direction and from B6 Å to B60 Å in the z-direction.
Motions over longer distances (small q x,z -values) are slower and require longer observation times (higher energy resolution or longer neutron wavelength, l) as compared to short-range motions at larger q x,z -values. Following Bragg's equation (q = 4p sin(y)/l or 2d sin(y) = l with d = 2p/q), the accessible q x and q z -range strongly depends on the incident neutron wavelength. The larger l, the larger the observation time and the higher the energy resolution, however, the smaller the q x,z -range covered. In order to cover the largest range of length and time scales, results from different energy resolutions must, therefore, be combined, as has been done in Fig. 7. The pattern of concentric rings of relaxation times in Fig. 7(b) is indicative of an isotropic motion of water molecules with longer relaxation times at longer length scales.
The values for b in Fig. 7(a) between 0.85 and 1 indicate that the corresponding nanoscale dynamics of water molecules is slightly sub-diffusive. Sub-diffusive dynamics is often found in supercooled water. Simulations of supercooled bulk water [68][69][70] show sub-diffusive behaviour with exponents in the range of b = 0.85-1 for our parameter ranges. Chen et al. 69,70 found that b stayed close to 1 for low Q values and low temperatures T o 225 K, however Sciortino et al. 68 reported that supercooled bulk water at 285 K at Q B 1.5 Å À1 reached b = 0.85. Experiments using light scattering 71 and neutron scattering 71,72 observed 1 Z b Z 0.6 for bulk water in our temperature range (293-313 K 72 ) for Q B 0.25-1.5 Å À1 .
The bulk water measurements prove that diffusion can be studied by analyzing different areas on the large 2-dimensional detector array. We note that it is not trivial to determine a diffusion constant in the sub-diffusive regime as the standard approach is to use eqn (7) and plot the quasi-elastic broadening as a function square of the scattering vector, Q 2 ; this approach fails in the case of sub-diffusion and non-Lorentzian fit functions.
3.2.2 Membrane hydration water. Fig. 8(a) and (b) depict b and t-values for the membrane hydration water as a function of q x and q z . The membranes were aligned horizontally in the spectrometer, such that q x denotes the scattering vector of the lateral motions and q z denotes motions transverse to the bilayers, as shown in Fig. 8(c). As for bulk water, five measurements with different energy resolution were combined to maximize the accessible q x,z -range. While concentric circles were observed for bulk water, the pattern of bs and ts is clearly anisotropic for hydration water in the lateral and transverse direction. b-values of B0.5 were observed at high q x and q z -values in part (a) corresponding to small distances, indicative of a sub-diffusive dynamics. b continuously increases towards longer lengths (smaller q x,z values) until Brownian diffusion with b = 1 is observed.
The relaxation times in Fig. 8(b) show fast hydration water dynamics at large q x , which continuously slows down towards longer length scales (smaller q-values). In contrast to the bulk water results in Fig. 7, the relaxation times along q x and q z are anisotropic with clear deviations from the circular pattern. The most striking feature is a localized region with slow relaxation times at q x B 0.15 Å À1 (equivalent to a length scale of 40 Å). The feature is also seen in the b data in part (a) and exhibits small values of b of less than 0.40. Shape and position of this feature agrees well with the feature observed in the computer model of Fig. 7 Values of (a) b and (b) log 10 (t) (t in ps), fit from each bulk water spectrum, are plotted as a function of q z and q x . Measurements were made with five different energy resolutions. The isotropic nature of the bulk water dynamics can be seen in the concentric rings in (b). (c) Sketch of the scattering geometry. We note that the 17th slice in each measurement was taken out because the detector array at this scattering angle was, unfortunately, malfunctioning during the experiment. the stacked bilayers and is a result of confinement, as will be discussed below. Fig. 8(d) shows q x -slices at different q z -values, which allow access to different water molecules. Based on the structural results in Fig. 5, the highest density of water molecules in the hydrated bilayer stack is found in the center of the hydration water layer. Different values of the transverse momentum transfer, q z , then correspond to the spread of z-values that the molecules explore when diffusing. q z = 0.586 Å À1 corresponds to a vertical distance of d z = 2p/0.586 = 10.7 Å, within the hydration water layer. The slice centered at q z = 0.272 Å À1 contains water molecules that penetrate the head group region of the bilayers; q z = 0.15 Å À1 molecules, which migrate to the lipid acyl chains. Based on Fig. 8(d), water molecules, which diffuse laterally in the hydration water layer show a sub-diffusive behaviour with a b exponent of B0.50. This b-value is significantly lower than the value that we measured for bulk water, indicative that hydration water motion is significantly anomalous. Water molecules in contact with lipid head groups and hydrophobic core show a sub-diffusive regime at small distances with a transition to Brownian motion at larger distances of B14 Å. We note that ballistic motion of water molecules with b-values of 2 is expected to occur at larger q x -values, outside of the experimentally accessible range in this experiment.
The results are in qualitative agreement with a recent computer simulation study by von Hansen, Gekle and Netz, 32 which reports a sub-diffusive regime for hydration water dynamics, however in the time domain, as shown in Fig. 8(e) (data adapted from ref. 32). The simulation data were evaluated for water molecules at different distances from the bilayer center, z 0 = 15 Å being located in the lipid head group region and z 0 = 35 Å, at the farthest distance from the bilayer. We note that simulations were conducted in a single hydrated bilayer, where water dynamics away from the bilayer are more likely Brownian, as compared to the confined water layer in our experiment. Water molecules diffusing along the head group region or within the lipid tails show a distinct sub-diffusive regime, with b values of 0.50, in good agreement with the experimental results in Fig. 8(d).
Corresponding cuts of the relaxation time, t, at q z -values of 0.15 Å À1 and 0.586 Å À1 (corresponding to z motions over 10.7 Å and 41.9 Å, respectively) are shown in Fig. 9. Again, slices were centered in the hydration water layer and q z -values are a measure of the corresponding slice width. The relaxation time for hydration water molecules diffusing laterally with small z motions (large q z ), inside the water layer, was found to be almost q x independent. A slight slowing down of the water molecules at a q x -value of B0.70 Å À1 is likely related to the presence of the bilayer interface as a disturbance caused by the lipid head groups. This feature corresponds to the slight decrease in b-values in Fig. 8(e) in the computer simulations by von Hansen, Gekle and Netz 32 at z 0 = 35 Å.
A hydration water molecule that moves about a range of z = 40 Å while diffusing laterally penetrates the bilayer to the lipid acyl chains. Its motion appears to slow down at a q x value corresponding to the distance between two lipid head groups indicated by the peak in the relaxation time. The relaxation times then increase significantly towards larger length scales.

Effect of experimental resolution.
A resolution effect is observed in Fig. 8(b): while the general trend is that motions become faster over smaller distances (larger Q-values), t values between the different resolution frames do not seem to fit seamlessly. When switching from a higher to a lower resolution (moving from a smaller into a larger rectangle in Fig. 8(b)), the relaxation time becomes slower again. Thus, it seems that different relaxation times can be obtained when the same length scale is measured at two different energy resolutions equivalent to different experimental observation times. Fig. 10 illustrates the effect: Fig. 10 (a) shows t distributions of the fit function to data measured at the same (q x , q z ) position in reciprocal space from the data in Fig. 8 (which should have the same relaxation time), however, measured at different experimental resolution times, t R . t R -values are given in the figure and marked by the shaded areas. For some resolutions, the peak in the time distribution is outside of the experimentally covered window and the corresponding relaxation time was extrapolated by the fit. As a consequence, when the quasi-elastic broadening was too small to be properly resolved by the instrumental resolution, the relaxation time was found to be systematically over-estimated.
Though counterintuitive, this behaviour is a systematic effect that can be modelled. If there is little or no noise or error, fit functions can typically be determined correctly even when the function is not fully covered by the experimentally accessible range of data. This is the case when the relaxation time is outside the experimental resolution window, as defined in eqn (20). The combination of experimental noise and partial coverage, however, may cause a systematic error in the determination of relaxation times in quasi-elastic incoherent neutron scattering experiments. In Fig. 10(b), the t-distribution from part (a) was generated and a window equivalent to the experimentally accessible time window was selected. Statistical noise was added to the generated data to mimic experimental conditions. The resulting data were then fitted again with the original distribution that was used to create the distribution (and should reproduce the original data). The results are shown as solid lines in part (b) and well reproduce the systematic deviation: the relaxation time tends to be overestimated when the distribution is outside of the experimental window and experimental noise is present.

Water dynamics from computer modelling
In order to consolidate the picture of anisotropy and anomaly of water dynamics obtained by the neutron scattering experiments, Fig. 9 Lateral diffusion of membrane water molecules as shown by relaxation time t (t in ps). Slices along q x are shown for two values of q z . The lower curve occurs at q z = 0.15 Å À1 corresponding to length scale of B40 Å. The lateral motion of a water molecule is slowed significantly for a lateral length scale of B10 Å denoted by an increase in t at q x = 0.6 Å À1 ; most likely due to the water's interaction with the head group, where hydrogen bonding occurs. The upper curve occurs at q z = 0.586 Å À1 corresponding to B10 Å. This signal is most likely dominated by water molecules diffusing in the hydration water layer. Dashed lines correspond to maximum viable t; lines and symbols have corresponding colors. we employed two sets of computer modelling. First, using fullatomic detail molecular dynamics, we studied the intrinsically anomalous character of isotropic, bulk water diffusion at picosecond time scales.
For hydration water dynamics, we used a coarse-grained Brownian dynamics simulation to investigate the effects of geometrical confinement and locally varying dynamical environments of water molecules interacting with lipid bilayers. The corresponding systems are inherently larger in size and the number of molecules involved, therefore, making full atomistic simulations challenging.
3.3.1 Bulk water. Fig. 11 summarizes the results for the atomistic bulk water simulations. b-values are shown in part (a) and the corresponding relaxation times in part (b). As expected for bulk water in isotropic and unconfined conditions, no angular dependence is observed in the relaxation time t and the stretch exponent b. Patterns of concentric circles are observed, in agreement with the experimental data. t increases monotonically with decreasing q x or q z , i.e., increasing length scale. In agreement with the experimental results bulk water dynamics on picosecond time scales and over Ångström distances (large q x,z -values) is sub-diffusive, with b exponents of B0.75. The diffusive nature changes to Brownian at low q x,z .
The MSD and the derived time-dependent stretch exponent b are shown in Fig. 11(c) and (d). The displacements along different directions (hdx 2 i, hdy 2 i and hdz 2 i) coincide, indicative of an isotropic motion. After a ballistic regime at sub-picosecond time scales, a sub-diffusive regime at time scales of several picoseconds follows, and, finally, normal diffusion is reached at nanosecond time scales, in agreement with computer simulations by von Hansen et al. 32 Thus, both experiments and simulations consistently indicate that water diffusion has an intrinsically subdiffusive character on picosecond time scales, also under ambient conditions.
3.3.2 Membrane hydration water. Fig. 12 shows the results for the simulations of water in proximity to membranes.
Before discussing the 2-dimensional maps, we stress that the computational framework used here involved Brownian dynamics, and can, therefore, not account for the intrinsic sub-diffusive character of water diffusion in the bulk. Consistently, the meansquared displacement and the time-dependent stretch exponent b should indicate normal diffusion at small times (with b = 1).
However, the stretch exponent b in Fig. 12(a) shows subdiffusive behaviour at all q x values, with a trend of stronger sub-diffusive motion at higher q x . Pronounced anisotropic behaviour is observed at low q x . The constant I N (Q) is zero except for q x = 0 Å À1 and q z t 0.2 Å À1 , indicating geometrical confinement due to the membrane on a length scale of B30 Å. The relaxation time t (plotted in units of 1 ps in Fig. 12(b)) shows anisotropic behaviour at low q x . Water motion is overall slower than in bulk water in Fig. 11(b).
The significantly longer t results from the slowing down of dynamics close to the membrane, and a longer average relaxation time. The value of b indicates a distribution of local dynamical environments, and thus a distribution of relaxation times. This finding is consistent with the observed trend of b approaching 1 with lowering q x,z , since the molecule explores a larger space and thus performs a temporal average of dynamical environments, finally reaching normal diffusion with the average diffusion coefficient.
With increasing correlation time, the molecules experience confinement due to the membrane, which induces an apparent sub-diffusive character for the diffusion displacement. We note that this factor is only caused by the z coordinate, as shown in the MSD in Fig. 12(c), while diffusion along x and y is normal diffusion, as expected from the simulation framework.
At low q x,z , distinct anisotropic features appear in t, b and I N , mirroring the features observed in experiments. In particular, when cutting along the q x axis for small q z , a non-monotonous behaviour is observed for b and t. The evolution of b as a function of time is plotted in Fig. 12(d). This non-monotonicity seems to be directly related to confinement effects, since its where the second line is the first term of the so-called cumulant expansion, the so-called Gaussian approximation. For a fluctuation with wave vector Q to relax, the MSD has to increase sufficiently. For the case of the membrane system, the z motion is confined on the experimental time scales, i.e. hdz 2 i settles at L layer B 40 Å for long times, in perfect agreement with the experiments. Consequently, for q x = q y = 0, correlation remains for long times, which is represented by I N (Q) from eqn (8). At large q x , the remaining long-time correlations are relaxed only through the x motion, i.e., diffusion in z and x direction play a fundamentally different role. Small q z corresponds to a larger remaining correlation in the z direction, which in turn means that a longer time is needed to relax the fluctuation via the other coordinates. We note that the argument is not dependent on the Gaussian approximation, but similarly holds for higher cumulants, since all z cumulants are constant and confined in the long time limit.

Motion of water molecules
Lipid diffusion has been studied for decades and it is commonly accepted that the Brownian motion of lipid molecules over long, micrometer length scales is characterized by a continuous diffusion process. [13][14][15][16][17]21 Although this is a well studied fundamental process, our understanding of nanoscale diffusion in membranes is still being challenged by new results from experiments and simulations.
Experimental evidence that lipids move coherently in loosely bound clusters, rather than as independent molecules was presented. 7,18-20 A ''hopping'' diffusion of lipids into nearest neighbour sites was observed in single supported bilayers. 6 Ballistic lipid motion in fluid membranes was recently reported. 5 It has also been suggested that there is a flow-like component to the motion of the lipid molecules over long length scales. 7 Significantly less work has been done to study dynamics of hydration water molecules, most likely due to the experimental and computational challenges. We find that lateral water diffusion is affected at all distances from the bilayers and motions of hydration water molecules differs significantly from the isotropic diffusion of bulk water.
3.4.1 Anomalous dynamics of hydration water. By labelling water molecules, lipid head groups, and lipid tails through selective deuteration, dimensions of hydration water layer, bilayer head group and tail region were determined by neutron reflectivity. While the highest water density is found in the center of the hydration water layer, the NSLDs give evidence for water molecules that are able to penetrate the hydrophobic membrane core; B1 water molecule per 10 lipid molecules were found in the center of the fully hydrated membranes.
From the quasi-elastic neutron scattering experiments and computer simulations, we obtain a consistent picture of the phenomenology of water diffusion in membrane systems, as well as different origins of anomaly and anisotropy. As a first origin of anomaly, water shows an intrinsic anomalous diffusion on picosecond time scales, as also observable in bulk water. The second origin of anomalous diffusion is related to the distribution of local dynamical environments of hydration water molecules, which leads to a distribution of relaxation times and sub-diffusive behaviour. The anisotropic effects of confinement and local dynamical environments on the motions in different directions give rise to directiondependent relaxation times and stretch exponents due to the varying sampling time of the system along the unconfined coordinates. Fig. 12 Results from coarse-grained Brownian dynamics simulation for water confined in a membrane system: (a) anomaly of water diffusion due to dynamical heterogeneity and confinement: the stretch exponent b shows sub-diffusive behaviour at all visible q x values, with a trend of stronger subdiffusive motion at higher q x . Anisotropic behaviour is observed at low q x due to geometrical confinement of the membranes. (b) The relaxation time t (plotted in units of 1 ps) shows a slight anisotropic behaviour at low q x . Water diffusion is overall slower than in bulk (see Fig. 11(a)). (c) and (d) Meansquared displacements hdr 2 i and time-dependent stretch exponents b(t) (see eqn (24)) for different coordinate directions evidence a confined motion in z-direction (blue dashed), while the MSD in xy-direction (red dash-dotted) resembles that of normal diffusion. The resulting overall MSD (black solid) has an apparent sub-diffusive character emerging solely from the confinement of hydration water molecules.
3.4.2 Anomalous diffusion models. Different models are typically used in the literature to describe membrane and hydration water dynamics in quasi-elastic neutron scattering experiments. Through a comparison between a single Lorentzian, two Lorentzians and the KWW function to fit the quasi-elastic broadening, the single Lorentzian was found to give the poorest agreement, while two Lorentzians and the KWW function both provided a satisfying fit. The fundamental difference between the two models is the description of relaxation times in the time domain. The two Lorentzian peaks typically correspond to a slow and a fast process, each with a single relaxation time. While the slow process is well captured, the observed fast process typically partially lies outside of the instrumentally accessible energy window and basically leads to a broad, curved background. This procedure is often used in the literature. While the KWW function fits 3 parameters, the model using two Lorentzian functions has 4 free parameters. Fitting more free parameters typically provides better fits with smaller w 2 -values.
A clear distinction of which model to use is often difficult in QENS data because of the limited statistics in the data, such that different models fit the data equally well. We note that FCS suffers from the same issue and cannot unambiguously determine the basic character of lipid dynamics in membranes. 11 However, we argue that motion of water molecules in anisotropic and heterogeneous environments is physically better described by a distribution of relaxation times (a KWW function) rather than a sum of individual Brownian processes (several Lorentzian functions). This interpretation is supported by the computer simulations, i.e., by the linear scaling of intermediate scattering function in the log-log representation over several decades which is indicative of a single stretched exponential.
We further tested the hypothesis of Lorentzian vs. KWW function in the experiment by using different instrumental resolutions. If the processes fitted by the 2 Lorentzian peak broadenings were actual physical processes, their corresponding width should be constant and they should be independent of the instrumentally available energy window. However, each energy resolution was well fit by a narrow and a broad Lorentzian peak, with no correlation between the broad peaks at different resolutions. This is strong evidence that the broad Lorentzian is likely a fitting artefact in the case of hydration water. We can, of course, not make a statement about lipid dynamics as this was not part of the current study.
Measurements at different hydrations could possibly be used to distinguish between Lorentzian and KWW models in QENS experiments in the future. One may expect that de-hydration leads to a more heterogeneous environment of water molecules and correspondingly smaller values of b, possibly up to the point where Lorentzians no longer describe the data.

Conclusion
Diffusion of hydration water molecules in hydrated phospholipid bilayers was studied using 2-dimensional incoherent quasielastic neutron scattering combined with computer modelling.
By aligning the membranes horizontally in the spectrometer, lateral and transverse molecular motions can be studied separately, but simultaneously, by analyzing different positions on the 2-dimensional detector array. The dynamics were described in the energy domain by KWW functions, corresponding to stretched exponentials in the time domain.
The nanoscale dynamics of hydration water molecules is anisotropic, and shows a sub-diffusive behaviour on nanometer length scales. As a first origin of the anomalous character of hydration water we have identified the intrinsic sub-diffusion of water at picosecond time scales, as also observed in bulk water. Secondly, a dynamically heterogeneous environment causes a distribution of relaxation times that finally appear as a subdiffusive motion. These 2-dimensional data show a distinct regime of very slow dynamics at small Q-values, related to the geometry of the stacked bilayers used for the experiment.
Diffusion experiments using 2-dimensional detectors will in the future be used to precisely measure penetration and permeability of water molecules across bilayers as a function of bilayer composition, 40,73 and to study the interaction of water with membrane inclusions, such as proteins and peptides.