## Tuesday, February 25, 2014

### Accuracy of order parameter measurements

The first comment in this blog after publication of the manuscript was written by Alexander Lyubartsev, who asked a very relevant question: how accurate are the measured (and simulated) order parameters? Troughout the blog and manuscript I have assumed the absolute values of the experimentally measured order parameters to be accurate within $$\pm$$ 0.02. I also stated in my guest post to the MARTINI group blog that the measured order parameters are "reliable and not sparse, and they do not depend on the technique used and the lab in which they were measured." In this post I justify my views. In addition, I will explain why the relative changes with changing conditions can be measured with higher accuracy than the absolute values.

### ²H NMR and ¹³C NMR: two fundamentally different experiments

The order parameter of a C-H bond vector can be assessed with two fundamentally different NMR techniques: measuring the quadrupolar splitting from ²H spectra or measuring the dipolar splitting from ¹³C spectra.

In the ²H experiment the quadrupolar splitting $$\Delta \nu_Q$$ is measured. The C-²H order parameter is
$$S_\mathrm{CD}=\frac{4}{3}\frac{h}{e²qQ} \Delta \nu_Q. (1)$$
Here $$e$$ is the elementary charge, $$Q$$ is the deuteron quadrupole moment and $$h$$ is the Planck's constant. The parameter $$q$$ is related to the largest electric field gradient and in practise its value is not known; therefore the static quadrupolar coupling constant $$\frac{e²qQ}{h}$$ is defined, and its value measured for different compounds in their solid state. In C-D order parameter measurements for lipids, it is typical to use the value  measured for different alkenes, $$\frac{e²qQ}{h}$$=170 kHz. The relation between order parameters and quadrupolar splittings then becomes
$$S_\mathrm{CD}=0.00784 \times \Delta \nu_Q.$$
This relation is useful as many publications report only the quadrupolar splittings. For a review and more accurate description see the work of Seelig.

In the ¹³C NMR experiment the effective dipolar coupling $$d_\mathrm{CH}$$ is measured. It must be stressed that $$d_\mathrm{CH}$$ is a fundamentally different physical observable than the $$\Delta \nu_Q$$ of the ²H experiment. The ¹³C-H bond order parameter is
$$S_\mathrm{CH}=\frac{4\pi\langle r_\mathrm{CH}^3 \rangle}{\hbar \mu_0 \gamma_h \gamma_c} d_\mathrm{CH}. (2)$$
Here $$r_\mathrm{CH}$$ is the C-H distance, $$\mu_0$$ is the vacuum permittivity, and $$\gamma_h$$ and $$\gamma_c$$ are the gyromagnetic constants for ¹H and ¹³C nuclei. In contrast to Eq. (1), all the parameters in Eq. (2) are in principle known. However, for the internuclear distance only the average $$\langle r_\mathrm{CH} \rangle$$ is known, not the third moment $$\langle r^3_\mathrm{CH} \rangle$$. For this reason values between 20.2-22.7 kHz are used for $$\frac{\hbar \mu_0 \gamma_h \gamma_c}{4\pi\langle r_\mathrm{CH}^3 \rangle}$$ [Hong et al.,Gross et al.,Dvinskikh et al.]

²H NMR and ¹³C NMR are fully independent experiments since the deuterium quadrupolar splitting $$\Delta \nu_Q$$  and the effective dipolar coupling $$d_{CH}$$ are different physical observables. In addition, the prefactors connecting the observables to the order parameter (Eqs. (1) and (2)) are independently measured. In ¹³C NMR further independent experiments can be performed by measuring the ¹³C dipolar couplings using different pulse sequences [Hong et al., Gross et al., Dvinskikh et al, Ferreira et al.].

The order parameter measurements are very reproducible. The results obtained for PC bilayers with different ¹³C pulse sequences and ²H measurements have been compared in Table 1 in Gross et al., Table 1 in Dvinskikh et al., and Table 1 and Fig. 3 in Ferreira et al.  The order parameters measured with different techniques by independent reasearch groups are in good agreement. The results for the headgroup and glycerol region measured from samples with purified PC lipid molecules are within a maximum difference of 0.02. The larger differences (maximum 0.04) in the tables arise from the samples with E. Coli or eggPC extracts. In the acyl chains there are some carbons for which the difference is sligtly larger even for purified lipid samples, but as we are focusing on the headgroup and glycerol order parameters of purified PC lipids, we can consider these to be known with absolute accuracy of $$\pm$$ 0.02.

