Structure of POPC Lipid Bilayers in OPLS3e Force Field

Abstract

It is crucial for molecular dynamics simulations of biomembranes that the force field parameters give a realistic model of the membrane behavior. In this study, we examined the OPLS3e force field for the carbon-hydrogen order parameters S_CH of POPC (1-palmitoyl-2-oleoylphosphatidylcholine) lipid bilayers at varying hydration conditions and ion concentrations. The results show that OPLS3e behaves similarly to the CHARMM36 force field and relatively accurately follows the experimentally measured S_CH for the lipid headgroup, the glycerol backbone, and the acyl tails. Thus, OPLS3e is a good choice for POPC bilayer simulations under many biologically relevant conditions. The exception are systems with an abundancy of ions, as similarly to most other force fields OPLS3e strongly overestimates the membrane-binding of cations, especially Ca²⁺. This leads to undesirable positive charge of the membrane surface and drastically lowers the concentration of Ca²⁺ in the surrounding solvent, which might cause issues in systems sensitive to correct charge distribution profiles across the membrane.

Figure 1

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

Figure 2

Carbon–hydrogen bond order parameters S_CH at full hydration for headgroup, backbone, and acyl chains in simulations and experiments. Experimental values for the POPC ¹H–¹³C NMR at 300 K are from ref (17) and for ²H NMR from ref (96). The ±0.02 error bars of the ¹³C experimental values represent also the range within which most of the published experimental data resides, see discussion in text. For naming of carbon segments, see Figure 1.

Figure 3

Illustration of forking. The θ-angle distribution of the α-carbon C–H bonds for R and S hydrogens (A) in CHARMM36 at full hydration (44 w/l) showing no forking and (B) in CHARMM36 at 5 w/l showing forking, cf. Figure 4. Distributions are calculated over 5 lipids as an example; error bars represent standard error of mean. θ is the angle between a C–H bond and the membrane normal; see Methods for more information. (C) A snapshot with the α-carbon hydrogens marked with arrows.

Figure 4

Response of the headgroup order parameters S_CH^β and S_CH^α to decreasing hydration level. Experimental values for POPC (²H NMR) at 296 K are from ref (99). Notably, small changes in temperature seem not to have a major effect on S_CH, see Figures S10 and S11.

Figure 5

Change of order parameters in the headgroup α (lower panels) and β (upper panels) segments in response to rising concentrations of NaCl (left panels) or CaCl₂ (right panels). Experimental values for DPPC (²H NMR) at 323 and 332 K are from ref (101) and for POPC (²H NMR) at 313 K from ref (102). The out-of-bounds ΔS_CH^β points of OPLS3e in response to CaCl₂ (top right panel) are −0.102 ± 0.0085 and −0.089 ± 0.0090 (500 mM), and −0.13 ± 0.011 and −0.11 ± 0.013 (1000 mM). Corresponding values for ΔS_CH^α (bottom right panel) are −0.10 ± 0.010 and −0.093 ± 0.010 (500 mM), and −0.073 ± 0.016 and −0.097 ± 0.015 (1000 mM). Full figure is shown as the Figure S3. Due to their very slow equilibration (see Figure S4), for the OPLS3e CaCl₂ 200, 500, and 1000 mM concentrations the last 100 ns of the 1 μs simulation was used here. Note that to show possible forking at [salt] = 0, best seen in the bottom left panel for the OPLS3e, the average of the C–H bond order parameters of the R and S hydrogens was used to set the baseline.

Figure 6

Distribution for Na⁺ (solid lines) and Cl^– (dashed) ions along the bilayer normal shown as percentage of salt concentration. Green represents OPLS3e and blue CHARMM36. The graphs were obtained by dividing the number densities with the total salt concentration. Note that both leaflets are plotted (two almost fully overlapping lines) to highlight the symmetry of the ion distributions.

Figure 7

Distribution for Ca²⁺ (solid lines) and Cl^– (dashed lines) ions along the bilayer normal shown as percentage of salt concentration. Green represents OPLS3e and blue CHARMM36. The graphs were obtained by dividing the number densities with the total salt concentration. Because of their very slow equilibration (see Figure S4), for the OPLS3e CaCl₂ 200, 500, and 1000 mM concentrations the last 100 ns of the 1 μs simulation was used here. Note that both leaflets are plotted (two mostly overlapping lines).

Figure 8

(A) Small-angle X-ray scattering (SAXS) form factors at full hydration (44 w/l); experimental data from ref (107). The responses of area per lipid A_L to lowering hydration level (B), added NaCl (C), and added CaCl₂ (D); error bars represent standard deviation. Experimental A_L as a function of lowering hydration for SOPC (NMR and X-ray) are from ref (110) and as a function of [NaCl] and [CaCl₂] for POPC from ref (111).

Figure 9

Responses of lipid lateral diffusion coefficients D_L to lowering hydration level (A), added NaCl (B), and added CaCl₂ (C). The experimental data in panel A are from ref (114), except for the 30 w/l values for 303 and 298 K, which are from ref (119). Note that panels B and C show D_L normalized with respect to the salt-free state, as the relevant experimental data in ref (118) were given as maximal percentage changes (indicated here with the gray rectangles).