Using Open Data to Rapidly Benchmark Biomolecular Simulations: Phospholipid Conformational Dynamics

Abstract

Molecular dynamics (MD) simulations are widely used to monitor time-resolved motions of biomacromolecules, although it often remains unknown how closely the conformational dynamics correspond to those occurring in real life. Here, we used a large set of open-access MD trajectories of phosphatidylcholine (PC) lipid bilayers to benchmark the conformational dynamics in several contemporary MD models (force fields) against nuclear magnetic resonance (NMR) data available in the literature: effective correlation times and spin-lattice relaxation rates. We found none of the tested MD models to fully reproduce the conformational dynamics. That said, the dynamics in CHARMM36 and Slipids are more realistic than in the Amber Lipid14, OPLS-based MacRog, and GROMOS-based Berger force fields, whose sampling of the glycerol backbone conformations is too slow. The performance of CHARMM36 persists when cholesterol is added to the bilayer, and when the hydration level is reduced. However, for conformational dynamics of the PC headgroup, both with and without cholesterol, Slipids provides the most realistic description because CHARMM36 overestimates the relative weight of ∼1 ns processes in the headgroup dynamics. We stress that not a single new simulation was run for the present work. This demonstrates the worth of open-access MD trajectory databanks for the indispensable step of any serious MD study: benchmarking the available force fields. We believe this proof of principle will inspire other novel applications of MD trajectory databanks and thus aid in developing biomolecular MD simulations into a true computational microscope-not only for lipid membranes but for all biomacromolecular systems.

Figure 1

C–H bond autocorrelation function g(τ). (A) Idealized illustration of the fast (white background) and the slow (green) mode of the correlation function in solid-state NMR experiments. The fast mode decays to a plateau on which g(τ) = S_CH², while the slow mode gives the final descent to zero. Oscillations at the slow mode region are due to magic-angle spinning. (B) Typical g(τ) obtained from an MD simulation, showing the decay toward S_CH². The gray area under the curve is equal to (1 – S_CH²)τ_e.

Figure 2

Effective correlation times (τ_e, upper panel) and R₁ rates (lower panel) do not markedly depend on the system size. Shown are two CHARMM36 POPC data sets that varied the size while keeping other simulation parameters fixed: ref (44) (blue, system sizes 72 and 648 lipids) and ref (46) (red, system sizes 200, 800, and 1800 lipids). Both data sets are shown normalized against their smallest system. The 15 datapoints shown for each system correspond, from left to right, to the carbon segments γ, β, α, g₃, ..., C17/C15′, C18/C16′, cf. Figure 3.

Figure 3

Effective correlation times (τ_e, top) and R₁ rates (bottom) in experiments (black) and MD simulations (colored) of POPC bilayers in the L_α phase under full hydration. Inset shows the POPC chemical structure and carbon segment labeling. Each plotted value contains contributions from all the hydrogens within its carbon segment; the data for segments 8–11 are only from the sn-2 (oleoyl) chain, whereas the (experimentally non-resolved) contributions of both tails are included for segments 2–3 (2′–3′ in the sn-1 chain) and 16–18 (14′–16′). Simulation results are only shown for the segments for which experimental data were available. For τ_e, a simulation data point indicates the average over C–H bonds; however, if τ_e could not be determined for all bonds, only the error bar (extending from the mean of the lower to the mean of the upper error estimates) is shown. The Berger data for segments γ, C18, and C16′ are left out as the protonation algorithm used to construct the hydrogens post-simulation in united atom models does not preserve the methyl C–H bond dynamics. Table 1 provides further simulation details, while information on the experiments is available at ref (39).

Figure 4

Contributions to the dynamics of the headgroup segments. (A) Zoom on the headgroup τ_e (left panel) and R₁ (right). (B) “Cumulative” R₁(τ) of the γ (top panel), β (middle), and α (bottom) segments. R₁(τ) is obtained, as detailed in Methods, by including in the sum of eq 12 only terms with τ_i < τ. Consequently, at τ → ∞, the R₁(τ) approaches the actual R₁. (C) Prefactor weights α_i from eq 11 of γ (top), β (middle), and α (bottom). Note that panels (B) and (C) show a sliding average over three neighboring data points.

Figure 5

Effect of cholesterol on POPC conformational dynamics. (A) Experimental effective correlation times τ_e (top panels) and R₁ rates (bottom) in 100/0 and 50/50 POPC/cholesterol bilayers at full hydration, see ref (39) for further details. (B) Change in τ_e (Δτ_e, top panels) and R₁ (ΔR₁, bottom), in NMR (black) and MD (color), when the bilayer composition changes from pure POPC to 50% cholesterol. Error estimates for the simulated Δτ_e are the maximal possible based on the errors at 0% and 50% cholesterol; for other data, regular error propagation is used. The Berger Δτ_e is not shown because the available open-access trajectories were too short to determine meaningful error estimates. Table 2 provides further simulation details; for segment labeling, see Figure 3.

Figure 6

Effect of drying on PC headgroup and glycerol backbone conformational dynamics. (A) Experimental effective correlation times τ_e for DMPC at low hydration (from ref (42)) do not significantly differ from the τ_e for POPC at full hydration (from ref (39)). (B) Calculated τ_e for POPC at decreasing hydration in three MD models. Symbols indicate the mean of segment hydrogens if τ_e could be determined for all of them; otherwise, only the error bar (extending from the mean of the lower to the mean of the upper uncertainty estimates) is drawn. The area limited by the error bars shaded for visualization. Note that four Berger data points (24, 16, 12, and 4 w/l) are from DLPC. (C) ¹³C NMR R₁ rates (at ω_C = 125 MHz) of the PC headgroup segments in experiments and simulations: experiments indicate an increasing trend upon dehydration. Experimental POPC (T = 298 K) data at 28 w/l is from ref (39) (solid boxes), POPC (298 K) at 20 and 5 w/l from ref (43) (solid diamonds), and DMPC (303 K) at 13 w/l from ref (42) (open boxes). See Table 3 for simulation details.