Toward Atomistic Resolution Structure of Phosphatidylcholine Headgroup and Glycerol Backbone at Different Ambient Conditions.

Phospholipids are essential building blocks of biological membranes. Despite a vast amount of very accurate experimental data, the atomistic resolution structures sampled by the glycerol backbone and choline headgroup in phoshatidylcholine bilayers are not known. Atomistic resolution molecular dynamics simulations have the potential to resolve the structures, and to give an arrestingly intuitive interpretation of the experimental data, but only if the simulations reproduce the data within experimental accuracy. In the present work, we simulated phosphatidylcholine (PC) lipid bilayers with 13 different atomistic models, and compared simulations with NMR experiments in terms of the highly structurally sensitive C-H bond vector order parameters. Focusing on the glycerol backbone and choline headgroups, we showed that the order parameter comparison can be used to judge the atomistic resolution structural accuracy of the models. Accurate models, in turn, allow molecular dynamics simulations to be used as an interpretation tool that translates these NMR data into a dynamic three-dimensional representation of biomolecules in biologically relevant conditions. In addition to lipid bilayers in fully hydrated conditions, we reviewed previous experimental data for dehydrated bilayers and cholesterol-containing bilayers, and interpreted them with simulations. Although none of the existing models reached experimental accuracy, by critically comparing them we were able to distill relevant chemical information: (1) increase of choline order parameters indicates the P-N vector tilting more parallel to the membrane, and (2) cholesterol induces only minor changes to the PC (glycerol backbone) structure. This work has been done as a fully open collaboration, using nmrlipids.blogspot.fi as a communication platform; all the scientific contributions were made publicly on this blog. During the open research process, the repository holding our simulation trajectories and files ( https://zenodo.org/collection/user-nmrlipids ) has become the most extensive publicly available collection of molecular dynamics simulation trajectories of lipid bilayers.

Introduction

Phospholipids containing various polar headgroups and acyl chains are essential building blocks of biological membranes. Lamellar phospholipid bilayer structures have been widely studied with various experimental and theoretical techniques as a simple model for cellular membranes.[1−8] Phospholipid molecules are composed of hydrophobic acyl chains connected by a glycerol backbone to a hydrophilic headgroup; see Figure for the structure of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC). The behavior of the acyl chains in a lipid bilayer is relatively well understood.[1−5,8,9] The conformations sampled by the glycerol backbone and choline in a fluid bilayer are, however, not fully resolved as even the most accurate scattering and Nuclear Magnetic Resonance (NMR) techniques give only a set of values that the structure has to fulfill, but there is no unique way to derive the actual structure from them.[9−18] Some structural details have been extracted from crystal structures, 1H NMR studies, and Raman spectroscopy,[19−25] but general consensus concerning the structures sampled in the fluid state has not been reached.[9−18,24,25] Importantly, the structural parameters for the glycerol backbone are similar for various biologically relevant lipid species (phosphatidylcholine (PC), phosphatidylethanolamine (PE), and phosphatidylglycerol (PG)) in various environments,[26] and the structural parameters for the choline headgroup are similar in model membranes and real cells (mouse fibroblast L-M cell).[27] Thus, resolving the PC-lipid glycerol and choline structures would be useful for understanding a wide range of different biological membranes.

Figure 1

Chemical structure of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC).

Classical atomistic molecular dynamics simulations have been widely used to study lipid bilayers.[2−7] As these models provide an atomistic resolution description of the whole lipid molecule, they have the potential to solve the glycerol backbone and headgroup structures. The experimental C–H bond order parameters (routinely compared between experiments and simulations for the acyl chains[2−6]) are also known for the glycerol backbone (g1, g2, and g3) and choline (α and β) segments (see Figure for definitions) and are among the main parameters used in attempts to derive lipid structures from experimental data.[10−13,15,16,18] Notably, the structures sampled in a simulation that reproduces these parameters will automatically comprise an interpretation of the experiments. In other words, such simulations can be considered as an accurate atomistic resolution description of the behavior of lipid molecules in a bilayer.

Only a few studies[28−37] have compared the glycerol backbone and choline headgroup order parameters between simulations and experiments. The main reason probably is that the existing experimental data for the glycerol backbone and choline headgroups are scattered over many publications and published in a format that is difficult to understand without some NMR expertise. In addition to the order parameters, dihedral angles for the glycerol backbone and headgroup estimated from experiments,[28,38−42]31P chemical shift anisotropy,[36] and 31P–13C dipolar couplings[43] have been used to assess the quality of a simulation model.

In this work, we first review the most relevant experimental data for the glycerol backbone and choline headgroup order parameters in a phosphatidylcholine lipid bilayer. Then, the available atomistic resolution lipid models are carefully compared to the experimental data. The comparison reveals that the CHARMM36,[31] GAFFlipid,[33] and MacRog[37] models have the most realistic glycerol backbone and choline structures. We also compare the glycerol backbone and choline structures between the most often used Berger-based lipid model[44] and the best performing models, to demonstrate that by using the order parameters we can distinguish the more reasonable structures from the less reasonable ones. However, none of the current models is accurate enough to properly resolve the atomistic resolution structures.

