Time-domain proton-detected local-field NMR for molecular structure determination in complex lipid membranes

Proton-detected local-field (PDLF) NMR spectroscopy, using magic-angle spinning and dipolar recoupling, is presently the most powerful experimental technique for obtaining atomistic structural information from small molecules undergoing anisotropic motion. Common examples include peptides, drugs, or lipids in model membranes and molecules that form liquid crystals. The measurements on complex systems are however compromised by the larger number of transients required. Retaining sufficient spectral quality in the direct dimension requires that the indirect time-domain modulation becomes too short for yielding dipolar splittings in the frequency domain. In such cases, the dipolar couplings can be obtained by fitting the experimental data; however ideal models often fail to fit PDLF data properly due to effects of radiofrequency field (RF) spatial inhomogeneity. Here, we demonstrate that by accounting for RF spatial inhomogeneity in the modeling of R-symmetry-based PDLF NMR experiments, the fitting accuracy is improved, facilitating the analysis of the experimental data. In comparison to the analysis of dipolar splittings without any fitting procedure, the accurate modeling of PDLF measurements makes possible three important improvements: the use of shorter experiments that enable the investigation of samples with a higher level of complexity, the measurement of C–H bond order parameters with smaller magnitudes $\left|S_{\mathrm{CH}}\right|$ and of smaller variations of $\left|S_{\mathrm{CH}}\right|$ caused by perturbations of the system, and the determination of $\left|S_{\mathrm{CH}}\right|$ values with small differences from distinct sites having the same chemical shift. The increase in fitting accuracy is demonstrated by comparison with ²H NMR quadrupolar echo experiments on mixtures of deuterated and non-deuterated dimyristoylphosphatidylcholine (DMPC) and with 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) membranes. Accurate modeling of PDLF NMR experiments is highly useful for investigating complex membrane systems. This is exemplified by application of the proposed fitting procedure for the characterization of membranes composed of a brain lipid extract with many distinct lipid types.

1 Introduction

The methodology for characterizing molecular structure in biological systems has been advancing rapidly over recent years (Cheng, 2018; Gauto et al., 2019; Wu and Lander, 2020). Among the various experimental techniques available, solid-state NMR spectroscopy provides the most powerful methods for investigating the molecular structure of smaller molecules (<30 kDa) undergoing anisotropic motion in membranes such as lipids, peptides, and membrane proteins (Andersson et al., 2017; Löser et al., 2018; Bacle et al., 2021; Zerweck et al., 2017; Aisenbrey et al., 2019; Cady et al., 2007; Park et al., 2012; Mandala et al., 2020; Umegawa et al., 2018; Mandala et al., 2018).

Of these techniques, researchers have consistently used ²H NMR spectroscopy since the 1970s for the molecular structural characterization of both lipids and peptides (Seelig, 1977; Davis, 1983; Strandberg and Ulrich, 2004; Leftin and Brown, 2011). The popularity of this method is due to both its simplicity and the high accuracy of the experiments. However, it requires specific ²H isotopic labeling for site resolution. This requirement severely limits its application to the investigation of biological extracts or complex model membranes that mimic biological systems. Alternatives to ²H NMR are separated local-field (SLF) 2D NMR experiments that make use of ¹H–¹³C heteronuclear dipolar recoupling during the indirect time (Hester et al., 1976). This technique is very advantageous compared to ²H NMR because SLF delivers the same type of information with ¹³C chemical shift selectivity and does not require isotopic labeling. These SLF ¹H–¹³C heteronuclear dipolar recoupling NMR techniques fall into two groups, carbon-detected local-field (CDLF) and proton-detected local-field (PDLF) experiments, with the latter having enhanced resolution (Nakai and Terao, 1992; Schmidt-Rohr et al., 1994; Bärenwald et al., 2016). Here, we discuss only PDLF techniques which provide C–H bond site resolution and are therefore directly comparable to ²H NMR.

