Mechanical properties of lipid bilayers from molecular dynamics simulation

Highlights

• Mechanical properties of 12 different bilayers are evaluated from simulation.

• The CHARMM 36 force field yields excellent area and compressibility.

• The relation of surface area and NMR order parameters does not fit simple models.

• Bending constants equal those from flicker experiments, but not others.

• Spontaneous curvatures of leaflets in bilayers are 30% lower than in hexagonal phases.

Abstract

Lipid areas (A_ℓ), bilayer area compressibilities (K_A), bilayer bending constants (K_C), and monolayer spontaneous curvatures (c₀) from simulations using the CHARMM36 force field are reported for 12 representative homogenous lipid bilayers. A_ℓ (or their surrogate, the average deuterium order parameter in the “plateau region” of the chain) agree very well with experiment, as do the K_A. Simulated K_C are in near quantitative agreement with vesicle flicker experiments, but are somewhat larger than K_C from X-ray, pipette aspiration, and neutron spin echo for saturated lipids. Spontaneous curvatures of bilayer leaflets from the simulations are approximately 30% smaller than experimental values of monolayers in the inverse hexagonal phase.

Introduction

The equilibrium surface area per lipid, A_ℓ, bilayer area compressibility, K_A, bilayer bending constant, K_C, and the monolayer spontaneous curvature, c₀, are critical mechanical properties of biological membranes. They determine the thickness of the membrane, its ability to compress, expand, or bend, and its propensity to curve. Perhaps surprisingly, there is considerable uncertainty in the experimental values of these properties even for single component bilayers. This is partly because experiments are often well tuned for certain lipids and not others. For example, X-ray methods (Nagle and Tristram-Nagle, 2000) work well for obtaining surface areas of lipids with phosphatidylcholine (PC) head groups, but not for those with phosphatidylethanolamine (PE) head groups. K_A are typically obtained by pipette aspiration methods (Evans and Rawicz, 1990), but these have been subject to revision both to take into account undulations (Rawicz et al., 2000), and more recently, “fast” stretching (Evans et al., 2013); they are also not available for many lipids. K_C determined by X-ray and pipette aspiration are comparable, but can differ by over a factor of two from K_C of the same lipid determined from vesicle flicker experiments (Marsh, 2006, Nagle et al., 2015); the source of the discrepancy is unclear. Lastly, experimental measurements of c₀ are not even obtained from bilayers. Rather, measurements are carried out on the inverse hexagonal (H_II) phase (Gruner et al., 1986, Marsh, 2006), and the results are frequently extrapolated to the lamellar phase.

The quality of all-atom molecular dynamics (MD) simulations of lipid bilayers has improved dramatically since their advent in the early 1990s (Pastor, 1994). Due to increases in computer power and algorithms, state-of-the-art trajectory lengths have increased from hundreds of picoseconds to hundreds of nanoseconds on standard laboratory clusters, to tens of microseconds on special purpose computers (Shaw et al., 2008). These computer advances have enabled rigorous refinements of force fields (FF). Lastly, there have been important formalistic developments for evaluating pressure profiles in systems with long range electrostatic interactions (Sonne et al., 2005), pressures in droplets and cylinders (Sodt and Pastor, 2012), and bending constants (Watson et al., 2012). These advances enable the determination of bilayer mechanical properties for all-atom MD models, and make it reasonable to propose that MD simulations may help to resolve some of the uncertainties and disagreements associated with experimental measurements.

This paper focuses on A_ℓ, K_A, K_C, and c₀ for the CHARMM (Chemistry at HARvard Macromolecular Mechanics) force field C36 (Klauda et al., 2010) for the lipids listed in Table 1. Table 1 also contains relevant nomenclature and abbreviations. Only fully hydrated single component bilayers are considered given the near absence of data for compressibility and bending moduli for those at low hydration or with more than one lipid type. Though values for some of these quantities have been previously published, simulations for most systems have been rerun or extended for uniformity of analysis. Specifically, relatively small systems (72 or 80 lipids) were all simulated for 420 ns, and larger ones (usually 648 lipids) were simulated for at least 120 ns (see Fig. 1 for representative snapshots of the DPPC bilayers). The set of polyunsaturated lipids, PDPC, PDPE, and SDPE, is entirely new. Additionally, previously published values of c₀ for DOPE and DOPC (Sodt and Pastor, 2013) relied on the experimental monolayer bending constant from measurements on H_II phases; c₀ for PSM (Venable et al., 2014) was estimated using the polymer brush model (Rawicz et al., 2000), a popular model that relates K_C to K_A and bilayer thickness. Here, all values of c₀ are evaluated directly from the simulations alone, without any input of experimental properties or assumption of empirical relations among the various mechanical properties.

By way of outline, the Background and Methods Section reviews the strategy used to parameterize C36 (Section 2.1), presents the critical formulae for the calculation of each of the preceding mechanical properties (Sections 2.2 Calculation of, 2.3 Calculation of K, 2.3.1 Membrane bending energetics, 2.3.2 Calculating K, 2.4 Calculation of), and provides the relevant details on the simulations (Section 2.5). The Results and Discussion presents the calculated mechanical values for the lipids in Table 1 (the first part of each section), and compares them with available experiments (the second part). Section 3.1 concerns A_ℓ and K_A, the relationship of the deuterium order parameter to area (Section 3.1.3), and the effects of time step and smoothing Lennard-Jones interactions (Section 3.1.4). Section 3.2 focuses on K_C, and includes a simulation-based test of the polymer brush model (Section 3.2.3), and further evidence that the pressure tensor-based method for calculating the Gaussian curvature modulus may be flawed (Section 3.2.4). Section 3.3 presents values of c₀, and comments on the notion of lipid shape (Section 3.3.3).