In addition to fully hydrated single component lipid bilayers, the glycerol backbone and choline order parameters have been measured under a large number of changing conditions: hydration level,[45−47] cholesterol content,[35,48] ion concentration,[49−53] temperature,[54] charged lipid content,[52,53] charged surfactant content,[55] drug molecule concentration,[30,56,57] and protein content[58,59] (listing only the publications most relevant for this work and the pioneering studies). Existence of these data allows the comparison of structural responses to varying conditions between simulations and experiments, in other words, validation of the simulation models and interpretation of the original experiments. Here we demonstrate the power of this approach in understanding the behavior of a bilayer as a function of hydration level and cholesterol content. Choline headgroup order parameters as a function of ion concentration, and their relation to the ion binding affinity, are discussed elsewhere.[60]

Methods

Open Collaboration

This work has been done as a fully open collaboration, using the nmrlipids.blogspot.fi blog[61] as a communication platform. Our approach is inspired by the Polymath project;[62] however, there are some essential differences. We started by publishing a manuscript[63] discussing the glycerol backbone and choline structures in a Berger-based model (the most used molecular dynamics simulation model for lipid bilayers). Simultaneously, we presented an open invitation for further contributions and discussion on the blog. All the scientific contributions were made publicly through the blog. Every contributor was offered coauthorship according to the guidelines defined in the beginning of the project;[64] the acceptance of the offer was based on authors’ self-assesment of their scientific contribution. These contributions are summarized in the Supporting Information.