A number of PDLF pulse sequences have been implemented to investigate lipid membranes and other systems (Griffin, 1998; De Paëpe, 2012; Molugu et al., 2017). The main difference between the distinct PDLF sequences reported are the heteronuclear dipolar recoupling blocks used to counteract the effect of magic-angle spinning (MAS) on the anisotropic Hamiltonian terms during the indirect dimension. Several recoupling pulse sequences have been designed for this purpose such as REDOR (Gullion, 1998), DROSS (Gross et al., 1997; Leftin et al., 2014), R-symmetry-based pulse sequences (Levitt, 2007; Dvinskikh et al., 2004; Hou et al., 2011), and others (De Paëpe, 2012). In recent years, an R-symmetry-based PDLF sequence designed by Dvinskikh and coworkers (Dvinskikh et al., 2004), dubbed R-PDLF (incorporating R18${}_{\mathrm{1}}^{\mathrm{7}}$ recoupling blocks), has been frequently used to successfully investigate lipid membranes (Dvinskikh et al., 2005; Ferreira et al., 2013; Löser et al., 2018; Bacle et al., 2021; Fridolf et al., 2022). The R-symmetry sequences are particularly advantageous, since they enable one to simultaneously recouple the heteronuclear and decouple the homonuclear dipolar interactions. Moreover, R-symmetry recoupling can be applied directly to an IS spin ensemble in thermal equilibrium with no need of generating transverse magnetization before the recoupling blocks and therefore avoiding potential T₂ losses during the recoupling period, i.e., increasing sensitivity. It is well-known that, as with other dipolar recoupling pulse sequences, R-symmetry recoupling is rather sensitive to radiofrequency field (RF) imperfections (Nishimura et al., 2001; Schanda et al., 2011) and therefore to the RF spatial inhomogeneity across the sample investigated – one unavoidable and unique feature of the MAS probe used to perform experiments (Tošner et al., 2017). To minimize this problem, windowed R-symmetry sequences have been proposed that are not as sensitive to RF inhomogeneity as the original windowless R-symmetry sequences (Gansmüller et al., 2013; Lu et al., 2016). These windowed versions, however, require using higher RFs which is a practical bottleneck due to hardware limitations.

In general, the strategy for developing better dipolar recoupling methods has been to minimize the sensitivity to RF spatial inhomogeneity, by using/designing dipolar recoupling sequences that are less sensitive to RF imperfections and reducing the sample studied to a narrow volume, and using recoupling pulse sequences with a high scaling factor that enable the measurement of smaller dipolar couplings (Dvinskikh et al., 2005; Chevelkov et al., 2009; Schanda et al., 2011; Gansmüller et al., 2013; Lu et al., 2016; Asami and Reif, 2017; Jain et al., 2019; Nimerovsky and Soutar, 2020). This last point is of significant importance, as higher scaling factors make it possible to obtain dipolar splittings in the indirect frequency domain with a reduced maximum time in the indirect dimension.

PDLF experiments have been mostly analyzed by “reading off” the dipolar splittings obtained in the indirect frequency domain of the 2D spectra. In contrast to the standard analysis of REDOR or DIPSHIFT dipolar evolutions (Gullion, 1998; deAzevedo et al., 2008), PDLF measurements have been rarely analyzed by employing fitting procedures. The few time-domain (TD) analysis attempts reported assumed samples under a homogeneous RF field (Gross et al., 1997; Gansmüller et al., 2013; Fridolf et al., 2022); however, it is well-known that dipolar recoupling sequences are highly sensitive to RF inhomogeneity (Schanda et al., 2011).

Here, we present an analysis procedure that enables one to increase the accuracy and applicability of R-symmetry-based PDLF NMR (R-PDLF) through R-PDLF numerical simulations that take into account the RF spatial inhomogeneity of the probe used explicitly, i.e., rather than trying to minimize the effect of RF spatial inhomogeneity, we include it in the NMR simulations used to fit the experimental measurements. The use of RF inhomogeneity for data analysis has been done recently for investigating protein molecular structure (Xue et al., 2022). Here we apply this strategy to analyze R-PDLF experiments that were proven suitable to investigate molecular structure in lipid membranes. The novelty of the procedure presented is solely in the analysis of the experimental data, no special hardware or different setup of the experiment is required. We tested the methodology on a DMPC/DMPCd sample to allow comparison of the ¹H-¹³C dipolar couplings determined by a R${\mathrm{18}}_{\mathrm{1}}^{\mathrm{7}}$-type PDLF sequence (Dvinskikh et al., 2004, 2005) with the C–D bond order parameters determined from ²H quadrupolar couplings. We then implemented such an analysis procedure to investigate a system of complex lipid membranes composed of a brain lipid extract. This enabled us to measure several dipolar couplings from the membrane hydrophilic layer composed of many types of distinct phospholipid headgroups. We believe that this strategy is highly useful for investigating both model systems and complex membrane systems with a low signal-to-noise ratio such as biological lipid extracts. In addition it has the potential for investigating in-cell molecular structures.