Note. In the work by Aussenac et al. (discussed in the context of the blog here, and compared to the simulations by Högberg et al. as pointed out in here) the order parameters are measured for bicelles, not vesicles or bilayer stacks. The order parameters for bicelles are known to be smaller than the ones for vesicles by roughly a constant factor, most likely due to the bilayer orientations with respect to the magnetic field [Raffard et al., Sanders et al.]. When corrected with this constant factor, the order parameters by Aussenac et al. are also close to the ones measured for vesicles.

### Absolute vs. Relative accuracy

The absolute accuracy of measured order parameters is largely affected by the accuracy of the prefactor connecting the measured splittings to the order parameters (Eqs. (1) and (2)). When considering the change of an order parameter with a changing external condition (hydration, ion concentration, etc.), the prefactor can be considered to be unchanging. Therefore, the change in a given order parameter, when measured by the same people with the same technique and equipment, can be determined in a much higher accuracy than the absolute value of the order parameter. We refer to this as the relative accuracy. It is determined by the accuracy of the splitting measurement. Especially the very high spectral resolution of the ²H-measurements allows the measurement of very small order parameter changes. Let us exemplify this with a classical experiment by Akutsu et al., where the effect of different ions on the headgroup $$\alpha$$- and $$\beta$$-quadrupolar splittings was measured, see Fig. 1. Figure 1. The measured quadrupolar splittings with the explanation in the original caption. Reprinted with permission from Akutsu and Seelig, Biochemistry, 20, 7366 (1981). Copyright (1981) American Chemical Society.

Clearly, the effects of different ions on the quadrupole splittings are differentiable within experimental accuracy in Fig. 1. (These changes were later shown to be consistent with the addition of different charges into the bilayer, and the electrometer concept was introduced to measure the amount of charge incorporated in the bilayer interface. [Scherer et al.]) Importantly, when the results of Fig. 1 are transformed to order parameters in Fig. 2, one can see that the numerical changes in $$S_\mathrm{CD}$$ are rather small: the y-axis only spans < 0.03 units for $$\beta$$ and 0.06 for $$\alpha$$. Figure 2. Same data as Fig. 1. but given as order parameters instead of splittings.

Indeed, the measured (upon addition of 1M CaCl$$_2$$ concentration, the changes due NaCl are negligible) changes in the $$\beta$$-and $$\alpha$$-carbons are only $$\sim$$ 0.02 and $$\sim$$ 0.05, respectively.

The resolution of ¹³C experiment depends somewhat on the pulse sequence used. For example Dvinskikh et al. measured the effect of dehydration in the lipid order parameters with APM-CP pulse sequence, and changes of 0.01 were measurable. Fig. 3 shows the order parameters for $$\alpha$$-and $$\beta$$-carbons as a function of dehydration independently measured with ¹³C and ²H NMR for different PC lipids. All the results are qualitatively similar, that is, they give the same relative change with dehydration. Figure 3. Dehydration changes on $$\alpha$$ and $$\beta$$ order parameters measured with different methods. The data taken from Dvinskikh et al., Ulrich et al. and Bechinger et al.

### Why Simulate?

From the experiments it is not clear how large conformational and energetical changes these small changes in the order parameters indicate. This could, however, be revealed with a molecular model that reproduces the measured order parameters. We have already succesfully applied this approach to explain the measured low order parameters (and their temperature dependence) in the lamellar phase of the nonionic surfactant C$$_{12}$$E$$_5$$ [Ferreira et al.]. For the lipid headgroup and glycerol region the problem is, however, that with the current force fields the order parameters from simulation models are both absolutely and relatively too far from the experimental values to be useful for interpretating the experiments.

1. Hey,

my name is Alexander Vogel and I'm working a lot on membranes, order parameters, solid-state NMR and MD simulations (actually comparing them a lot). I didn't read everything yet but here are a few quick thoughts:

1. You might want to check out one of my papers ;-) It is about the comparison between MD simulation (CHARMM C32 forcefield) and experiment of DPPC headgroup conformations...might be interesting for you:

Headgroup Conformations of Phospholipids from Molecular Dynamics Simulation: Sampling Challenges and Comparison to Experiment
Alexander Vogel, Scott E. Feller
J Membrane Biol (2012) 245:23–28
DOI 10.1007/s00232-011-9411-5