Almost all simulation data, including input files for reproduction and trajectories for further analysis, are collected on our CERN-hosted Zenodo file repository (https://zenodo.org/collection/user-nmrlipids). Thus, in addition to the main topic of this manuscript, we present the most extensive publicly available collection of simulation trajectories for lipid bilayers, opening up numerous possiblities for different analyses with much less effort than previously required. Further information, such as scripts, figures, and manuscript text files, are available through our GitHub repository.[65]

Order Parameters from Experiments

The order parameter of a hydrocarbon C–H vector is defined aswhere the angle brackets denote an ensemble average over the sampled conformations, and θ is the angle between the C–H bond and the membrane normal. The absolute values of order parameters can be measured by detecting quadrupolar splitting with 2H NMR[66] or by detecting dipolar splitting with 1H–13C NMR.[35,67−69] The measurements are based on different physical interactions, and also the connections between order parameter and quadrupolar or dipolar splitting are different. The absolute values of order parameters from the measured quadrupolar splitting ΔνQ (2H NMR) are calculated using the equation , where the value for the static quadrupole splitting constant is estimated from various experiments to be 170 kHz leading to a numerical relation |SCD| = 0.00784 × ΔνQ.[66] The absolute values of order parameters from the effective dipolar coupling dCH (1H–13C NMR) are calculated using the equation , where values between 20.2–22.7 kHz are used for , depending on the original authors.[35,67−69] The effective dipolar coupling dCH is related to the measured dipolar splitting ΔνCH through a scaling factor that depends on the pulse sequence used in the 1H–13C NMR experiment.[35,67−69] It is important to note that the order parameters measured with different techniques based on different physical interactions are in good agreement with each other (see Results and Discussion), indicating very high quantitative accuracy of the measurements. For a more detailed discussion, see ref (70).

The absolute values of order parameters are accessible with both 2H NMR and 1H–13C NMR techniques. However, only 1H–13C NMR techniques also allow the measurement of the sign of the order parameter.[16,67,68] The measured sign is negative for almost all the carbons discussed in this work, except for α which is positive.[16,67,68] For more detailed discussion about the determination of the sign of the order parameters, see ref (71).

For most CH2 segments in a fluid phospholipid bilayer, the order parameters of both hydrogens are equal. However, in some cases (e.g., g1, g3, and C2 carbon in the sn-2 chain) the two order parameters are not equal; this can be observed with both 2H NMR and 1H–13C NMR techniques. In the present work, to avoid confusion with the dipolar and quadrupolar splittings in NMR terminology, we call the phenomenon of unequal order parameters for hydrogens attached to the same carbon forking. Forking has been studied in detail with 2H NMR techniques by deuterating the R or S position in a CH2 segment, and the studies show that it arises from differently sampled orientations of the two C–H bonds, not from two separate populations of lipid conformations.[26,72]

Order Parameters from Simulations

The order parameters from simulations were calculated directly using the definition of eq . For the united atom models, the hydrogen positions were generated post-simulationally from the positions of the heavy atoms and the known hydrocarbon geometries. The statistical error was estimated on the basis of the assumption that different lipids are statistically independent entities, which should be the case in fluid phase: The time average of a given order parameter was first calculated separately for each lipid, then the standard error of the mean over the time averages was taken as the error bar for this order parameter.

It has been pointed out that the sampling of individual dihedral angles might be very slow compared to the typical (100 ns) simulation time scales.[73] After 200 ns, however, even the slowest rotational correlation function of a C–H bond (g1) reaches a plateau (SCH2) in the Berger-POPC-07 model,[74] and notably, the dynamics of this segment have been shown to be significantly slower in simulations than in experiments.[75] In practice, due to averaging over different lipid molecules, less than 200 ns of simulation data should be enough for the order parameter calculation; if the sampling within typical simulation times is not enough for the convergence of the order parameters, then the simulation model in question has unphysically slow dynamics.

Simulated Systems

All simulations are run with a standard setup for a planar lipid bilayer in zero tension and constant temperature with periodic boundary conditions in all directions by using the GROMACS software package[76] (version numbers 4.5.X–4.6.X), LAMMPS,[77] MDynaMix,[78] or NAMD.[79] The number of molecules, simulation temperatures, and the length of simulations of all the simulated systems are listed in Tables , 2 and 3. Full simulation details are given in the Supporting Information (SI) or in the original publications in case the data is used previously therein. All simulation parameters were set as close to the original parametrization works as possible. Additionally, the files related to the simulations and the resulting trajectories are publicly available for almost all systems in the Zenodo collection https://zenodo.org/collection/user-nmrlipids. The references pointing to simulation details and files are also listed in Tables –3.

Table 1

Fully Hydrated Single Component Lipid Bilayer Systems Simulated for Figure : 1,2-Dimyristoyl-sn-glycero-3-phosphocholine (DMPC), Dilauroylphosphatidylcholine (DLPC), Dipalmitoylphosphatidylcholine (DPPC), and 1-Palmitoyl-2-oleoylphosphatidylcholine (POPC)a

force field	lipid	Nlb	Nwc	Td (K)	tsime (ns)	tanalf (ns)	filesg	detailsh
Berger-DMPC-04[80]	DMPC	128	5097	323	130	100	(81)*	(82)
Berger-DPPC-98[83]	DPPC	72	2864	323	60	30	(84)*	SI
Berger-POPC-07(74)	POPC	128	7290	298	270	240	(85)*	(75)
CHARMM36[31]	DPPC	72	2189	323	30	25	(86)*	SI
CHARMM36[31]	DPPC	72	2189	323	130			(31)i
CHARMM36[31]	POPC	72	2242	303	30	20	(87)*	SI
CHARMM36(31)	POPC	128	5120	303	200	100	(88)*	SI
MacRog[89]	POPC	288	12 600	310	100	80	(90)*	SI
MacRog(89)	POPC	128	6400	310	400	200	(91)*	SI
MacRog[89]	POPC	288	14 400	310	90	40	(92)*	SI
GAFFlipid[33]	DPPC	72	2197	323	90	50	(93)*	SI
GAFFlipid[33]	DPPC	72	2167	323	250	250		(33)j
GAFFlipid(33)	POPC	126	3948	303	137	32	(94)*	SI
Lipid14[95]	POPC	72	2234	303	100	50	(96)*	SI
Poger[97]	DPPC	128	5841	323	2 × 100	2 × 50	(98, 99)*	SI
Slipids[100]	DPPC	128	3840	323	150	100	(101)*	SI
Slipids[102]	POPC	128	5120	303	200	150	(103)*	SI
Kukol[104]	POPC	512	20 564	298	50	30	(105)*	SI
Chiu[106]	POPC	128	3552	298	56	50	(107)*	SI
Högberg08[29]	DMPC	98	3840	303	75	50	(108)*	(29)
Högberg08[109]	POPC	128	3840	303	100	80	(110)*	(109)
Ulmschneiders[111]	POPC	128	3328	310	100	50	(112)*	SI
Tjörnhammar14[113]	DPPC	144	7056	323	200	100	(114)*	(113)
Botan-CHARMM36-UA[115]	DLPC	128	3840	323	30	20	(116)	SI
Lee-CHARMM36-UA[117]	DPPC	72	2189	323	70	50	(118)*	SI

The bolded systems were used also for Figure .

Number of lipid molecules.

Number of water molecules.

Temperature.

Total simulation time.

Time used for analysis.

Reference link for the downloadable simulation files; the data sets marked with an asterisk also include a part of the trajectory.

Reference for the full simulation details. The original publication is cited if simulation data from previously published work has been directly used; for other systems the simulation details are given in the Supporting Information.

Magnitudes from Figure S4 of Klauda et al.;[31] signs matched to our simulations.

Magnitudes from Figure 9 of Dickson et al.;[33] signs matched to our simulations.

Table 2

Simulated Single Component Lipid Bilayers with Varying Hydration Levelsa

force field	lipid	nb (w/l)	Nlc	Nwd	Te (K)	tsimf (ns)	tanalg (ns)	filesh	detailsi
Berger-POPC-07[74]	POPC	57	128	7290	298	270	240	(85)*	SI
	POPC	7	128	896	298	60	50	(119)*	SI
Berger-DLPC-13[120]	DLPC	28	72	2016	300	80	60	(121)*	(120)
	DLPC	24	72	1728	300	80	60	(122)*	(120)
	DLPC	20	72	1440	300	80	60	(123)*	(120)
	DLPC	16	72	1152	300	80	60	(124)*	(120)
	DLPC	12	72	864	300	80	60	(125)*	(120)
	DLPC	8	72	576	300	80	60	(126)*	(120)
	DLPC	4	72	288	300	80	60	(127)*	(120)
CHARMM36[31]	POPC	40	128	5120	303	150	100	(88)*	SI
	POPC	31	72	2242	303	30	20	(87)*	SI
	POPC	15	72	1080	303	59	40	(128)*	SI
	POPC	7	72	504	303	60	20	(129)*	SI
MacRog[89]	POPC	50	288	14 400	310	90	40	(92)*	SI
	POPC	44	288	12 600	310	100	80	(90)*	SI
	POPC	25	288	7200	310	100	50	(92)*	SI
	POPC	20	288	5760	310	100	50	(92)*	SI
	POPC	15	288	4320	310	100	50	(92)*	SI
	POPC	10	288	2880	310	100	50	(92)*	SI
	POPC	5	288	1440	310	100	50	(92)*	SI
GAFFlipid[33]	POPC	31	126	3948	303	137	32	(94)*	SI
	POPC	7	126	896	303	130	40	(130)*	SI

The simulation file data sets marked with an asterisk also include part of the trajectory.

Water/lipid molar ratio.

The number of lipid molecules.

The number of water molecules.

Simulation temperature.

The total simulation time.

Time frames used in the analysis.

Reference link for the downloadable simulation files.

Reference for the full simulation details.

Table 3

Simulated Lipid Bilayers Containing Cholesterola

force field	lipid	Nlb	Ncholc	CCHOLd	Nwe	Tf (K)	tsimg (ns)	tanalh (ns)	filesi	detailsj
Berger-POPC-07[74]/Höltje-CHOL-13[35,131]	POPC	128	0	0%	7290	298	270	240	(85)*	(75)
	POPC	120	8	6%	7290	298	100	80	(132)*	(35)
	POPC	110	18	14%	8481	298	100	80	(133)*	(35)
	POPC	84	44	34%	6794	298	100	80	(134)*	(35)
	POPC	64	64	50%	10 314	298	100	80	(135)*	(35)
	POPC	50	78	61%	5782	298	100	80	(136)*	(35)
CHARMM36[31,137]	POPC	128	0	0%	5120	303	150	100	(88)*	SI
	POPC	512	0	0%	23 943	298	170	100	(138)*	SI
	POPC	460	52	10%	23 569	298	170	100	(138)*	SI
	POPC	436	76	15%	23 331	298	170	100	(138)*	SI
	POPC	100	24	19%	4960	303	200	100	(139)*	SI
	POPC	410	102	20%	20 972	298	170	100	(138)*	SI
	POPC	384	128	25%	22 327	298	170	100	(138)*	SI
	POPC	332	180	35%	21 340	298	170	100	(138)*	SI
	POPC	256	256	50%	20 334	298	170	100	(138)*	SI
	POPC	80	80	50%	4496	303	200	100	(140)*	SI
MacRog[89]	POPC	128	0	0%	6400	310	400	200	(91)*	SI
	POPC	114	14	11%	6400	310	400	200	(91)*	SI
	POPC	72	56	44%	6400	310	400	200	(91)*	SI
	POPC	64	64	50%	6400	310	400	200	(91)*	SI
	POPC	56	72	56%	6400	310	400	200	(91)*	SI

The simulation file data sets marked with an asterisk also include part of the trajectory.

The number of lipid molecules.

The number of cholesterol molecules.

Cholesterol concentration (mol %).

The number of water molecules.

Simulation temperature.

The total simulation time.

Time frames used in the analysis.

Reference link for the downloadable simulation files.

Reference for the full simulation details.

Figure 3

Order parameters for POPC glycerol and choline groups from simulations with Berger-POPC-07, MacRog, GAFFlipid, and CHARMM36 force fields (the bolded systems in Table ) together with experimental values. The error bars of simulation data are standard errors of mean (see Methods section for details). The magnitudes for experimental order parameters are taken from Ferreira et al.,[35] the signs are based on the measurements by Hong et al.[16,67] and Gross et al.,[68] and the R/S labeling is based on the measurements by Gally et al.[26]

Results and Discussion

Full Hydration: Experimental Order Parameters for the Glycerol Backbone and Headgroup

The specific deuteration of α, β, and g3 segments of DPPC has been successful, allowing the absolute values of the order parameters for these segments to be measured by 2H NMR.[48−50,54] In addition, the absolute values of order parameters for all glycerol backbone and choline headgroup segments in egg yolk lecithin,[67] DMPC,[16,68,69] DOPC,[141] and POPC[35,141] have been measured with several different implementations of 1H–13C NMR experiments. Furthermore, for some systems the signs of the order parameters have been measured with 1H–13C NMR techniques.[16,67,68] The experimental values of the glycerol backbone and choline order parameters from various publications,[35,50,54,68,69] with the signs measured in refs (16, 67, and 68), are shown in Figure .

Figure 2

Order parameters from simulations listed in Table and experiments for glycerol and choline groups. The experimental values were taken from the following publications: DMPC 303 K from ref (68), DMPC 314 K from ref (69), DPPC 322 K from ref (54), DPPC 323 K from ref (50), POPC 296 K from ref (45), and POPC 300 K from ref (35). The vertical bars shown for some of the computational values are not error bars, but demonstrate that for these systems we had at least two data sets (see Table ); the ends of the bars mark the extreme values from the sets, and the dot marks their measurement-time-weighted average. An interactive version of this figure is available at https://plot.ly/~HubertSantuz/72/lipid-force-fieldcomparison/.

In general there is a good agreement between the order parameters measured with different experimental NMR techniques. Almost all the reported values are within a variation of ±0.02, which is also the error estimate given by Gross et al.[68] for all fully hydrated PC bilayers, regardless of variation in their acyl chain composition and temperature. Exceptions are the somewhat lower order parameters reported from some measurements using 1H–13C NMR.[16,67,141] In these experiments, however, either the reported error bars are relatively large,[16,67] or the spectral resolution is quite low and numerical line shape simulations have not been used in the analysis.[141] As it, therefore, is highly likely that the reported lower order parameters are due to lower experimental accuracy, we exclude them from the present discussion; for more details, see ref (70). Motivated by the high experimental reproducibility, we have highlighted in Figure subjective sweet spots (light blue areas spanning 0.04 units around the average of the extremal experimental values), within which we expect the calculated values of the order parameters of a well-performing force field to fall.

In addition to the numerical values, forking (see Order Parameters from Experiments section) is an important feature of the order parameters. In contrast to the lack of forking in the choline segments α and β, both CH2 segments of the glycerol backbone fork. In the g3 segment, forking is small (≈0.02), and some experiments only report the larger or the average value.[35,50] However, forking is significant for the g1 segment, whose lower order parameter is close to zero, and the larger one has an absolute value of approximately 0.13–0.15. Forking was studied in detail by Gally et al.,[26] who used E. coli to stereospecifically deuterate the different hydrogens attached to the g1 or g3 groups in PE lipids, and measured the order parameters from the lipid extract. This experiment gave the lower order parameter when deuterium was in the S position of g1 or R position for g3. Since the glycerol backbone order parameters are very similar irrespective of the headgroup chemistry (PC, PE, or PG) or lipid environment,[26] it is reasonable to assume that the stereospecifity measured for the PE lipids holds also for the PC lipids.

The most detailed experimentally available order parameter information for the glycerol backbone and choline segments of POPC bilayer is collected by taking the absolute values from ref (35), the signs from refs (16, 67, and 68), and the stereospecific labeling from ref (26), and is shown in Figure .

Full Hydration: Comparison between Simulation Models and Experiments

The order parameters of the glycerol backbone and headgroup calculated from different force fields for various lipids have been previously compared to those from experiments.[28−37] The general conclusion from these studies seems to be that the CHARMM-based,[29,31] GAFFlipid,[33] and MacRog[37] force fields performs better for the glycerol backbone and headgroup structures than the GROMOS-based models.[30,32,34,35] However, none of the studies exploits the full potential of the available experimental data discussed in the previous section, i.e. the quantitative accuracy, known signs, and stereospecific labeling of the experimental order parameters.

To get a general idea of the quality of the glycerol backbone and choline headgroup structures in different models, we calculated the order parameters for these parts from 13 different lipid models (Table ) and plotted the results together with experimental values in Figure . Two criteria were used to judge the quality of the model. (1) There must not be significant forking in the α and β carbons, there must be only moderate forking in the g3 carbons, and there must be significant forking in the g1 carbon. (2) The magnitude should be preferably inside the subjective sweet spots determined from experiments (blue shaded regions in Figure ). The results for each force field with respect to the above criteria are summarized in Figure .

Figure 4

Rough subjective ranking of force fields based on Figure . Here “M” indicates a magnitude problem, “F” a forking problem; letter size increases with problem severity. Color scheme: “within experimental error” (dark green), “almost within experimental error” (light green), “clear deviation from experiments” (light red), and “major deviation from experiments” (dark red). The Σ-column shows the total deviation of the force field, when individual carbons are given weights of 0 (matches experiment), 1, 2, and 4 (major deviation). For full details of the assessment, see Supporting Information.

None of the studied force fields fulfills these criteria completely; however, CHARMM36 is close. This is not surprising since the dihedral potentials in this model are tuned against experiments to better reproduce these order parameters.[31] The next models in the list are CHARMM36-UA[115,117] and Högberg08,[29] which is also not surprising since these models are using the CHARMM bonded potentials for glycerol backbone and choline. The fourth and the fifth models in the list, MacRog[37] and GAFFlipid,[33] have independently determined dihedral potentials. All the models based on GROMOS potentials and Slipids perform less well. In the following sections we subject CHARMM36, MacRog, GAFFlipid, and Berger-POPC-07 to a more careful comparison including the stereospecific labeling (Figure ), atomistic level structure, and responses to dehydration and cholesterol content. These models are selected for more detailed study since they are the best representatives of different dihedral potential parametrization techniques (CHARMM36, MacRog, GAFFlipid), and the Berger-based models are the most used lipid models in the literature.

Full Hydration: Atomistic Resolution Structures in Different Models

The results in the previous section revealed significant differences of the glycerol backbone and choline headgroup order parameters between different molecular dynamics simulation models. However, it is not straightforward to conclude which kind of structural differences (if any) between the models the results indicate, because the mapping from the order parameters to the structure is not unique. In this section we demonstrate that (1) the differences in order parameters indicate significantly different structural sampling, which is strongly correlated with the dihedral angles of the related bonds, and that (2) the comparison between experimental and simulated order parameters can be used to exclude nonrealistic structural sampling in molecular dynamics simulations. The demonstration is done for the dihedral angles defined by the g3–g2–g1–O(sn-1) segments in the glycerol backbone and the N-β-α-O segments in the headgroup. These dihedrals were chosen for demonstration, because significant differences between the models are observed around these segments in Figure . We note that performing a similar comparison through all the dihedrals in all the 13 models would probably give highly useful information on how to improve the accuracy of the models, yet this is beyond the scope of the current report.

The dihedral angle distributions for the g3–g2–g1–O(sn-1) dihedral calculated from different models are shown in Figure . The distribution is qualitatively different for the Berger-POPC-07 model, showing a maximum in the gauche+-conformation (60°) compared to all the other models showing a maximum in the anti conformation (180°). The distributions in all the other models have the same general features, with the main difference being that the fraction of configurations in the gauche–-conformation (−60°) is zero for the MacRog, detectable for the CHARMM36 and equally large to the gauche+ fraction in GAFFlipid. From the results we conclude that most likely the wrongly sampled dihedral angle for the g2–g1 bond explains the significant discrepancy to the experimental order parameters for the g1 segment in the Berger-POPC-07 model (Figure ). In conclusion, models preferring the anti conformation for this dihedral give more realistic order parameters; this is in agreement with previous crystal structure and 1H NMR studies.[19−21,23−25]

Figure 5

Dihedral angle distributions for g3–g2–g1–O(sn-1) dihedral from different models (POPC bilayer in full hydration).

The dihedral angle distribution for the N-β-α-O dihedral calculated from the same four models is shown in Figure . Also for this dihedral there are significant differences in the gauche–anti fractions. The gauche conformations are dominant in CHARMM36, only anti conformations are present in MacRog, while in Berger-POPC-07 and GAFFlipid the gauche and anti conformations have equal probabilities. On the other hand, comparison of α and β order parameters in Figure reveals that for these carbons the CHARMM36 is closest to the experimental results, and it is also the only model that has the correct sign (negative) for the β order parameter. This result is again in agreement with previous crystal structure, 1H NMR, and Raman spectroscopy studies,[19−22] which suggest that this dihedral is in the gauche conformation in the absence of ions.

Figure 6

Dihedral angle distributions for N-β-α-O dihedral from different models (POPC bilayer in full hydration).

These examples show that the glycerol backbone and headgroup order parameters reflect the atomistic resolution structure and that the comparison with experiments allows the assessment of the quality of the suggested structure. We were able to pinpoint specific problems in the structures in different models and suggest potential improvement strategies. If an improved atomistic molecular dynamics simulation model would reproduce the order parameters and other experimental observables (like 31P chemical shift anisotropy[36] and 31P–13C dipolar couplings[43]) within experimental accuracy, it would give an interpretation for the atomistic resolution structure of the glycerol backbone and choline.[10−13,15,16,18] The research along these lines is left, however, for future studies.

Response to Dehydration and Cholesterol Content

In addition to pure phosphatidylcholine bilayers at full hydration, the choline headgroup order parameters have been measured under various different conditions.[30,32,35,45−51,54,55] Also, the order parameters for the glycerol backbone have been measured with 1H–13C NMR in dehydrated conditions,[47] and as a function of anesthetics[30] and glycolipids[32] for DMPC, and as a function of cholesterol concentration for POPC.[35] Due to the high resolution in the NMR (especially 2H NMR) experiments, even very small order parameter changes resulting from the varying conditions can be measured (see ref (70) for more discussion), but as already discussed above, it is not simple to deduce the structural changes from order parameter changes.[15,18] However, comparison of the order parameters between simulations and experiments in different conditions can be used to assess the quality of the force field in different situations, and, if the quality is good, to interpret the structural changes in experiments. Here we exemplify such a comparison for a lipid bilayer under low hydration levels and when varying amounts of cholesterol are included in the bilayer. The interaction between ions and a phosphatidylcholine bilayer will be discussed in a separate study.[60]

Phospholipid Bilayer with Low Hydration Level

Figure shows the published[45−47] experimental order parameters for the glycerol backbone and choline as a function of hydration level. Despite slight differences in temperature and acyl chain composition, the three independently reported data sets for the choline (β and α) segments agree well with each other: Both order parameters increase with decreasing hydration level. The glycerol backbone order parameters (g3, g2, g1), in contrast, have been observed[47] to slightly decrease with dehydration. Note that the original experiments[45−47] measured only absolute values, but here we included the signs measured in separate studies.[16,67,68] Consequently, the negative β order parameter actually increases with dehydration as its absolute value decreases.[45−47]

Figure 7

Effect of dehydration on glycerol and choline order parameters in experiments. The magnitudes of order parameters are measured for DMPC (1H–13C NMR) at 314 K,[47] for POPC (2H NMR) at 296 K,[45] and for DOPC (2H NMR) at 303 K.[46] The signs are based on the measurements by Hong et al.[16,67] and Gross et al.[68] Note that to elucidate the relative change as a function of hydration level, the simulation results were vertically shifted; the shift magnitudes for each of the force fields are listed (SCH + shift) in the y-label.

Lipid bilayer dehydration has been studied also with molecular dynamics simulations,[142−147] typically motivated by the discussion concerning the origin of the “hydration repulsion”.[148−150] Only one[142] of these studies, however, compared their simulation model to the experimental choline and glycerol backbone order parameters. Figure shows these order parameters as a function of hydration level for the CHARMM36, MacRog, and GAFFlipid models (having the most realistic atomistic resolution structures) and a Berger-based model (which is the most used lipid model); note that the simulation results have been vertically shifted to ease the comparison with experimental response to dehydration. Despite some fluctuations, the increase of the choline (β and α) order parameters is seen in all four models. The response of the choline order parameters to dehydration can, therefore, be interpreted to qualitatively agree with experiments. The situation is significantly more complicated for the glycerol backbone: None of the four models produced the experimentally seen trends in all the (g3, g2, g1) segments.

The qualitative agreement of the α and β order parameters with experiments in all four simulation models (Figure ) indicates that, despite the unrealistic structures at full hydration (Figures and 4), the structural response of the choline headgroup to dehydration is somewhat realistic. A likely explanation is that as the interlamellar space shrinks with dehydration, the whole choline group orients more parallel to the membrane. Indeed, upon dehydration the angle between P–N (phoshate phosphorus to choline nitrogen) vector and membrane normal increases for all four models (Figure ). However, the amount of increase depends on the model. Especially, the DLPC simulations with Berger model predict a significantly stronger P–N vector tilt than the other models. The Berger model also has generally larger P–N vector angles, and its choline order parameters are more off from experiments than the other three models (Figure ). Thus, the relatively modest tilting with dehydration predicted by MacRog, CHARMM36, and GAFFlipid is probably more realistic.

Figure 8

Average angle between membrane normal and P–N vector as a function of hydration level calculated from different simulations.

It must be stressed that in the models incapable of reproducing the experimental order parameters the free energy landscape is not correct. Thus, even though the order parameter response to dehydration is qualitatively correct, the energetic response is likely to be incorrect. This may have some influence on dehydration energetic calculations made using the Berger model.[145,147]

The response of the glycerol backbone seems to be more subtle than that of the choline headgroup; none of the four models reproduced the experimental trends upon dehydration with enough accuracy to invite a structural interpretation.

Cholesterol-Containing Phospholipid Bilayer

As cholesterol is abundant in biological membranes and has been suggested to be an important player, for example, in domain formation,[151,152] phospholipid–cholesterol interactions have been extensively studied with theoretical[153−156] and experimental[8,35,48,157] methods. It is widely agreed that cholesterol orders lipid acyl tails and thus decreases the area per molecule (condensing effect), but its influence on the lipid headgroup and glycerol backbone remains debated.[151−153] It has been suggested, for example, that the surrounding phospholipids shield cholesterol from exposure to water by reorienting their headgroups (“umbrella model”)[153] or that cholesterol acts as a spacer between the headgroups to increase their entropy and dynamics (“superlattice model”).[152] Molecular dynamics simulations have supported both the umbrella[156] as well as the superlattice[154] model, in addition to suggesting specific interactions of cholesterol with the glycerol backbone.[155] In these studies[154−156] the responses of the glycerol backbone and choline headgroup to increasing cholesterol content were not, however, compared to experiments.

Figure shows the responses of the choline headgroup (β and α) order parameters of POPC (measured by 1H–13C NMR[35]) and DPPC (2H NMR[48]) to increasing cholesterol content. Again, the two independent data sets agree very well: Only very modest (ΔCH < 0.03) changes occur in the choline order parameters as cholesterol content increases from 0 to 60%. The extreme sensitivity of the high resolution 2H NMR experiments is beautifully demonstrated by the measurable[48] (but barely visible on the scale used in Figure ) cholesterol-induced forking of the α order parameter.

Figure 9

Effect of cholesterol content on the glycerol backbone and choline order parameters in experiments[35,48] and simulations with the Berger-POPC-07/Höltje-CHOL-13, CHARMM36, and MacRog force fields. The signs in the experimental values are based on the measurements by Hong et al.[16,67] and Gross et al.[68] In order to elucidate the relative change as a function of cholesterol content, the simulation results were vertically shifted; the shift magnitudes for each of the force fields are listed (SCH + shift) in the y-label.

We note that the modest (ΔCH < 0.02 for g1; <0.04 for g2, g3; see Figure ) effects of cholesterol on the glycerol backbone order parameters of POPC measured by 1H–13C NMR[35] agree well with the results for phosphatidylethanolamine (PE) measured by 2H NMR.[158] This further supports the ideas that the glycerol backbone structural behavior is independent of the headgroup composition[26] and that the headgroup structure is largely independent of the acyl chain region content unless charges are present.[27]

In addition to the experimental data, Figure shows our results for the CHARMM36 and MacRog force fields and the previously published[35] Berger-POPC-07/Höltje-CHOL-13 results. Note that the simulation data are shifted vertically to ease comparison with experimental responses. As previously pointed out,[35] the Berger-based model seriously exaggerates the effect of cholesterol on the phospholipid glycerol backbone and choline headgroup. In comparison, the choline and glycerol backbone responses of CHARMM36 and MacRog are in better qualitative agreement with experiments. Therefore, to resolve the nature of cholesterol-induced structural changes, we calculated from CHARMM36 the glycerol backbone orientation and dihedral angle distributions at various cholesterol contents (Supporting Information). The only detectable changes are the small decrease of gauche(−) and increase of gauche(+) probability of the g3–g2–g1–O(sn-1) dihedral and slight (less than 5°) change in the glycerol backbone orientation. In conclusion, our results suggest that the significant effects of cholesterol on lipid conformations observed in simulations[154−156] are overestimated by the computational models used; cholesterol only induces very small structural changes in the glycerol backbone.

Finally, it is important to note that the CHARMM36 force field parameters (glycerol backbone dihedral potentials) have been tuned to reproduce the experimental order parameters at full hydration.[31] This approach introduces a risk of overfitting, which would manifest itself as wrong responses to changing conditions. Interestingly, according to our results, tuning did not lead to overfitting problems as far as dehydration or cholesterol content is considered.

Conclusions

The atomistic resolution structures sampled by the glycerol backbone and choline headgroup in phoshatidylcholine bilayers are not known despite vast amounts of accurate experimental data. An atomistic resolution molecular dynamics simulation model that would reproduce the experimental data would automatically resolve the structures, thus giving an unprecedently detailed interpretation of the experimental data. In this work we have collected and reviewed the experimental C–H bond vector order parameters available in the literature. These accurate experimental data are then compared to 13 different atomistic resolution simulation models for a fully hydrated lipid bilayer system, followed by bilayers dehydrated to different extents, and finally bilayers containing various amounts of cholesterol. We are led to the following four main conclusions. (1) The C–H bond order parameters measured with different NMR techniques are consistent. By combining the experimental results from various sources, we concluded that the order parameters for each C–H bond are known with a quantitative accuracy of ±0.02. (2) Comparison of order parameters between experiments and different atomistic resolution models together with structural analysis showed that the order parameters can be used to judge the structural accuracy of a model. Thus, the combination of atomistic resolution molecular dynamics simulations and NMR experiments can be used to resolve the atomistic resolution structures of biomolecules in biologically relevant conditions. This approach can be extended from lipids to, for example, membrane proteins. (3) The review of previous experimental results revealed that when a bilayer is dehydrated the choline order parameters increase. Our simulations suggested that this can be explained by the P–N vector tilting more parallel to the membrane. This strongly supports and complements the idea that charge-induced choline tilting can be measured using order parameter changes.[55,60] (4) Only modest changes of glycerol backbone and choline order parameters are observed experimentally with increasing cholesterol content. When interpreted using the computational lipid model that we found to have the most realistic response to cholesterol, this observation means that cholesterol induces only minor changes in the g3–g2–g1–O(sn-1) dihedral of the glycerol backbone; in other words, there is no major conformational change of the lipid. Besides these four main conclusions, we note that we have created the most extensive publicly available collection of molecular dynamics simulation trajectories of lipid bilayers (https://zenodo.org/collection/user-nmrlipids). The mere existence of this collection opens up numerous possibilities for unforeseen analyses, such as data mining, and rapid testing of ideas with much less computational effort than previously required.

In general, we conclude that, in order to fully utilize the potential of atomistic-resolution classical molecular dynamics simulations in the structural interpretation of high resolution NMR data[159] for lipid bilayers, one must improve the phoshatidylcholine glycerol backbone and choline headgroup parameters of the existing force fields.

This work has been done as a fully open collaboration, using nmrlipids.blogspot.fi as the communication platform. All the scientific contributions have been communicated publicly through this blog.[61]