2 Methodological framework

The typical approach for determining C–H bond order parameter magnitudes with PDLF dipolar recoupling NMR is by performing a Fourier transform,

$$\begin{array}{}\text{(1)}& I\left({\mathit{\omega }}_{\mathrm{1}}\right)=\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}s\left(t_{\mathrm{1}}\right)\exp (-i{\mathit{\omega }}_{\mathrm{1}}{t}_{\mathrm{1}})\mathrm{d}{t}_{\mathrm{1}},\end{array}$$

of the time-domain recoupled dipolar modulation over the indirect dimension, s(t₁), and reading-off the splitting(s) of the resulting spectral dipolar line shape (Dvinskikh et al., 2004; Gross et al., 1997). If the scaling factor, κ, for the particular pulse sequence under use is known, determining the magnitude of a C–H bond order parameter from the splitting(s) observed is straightforward by using

$$\begin{array}{}\text{(2)}& \left|S_{\mathrm{CH}}\right|={\mathit{\kappa }}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\nu }}{d_{\mathrm{CH}}}}=|{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle \mathrm{3}{\cos }^{\mathrm{2}}\mathit{\theta }-\mathrm{1}\rangle |,\end{array}$$

where d_CH is the magnitude of the rigid dipolar coupling for a static C–H bond equal to ≈ 22 kHz, the angle θ is between the C–H bond and the symmetry axis of its uniaxial motion (in case of lipid membranes, this is the bilayer normal), and the angular brackets denote a time average over a time interval up to approximately the inverse of the splitting measured.

For a dipolar modulation corresponding to a single C–H bond order parameter and acquired under optimal experimental conditions, the use of Eq. (2) allows accurate values to be obtained provided that the digitized signal, s(t₁), extends to a time longer than the reciprocal of the scaled dipolar coupling frequency (e.g., four times the reciprocal gives approximately 10 % error). If s(t₁) is instead a superposition of two distinct contributions with different S_CH, for resolving the two components in the spectrum, the signal acquisition must be, at least, longer than half of a beat period defined by the frequencies of the crystallite orientations responsible for the two splittings (Lindon and Ferrige, 1980). This is illustrated in Fig. 1 for a superposition of two Pake patterns. Therefore, both for single-component dipolar modulations with low order parameters and for two-component dipolar modulations with a small frequency difference between the two components, a large number of indirect dimension points is needed in 2D dipolar recoupling experiments which severely increases the experimental time. Such a requisite is especially difficult to achieve for samples that require acquisition of a large number of transients in the direct dimension. Moreover, there is a maximum limit to the number of indirect dimension points acquired both due to RF heating of samples and to maintain the integrity of the NMR probe.

Figure 1 The effect of using a limited number of time-domain points to describe a superposition of two Pake patterns with distinct splittings.

Figure 2 Numerical simulations of the R-PDLF pulse sequence for a fixed dipolar coupling displaying the effects of RF nutation frequency deviation from ideal settings (a) and ¹H chemical shift offset (b) given in parts per million units. The solid red and faint red curves correspond to the dipolar modulation with ideal settings which yield a dipolar splitting equal to 0.315S_CHd_CH. Negative and positive chemical offset shifts with the same magnitude have exactly the same effect. All simulations in this work were performed for a ¹H Larmor frequency of 400 MHz.

Analyzing the measured data with a fitting model (in either the time or frequency domains) can be used to circumvent the limitations outlined above. The practical bottleneck for such an analysis is the RF inhomogeneity across the sample, intrinsic to the majority of experimental setups (Tošner et al., 2017), combined with the sensitivity of recoupling pulse sequences to RF strength accuracy and the effect of chemical shift offsets (Schanda et al., 2011; Lu et al., 2016). Figure 2 illustrates such dependencies for the R-PDLF pulse sequence (Dvinskikh et al., 2004), showing a non-linear dependence of the dipolar modulation on the RF pulse power level deviation from the ideal value and on the ¹H chemical shift offset. The effect of RF miscalibration becomes rather pronounced at RF frequencies higher or lower than 6 % of the ideal frequency. This means that in typical experiments, for which the RF spatial inhomogeneity profile may reach 80 % of the maximum value in the outer parts of the sample (Tošner et al., 2017), one should expect deviations from the ideal dipolar lineshape. The effect of the ¹H chemical shift offset is almost negligible up to 4 ppm above which large deviations from the ideal behavior occur. For achieving optimal resolution in an R-PDLF experiment, the chemical shift offset should therefore be no more than 4 ppm at a magnetic field inducing a Larmor frequency for ¹H equal to 400 MHz.

The outcome of an experiment is always the sum over all the detectable spatially distributed sample volumes. Therefore, the experimental data measured with an R-PDLF experiment are a sum over dipolar modulations, each modulation with a characteristic lineshape that depends on the local RF field. To simulate realistic R-PDLF experiments accounting for RF inhomogeneity, we first measured the RF inhomogeneity in the probe used in this work by the method suggested by Odedra and Wimperis (2013) (details given in the Supplement). The simulation data were then generated by integrating over R-PDLF NMR simulations that covered the RF spatial inhomogeneity measured. This was done by composing a MATLAB script with SIMPSON call-outs enabling to build an R-PDLF simulation database with realistic modulations for a range of S_CH values. The database generated included the effect of chemical shift offsets and RF miscalibration with or without a built-in RF spatial inhomogeneity profile. The MATLAB/SIMPSON files and the NMR database generated are available as open data (https://github.com/tfmFerreira/inhomogeneous-rf-nmr.git).

Figure 3 Simulation of the effect of RF spatial inhomogeneity on R-PDLF dipolar modulation and illustration of the advantage of time-domain analysis of the data over reading-off the dipolar splitting(s) in the frequency domain. (a) R-PDLF pulse sequence used in this work (Dvinskikh et al., 2004). (b) Ideal R-PDLF dipolar modulation with two components having C−H bond order parameters equal to 0.03 and 0.04 (gray lines display the individual components used). In the frequency domain, the two splittings cannot be distinguished, due to the limited number of points used in the time domain. (c) The same as in (b) but including the effect of RF spatial inhomogeneity taken as a Gaussian distribution that reflects the RF inhomogeneity of the coil used in this work (details in the SI). (d) Simulation of an experimental result with low signal-to-noise ratio by addition of random noise to the simulated curve in (c). On the bottom plot, the curves show the result of fitting the data with two components using numerical simulations performed with ideal settings (red) and by accounting for the RF spatial inhomogeneity profile (blue).

The result of accounting for RF inhomogeneity in the R-PDLF simulations is exemplified in Fig. 3 for a two-component dipolar modulation, consisting of two distinct C–H bond order parameters of 0.03 and 0.04. The main effect of including RF spatial inhomogeneity in the simulation is the damping of the signal and the non-zero long time average of the modulation which gives rise to a middle peak in the dipolar spectrum as also reported previously by Polenova and coworkers (Lu et al., 2016). Due to the limited number of points used, the distinct dipolar couplings are not visible in either dipolar spectrum of the two simulations. Figure 3d shows the result of adding random noise to the simulated data in Fig. 3c and fitting the data either with an ideal model or by taking the RF spatial inhomogeneity into account. The fit performed with NMR simulations having ideal settings gave values of 0.034 and 0.094, thus precluding the use of such fitting procedure for accurate analysis under such conditions. In contrast, the fit performed based on NMR simulations accounting for RF spatial inhomogeneity matches the data almost perfectly – as expected since these simulations were used to generate the original data – giving order parameters values of 0.028 and 0.040, very close to the original values of 0.03 and 0.04 used. This suggests that, if the (coil dependent) RF spatial inhomogeneity across the sample is measured with reasonable accuracy, it would be possible to simulate a set of realistic data for the complete range of dipolar couplings (0–22 kHz for ¹H-¹³C). This would better describe the experimental measurements and therefore enable a more accurate analysis.

3 Results and discussion

3.1 R-PDLF fits including RF inhomogeneity enable higher accuracy and shorter experiments

From the description shown in the preceding section, it is clear that RF spatial inhomogeneity affects the dipolar modulation in R-PDLF experiments. Therefore, for accurate fits of experimental data, this effect needs to be taken into account. In this and the following section we demonstrate that accounting for the effect of RF inhomogeneity in time-domain fits of experimental data indeed leads to highly accurate fits of the experimental data and consequently to a considerable improvement of the accuracy of the C–H bond order parameters determined. Instead of in the time domain, the fitting could also be performed in the frequency domain, i.e., by performing the Fourier transform of the experimental data in the indirect dimension and by then fitting these data with the Fourier transform of the numerical simulations. Here, we simply use time-domain fits and avoid performing the additional step of the Fourier transform which is unnecessary and has an extra computational time cost.

Figure 4 shows the application of this time-domain analysis to a sample of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) multilamellar vesicles (MLVs), illustrating the analysis method on the β carbon of the POPE headgroup. POPE is one of the most abundant lipids in cellular membranes. It has recently been investigated with R-PDLF NMR in the context of the NMRlipids project. This project aims to resolve the atomistic molecular structure of the most abundant lipids in nature. To this end, the project uses a combination of solid-state NMR experiments and all-atom molecular dynamics simulations (Bacle et al., 2021). Using only the first eight points along the indirect dimension, the time-domain curve fitted to the experimental dipolar modulation measured for the β carbons gives an estimation of S_CH=0.028${}_{\mathrm{0.009}}^{\mathrm{0.034}}$, where the superscript and subscript are upper and lower limits as defined in the Supplement. Such a fit enables us to predict the experimental dipolar spectrum that one would obtain using a higher number of indirect dimension points as shown in Fig. 4b and c. Figure 4d shows how the time-domain fit procedure enables the estimation of an accurate order parameter using much shorter experiments in comparison to the use of the typical Fourier transform methodology (gray circles). As it will be shown in a later section, this is particularly useful for complex systems for which only a limited number of indirect dimension points can be measured. Moreover, the fit procedure enables the estimation of the uncertainty of the value measured (the detailed description for defining the errors is given in the Supplement), which is more difficult to define with the Fourier transform “read-off” method alone. Figure 4d also shows the result of using fits with the ideal model (i.e., that does not include the effect of RF inhomogeneity). One observes that these fits become less accurate with an increase of the number of points used (a comparison with the model including RF inhomogeneity is also given in the Supplement in Fig. S2). This is a consequence of the effect of RF inhomogeneity which becomes more prominent at longer t₁ times. However, Fig. 4d also suggests that fitting with the ideal model may provide good results when analyzing short t₁ intervals. Based on previous ²H NMR measurements, the β methylene carbon is expected to have the equivalent order parameters for the two C−H bonds (Gally et al., 1975). In the next section we demonstrate the usefulness of R-PDLF time-domain fits to analyze cases of two-component dipolar modulations and estimate the accuracy of the order parameters determined with the proposed method by comparison with ²H NMR quadrupolar splittings.

Figure 4 The proposed method applied to a sample of POPE MLVs. (a) ¹³C rINEPT spectrum of the POPE MLVs. (b) Time-domain (TD) fit accounting for RF inhomogeneity of the dipolar dimension of an R-PDLF experiment with a total of eight points (thick solid line). The experimental dipolar spectrum is shown on the right together with the fit prediction of the experimental dipolar spectrum (using a total of 32 time-domain points). The thin lines show the time-domain curves that correspond to the error limits as defined in the Supplement. (c) The same as in (b) but fitting a total of 32 experimental points in the indirect dimension. A single-exponential decay of ∼ 14 Hz was multiplied to the simulated data. (d) Order parameter estimated by the time-domain fit including RF inhomogeneity (blue), with the ideal model (red), and calculated from the dipolar splitting in the Fourier transform of the experimental data (gray) as function of the total number of points in the indirect dimension used.

Figure 5 The proposed methodology applied to a sample of a DMPC/DMPCd54 liquid crystalline system. (a) ¹³C rINEPT spectrum and labels used to identify different carbons. (b) Time-domain fits neglecting the RF inhomogeneity across the sample. (c) Time-domain fits accounting for the RF inhomogeneity across the sample using one single dipolar coupling as fit parameter and neglecting any possible relaxation during t₁. (d) The same as in (c) but using two dipolar couplings as fit parameters assuming a two-component dipolar modulation.

Figure 6 ²H NMR spectrum of the DMPC/DMPCd54 liquid crystalline system and symbols showing the positions used to determine order parameters.

3.2 Comparison with ²H NMR quadrupolar splittings

To further test the accuracy of using a model of R-PDLF experiments that includes RF inhomogeneity, we used a water/DMPC/DMPCd54 liquid crystalline system (L_α phase), such that it was possible to measure the ²H NMR spectrum of the perdeuterated acyl chains of the DMPCd54 molecules and compare it with the ¹H-¹³C dipolar couplings determined from the acyl chains of the molecules with natural abundance of isotopes. Figure 5 shows the R-PDLF dipolar modulation for a number of resolved carbons from the liquid crystalline system and their corresponding fits with and without accounting for RF inhomogeneity. The improvement of the fits by accounting for RF inhomogeneity is visually evident, especially for the smaller coupling. We stress here that the RF inhomogeneity profile used in the simulations was experimentally measured and modeled by a Gaussian profile. The RF inhomogeneity is therefore an objectively quantified a priori constraint that we have obtained from an independent experiment. Moreover, since each of the three carbon peaks analyzed is a sum of components from the two individual acyl chains, sn-1 and sn-2, we also performed the time-domain fits of the dipolar modulations using two components. In the case of the acyl chain methyl groups ω, equivalent dipolar couplings minimize χ². On the other hand, for the ω−1 and ω−2 methylenes, two distinct dipolar couplings enable a better fit than a single value which suggests a distinct motional geometry in the sn-1 and sn-2 acyl chains for these methylenes (a detailed analysis is given in the Supplement). The non-equivalence of these order parameters is in agreement with previous investigations of specifically deuterated dipalmitoylphosphatidylcholine (DPPC) molecules (Seelig and Niederberger, 1974).

Table 1 Comparison of the order parameter magnitudes, $\left|S_{\mathrm{CH}}\right|$, determined from the DMPC/DMPCd54 liquid crystalline sample for the ω, ω−1, and ω−2 carbons from using the observed splitting in the frequency domain (FT) of an R-PDLF experiment with a total of 32 points in the indirect dimension, from time-domain (TD) fits with numerical simulations either accounting or not for the effect of RF inhomogeneity on R-PDLF experiments and using a single or two order parameters, and from quadrupolar splittings observed with a ²H quadrupolar echo.

The ²H NMR spectrum measured from DMPCd54 molecules is shown in Fig. 6. The five lowest quadrupolar splittings were used to calculate C–D bond order parameters. The results are included in Table 1 together with the results from the time-domain fits of the dipolar modulations from the molecules with natural abundance of isotopes. The agreement between the values obtained by using the different techniques is striking. The maximum difference between the order parameter magnitudes determined with the methodology proposed and the values determined from the ²H NMR spectrum is below ±0.005, much below the errors reported previously (up to ±0.02) when doing similar comparisons between dipolar recoupling experiments with ²H NMR spectroscopy (Gross et al., 1997; Dvinskikh et al., 2005; Ferreira et al., 2013). Here, we used values taken from previous studies for the rigid ¹H-¹³C dipolar and ²H quadrupolar couplings, namely 22 and 126 kHz, respectively (Dvinskikh et al., 2004; Davis et al., 2009). We may therefore exclude a significant isotope effect on the molecular dynamics of the acyl chains in the water/DMPC/DMPCd54 liquid crystalline system. Moreover, we can unambiguously assign the C–D order parameters measured with ²H NMR based on the two-component fits of the R-PDLF data. The measurement of the ²H NMR spectrum alone (from perdeuterated acyl chains) does not provide information concerning the origin of each of the distinct quadrupolar splittings observed. However, R-PDLF enables chemical shift selectivity in the direct dimension. By fitting the indirect dimension time domain of the R-PDLF experiment for the ω, ω-1, and ω-2 carbons with two components (one per acyl chain), one readily enables the assignment of the five lowest quadrupole couplings in the ²H NMR spectrum.

The agreement between the ²H NMR and R-PDLF values highlights the usefulness of modeling R-PDLF data in contrast to the conventional dipolar splitting read-off from the frequency-domain analysis. Such distinct dipolar couplings cannot be observed in the dipolar spectrum due to the limited number of points measured as described previously. In fact, in this case, from the read-off method one obtains for both the ω−1 and ω−2 carbons the lowest order parameters with great precision. The only fitting parameters used to make all the plots shown in Fig. 5 were the dipolar couplings (i.e., order parameters), namely just one fitting parameter in case of a one-component fit and two fitting parameters in case of a two-component fit. In the case of the two-component fits we have assumed that the intensity for each component was equal. This is a reasonable assumption in this case, because the expected differences between the order parameters and the dynamics of the two acyl chains are very small. Therefore, the rINEPT intensity of each component should be very similar according to the theoretical model for the dependence of ¹³C rINEPT intensity on CH order parameter and correlation times described by Topgaard and coworkers (Nowacka et al., 2010, 2013). In general though, two-component fits may require fitting also the relative magnitude of each component for obtaining a good fit. Here we opted for a minimal amount of fitting parameters, i.e., the simplest model enabling to match the experimental data within experimental uncertainty.

For determining the numbers shown in Table 1, we have additionally used a phenomenological single exponential relaxation contribution as described in the Supplement. This did not affect the order parameters of the two-component fits determined for the carbons ω-1 and ω-2. On the other hand, for the ω carbon, without inclusion of the single exponential relaxation, the order parameters calculated were 0.026 and 0.020. The line broadening corresponding to the minimum of χ² for the single-component fits shown in Table 1 and in Fig. 4 using the model that includes RF inhomogeneity was in the range of 4–10 Hz. On the other hand, for the model that does not account for RF inhomogeneity, the resulting line broadening from the minimization of χ² increases significantly in some cases, especially for low frequency dipolar modulations, to as much as 70 Hz. Such an increased relaxation rate is due to the significant deviation of the ideal model from the experimental data as shown in Fig. 5 for the ω carbon. For fits with the ideal model that do not lead to such large relaxation rates, although the minimum of χ² is higher than for the model including RF inhomogeneity (as shown in the Supplement in Figs. S2–S8), the fitting outcome is also good (which is the case for the methylene carbons in Table 1 and for fits of short t₁ intervals as shown in Fig. 4).

Figure 7 Phospholipid headgroup C–H bond order parameter magnitudes determined from a brain lipid extract using time-domain fits accounting for B₁ inhomogeneity. The polarization transfer method used was rINEPT, recycle delay was 5 s, MAS rate was 5 kHz, and a total of 4096 scans were acquired for each point in the indirect dimension which amounts roughly to 2 d of experimental time. (a) Chemical structures and carbon labels of the most abundant phospholipids in the brain lipid extract are phosphatidylethanolamine (PE), phosphatidylcholine (PC), and phosphatidylserine (PS). (b) ¹H-¹³C rINEPT spectrum of the brain lipid extract membranes together with reference spectra for the most abundant lipids. The additional labels c and s were used to identify cholesterol and galactose peaks, respectively. (c) Application of the proposed methodology for fitting the dipolar modulations and determining the corresponding order parameter magnitudes for selected carbons.

3.3 Application to a brain lipid extract

To showcase the applicability of the presented methodology, we investigated the molecular structure of a complex membrane system composed of a brain lipid extract with several lipid types present. The lipid mixture consists of the chloroform : methanol extract of porcine brain tissue as purchased from Avanti Polar Lipids.

Figure 7 shows the ¹³C spectrum of the brain lipid extract at excess hydration together with spectra acquired from single phospholipid model systems. Due to the much higher complexity of the brain lipid extract sample in comparison to model membrane systems (usually composed of one to three distinct components), the number of transients required to achieve a reasonable signal-to-noise ratio is considerably higher. Here, 4096 transients were acquired for each point in the indirect dimension with a total experimental time cost of approximately 48 h (in comparison to the 256 transients used for the simple model systems). The comparison of the ¹³C spectra measured from the brain extract with the spectra from the simple models enables one to identify a number of peaks from different components known to exist in the extract (according to the chemical composition provided by Avanti Polar Lipids for 41.3 wt % of the lipid extract), namely the α, β and γ segments of the headgroup of phosphatidylcholine (PC), the α and β segments of phosphatidylethanolamine (PE), and the α and β segments of phosphatidylserine (PS). From inspection of the ¹³C rINEPT spectrum, and knowing the lipid composition of myelin that makes up most of the dry weight of brain tissue (Norton and Autilio, 1966), it is clear that most of the remaining unknown 58.7 wt % content of the lipid extract is cholesterol plus a considerable fraction of galactose cerebrosides (GalCer) labeled c and s in Fig. 7b, respectively.

(Brown and Seelig, 1978)(Ulrich and Watts, 1994)(Bechinger and Seelig, 1991)(Brown and Seelig, 1978)(Ulrich and Watts, 1994)(Bechinger and Seelig, 1991)(Brown and Seelig, 1978)(Brown and Seelig, 1978)(Brown and Seelig, 1978)(Roux and Bloom, 1990)(Roux and Bloom, 1990)(Roux and Bloom, 1990)

Table 2 Order parameter magnitudes determined from a brain lipid extract using the methodology proposed in this work as shown in Fig. 7. The superscript and subscript for each value denote upper and lower 90 % confidence bounds as shown in the Supplement (Fig. S9). The values are compared to ²H NMR values measured from phospholipid model systems reported previously.

^a $T=\mathrm{40}{}^{\circ }$C, ^b $T=\mathrm{30}{}^{\circ }$C and ^c $T=\mathrm{23}{}^{\circ }$C.

Figure 7c shows the indirect time-domain modulation for a number of identified peaks in the spectra measured by R-PDLF spectroscopy. In this case, due to the large number of transients required to achieve a reasonable signal-to-noise ratio (4096 scans), only very short indirect dimension modulations can be measured. Namely, only eight points in the indirect dimension were acquired for the brain lipid extract sample. Such a reduced total indirect dimension time is not sufficient to obtain dipolar splittings after Fourier transformation. However, application of time-domain fitting including RF inhomogeneity enables the determination of the dipolar couplings and, therefore, the estimation of C–H bond order parameters. The estimated order parameters are presented in Table 2.

The headgroup order parameters determined from the set of different phospholipid types in the brain lipid extract are very much in line with what is observed in corresponding lipid model membranes. The only exception is the PS β order parameter with a rather small order parameter equal to ${\mathrm{0.04}}_{\mathrm{0.031}}^{\mathrm{0.052}}$. This suggests that the brain lipid extract may contain calcium ions that are known to significantly affect the PS headgroup orientation. A thorough analysis of the order parameters in terms of the molecular structure in the complex brain lipid extract system will not be given here. Such analysis should be done in combination with experiments and MD simulations on PC/PE/PS/cholesterol/GalCer model membranes which resemble the lipid extract composition, to investigate potential phase separation behavior, effects of distinct ions present in the system and pH, or ultimately the effect of macromolecules such as basic myelin protein, the most abundant protein in myelin. We believe that such a combined approach will lead to a more accurate description of the lipid/protein interactions in myelin, which may be important in the context of a number of diseases related to demyelination. The methodology proposed here will enable one to perform such investigations with higher accuracy using shorter experiments and without the requirement of isotopic labeling. Moreover, we envision applications to many other complex membranes such as bacterial membrane models, eukaryotic membrane models with cell type or organelle-specific compositions, technological systems such as lipid–nanoparticles, and to address the molecular structure/orientation of drugs, peptides, or other molecules with anisotropic motion in lipid membranes.

4 Conclusions

In summary, we have described a PDLF NMR methodology that consists of performing time-domain fits of the dipolar modulation with numerical simulations that account for the RF spatial inhomogeneity of the probe used. The RF inhomogeneity is experimentally measured and is an a priori constraint. A dipolar modulation of a single component is therefore fitted with only two parameters, the dipolar coupling and a phenomenological relaxation time to account for exponential relaxation. The proposed methodology enables one to determine C–H bond order parameters from simple model systems with a much higher accuracy than previously and enables the investigation of complex lipid membrane systems that were so far inaccessible using the conventional read-off PDLF NMR methodology. We believe that the method presented will be extremely useful in the future concerning molecular structural investigations of complex systems such as multi-component models, lipid extracts, and lipid membranes with drugs, peptides, or other molecules incorporated. Moreover, the higher accuracy and farther reach of the method will also be fundamental to validate molecular dynamics simulations for which structural experimental data were not accessible so far.