I can send you the pdf if you don't have access.

2. I have a recent 200ns POPC simulation lying around (relatively small system...just 50 POPC, 301K). This time with the CHARMM C36 forcefield...which might be the best lipid forcefield at the moment. It would be interesting how that compares to the experimental data. If the results are interesting and you are interested I might even run a simulation more similar to yours (128 lipids, didn't find the temperature you were using) just too see if there are differences between the forcefields and if CHARMM C36 does a better job.

If you are interested send me an email: alexander.vogel {at} medizin.uni-leipzig.de

Best regards,

Alexander

1. Hi Alexander, thank you for your interest!

In reply to your second point. Currently we have fully hydrated CHARMM36 simulation data for DPPC at 323K (one short 25 ns data set) and for POPC at 303K (one short 20 ns, and one 100 ns). These actually seem to match the experimental data quite well, as discussed in detail in this post that was published today:

http://nmrlipids.blogspot.fi/2014/03/the-lipid-forcefield-comparison.html

The plot (Fig. 1) in that posting sums up all the full hydration data we have gathered so far.

Considering your data set. It would definitely be interesting, as it is an independent data set, and as CHARMM36 appears to be among the best performing force fields we have tested. If you could run the order parameter analysis on your trajectory and post the results here, it would be very nice. The script Samuli used for CHARMM36 is available here (it expects as input the trajectory as a Gromacs gro-file called runPROT.gro):

https://www.dropbox.com/s/ub2gitr20568qed/calcORDPcharmmPOPC.sh

If there are any problems, please do not hesitate to ask for clarification!

2. Thanks for pointing out the publication! I was not aware of it, and found it a useful read. The possible sampling issue is very relevant.

We have, actually, recently compared the rotational dynamics of C-H vectors in MD to NMR experiments. This work was partially published in the thesis of Tiago Ferreira:

In this work we suggested that the effective correlation time (integrated area below the rotational correlation function) can be estimated by combining different NMR observables (see e.g. thesis sections 5.5 and 6.3.2). The comparison between the effective correlation times from experiments and Berger MD simulations gave the following result:

https://www.dropbox.com/s/u8cxov46non1pk9/effCTresult.pdf
(we are still optimizing the experimental methodology, but the results should be within +-30% error)

This indicates that dynamics of the glycerol region are significantly too slow in the Berger model. I am not sure how it is in CHARMM, but based on the dispersion results reported by Klauda et al. (http://dx.doi.org/10.1021/jp101759q, Fig. 8) it seems possible that dynamics might not be so well reproduced closer to the headgroup. It is not clear if the dynamics are too fast or too slow, though.

In conclusion, I think that the problem in sampling might arise from too slow dynamics in the model compared to experiments.

To this end, I think it would be highly interesting to calculate the effective correlation times also for the CHARMM model.

3. The above mentioned manuscript, where the experimental effective rotational correlation times are compared to the Berger simulation, is now published:

Model-free estimation of the effective correlation time for C–H bond reorientation in amphiphilic bilayers: 1H–13C solid-state NMR and MD simulations
T. M. Ferreira, O. H. Samuli Ollila, Roberta Pigliapochi, Aleksandra P. Dabkowska and Daniel Topgaard
J. Chem. Phys. 142, 044905 (2015)

http://dx.doi.org/10.1063/1.4906274

2. About the accuracy in simulations: We now calculate the error bar from the simulations by calculating the error of the mean for the average over individual lipids in the system. This is done since all the lipids should sample the same conformations with equal weights in a long simulation enough, i.e., in a long enough simulation the order parameters from each lipid should be the same.

There are still some theoretical arguments agains this choice, e.g. the lipid structures may be still affected by the initial structure, or there may be too long correlations between neighboring lipids compared to simulation time. These arguments are presented from time to time, however, it is very unclear to me how well people really understand the strengths of the correlations. Anyway, I think that current choice is least affected by the above arguments.

To compare the results and error bars in two indepent simulations, I took the CHARMM36 POPC simulations ran by myself and Hubert Santuz. These are indepently done and with different amount of data and lipids. The results with error bars are shown here: https://www.dropbox.com/s/tq3nv68c1tmw04l/ERRORbarCOMP.jpg?dl=0

The error bars from different results are overlapping (with very small overlap in other beta hydrogen, though), so I think that the choice of error bar calculation method supports the current purpose of this work.