A Critical Comparison of Biomembrane Force Fields: Structure and

Subscriber access provided by SUNY DOWNSTATE

Article

A Critical Comparison of Biomembrane Force Fields: Structure and Dynamics of Model DMPC, POPC, and POPE Bilayers Kristyna Pluhackova, Sonja Kirsch, Jing Han, Liping Sun, Zhenyan Jiang, Tobias Unruh, and Rainer A Böckmann J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b01870 • Publication Date (Web): 01 Apr 2016 Downloaded from http://pubs.acs.org on April 2, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A Critical Comparison of Biomembrane Force Fields: Structure and Dynamics of Model DMPC, POPC, and POPE Bilayers Kristyna Pluhackova,∗,† Sonja A. Kirsch,† Jing Han,† Liping Sun,† Zhenyan Jiang,† Tobias Unruh,‡ and Rainer A. Böckmann∗,† †Computational Biology, Department of Biology, Friedrich-Alexander-University of Erlangen-N¨ urnberg, Staudtstr. 5, 91058 Erlangen ‡Lehrstuhl f¨ ur Kristallografie und Strukturphysik, Department Physik, Friedrich-Alexander-University of Erlangen-N¨ urnberg, Staudtstr. 3, 91058 Erlangen E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Atomistic molecular dynamics simulations have become an important source of information for the structure and dynamics of biomembranes at molecular detail difficult to access in experiments. A number of force fields for lipid membrane simulations have been derived in the past, the choice of the most suitable force field is, however, frequently hampered by the availability of parameters for specific lipids. Additionally, the comparison of different quantities among force fields is often aggravated by varying simulation parameters. Here, we compare four atomistic lipid force fields, namely the united-atom GROMOS54a7 and the all-atom force fields CHARMM36, Slipids, and Lipid14, for a broad range of structural and dynamical properties of saturated and monounsaturated phosphatidylcholine bilayers (DMPC and POPC) as well as for monounsaturated phosphatidylethanolamine bilayers (POPE). Additionally, the ability of the different force fields to describe the gel-liquid crystalline phase transition is compared and their computational efficiency estimated. Moreover, membrane properties like the water flux across the lipid bilayer and lipid acyl chain protrusion probabilities are compared.

Introduction Biomembranes, thin layers composed of phospho- and glycolipids, sterols and membrane proteins, form boundaries between cells and allow for their compartmentalization. They are also the stage of action for many processes like signalling, protein dimerization, and enzymatic reactions. Additionally, selective both active and passive transport of metabolites across biomembranes is enabled by integral membrane proteins. Lipid bilayers have been studied in experiments for already more than 90 years 1 and a broad range of membrane properties have been determined. Over the years, the information obtained by standard experimental techniques like calorimetry (applied to study membrane phase behavior) has been extended by sophisticated methods like FRET, 2 AFM, 3 2 H-NMR, 4–6 2

13


C-NMR, 7–9 SANS, 10 SAXS, 11

Page 2 of 44

Page 3 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


QENS, 12 and mass spectrometry (lipidomics). 13 Recently, the development of the orientedsample solid-state

31

P NMR technique 14 enabled an extended view on lipid orientation and

altered lipid dynamics under the influence of charged components 15 as well as in the vicinity of membrane-bound peptides. 16–19 Subsequently, molecular dynamics (MD) simulations have been shown to be able to reproduce the spectral features of those experiments 20 and to identify the molecular interactions underlying the experimental observations. Complementary with the knowledge of the complex biomembrane composition 21 force fields suitable for MD simulations of biomembranes have experienced a boost both in terms of the number of different force fields and of the diverse components of biomembranes (for a recent review on atomistic and coarse-grained simulations of biomembranes see Pluhackova and Böckmann 22 ). Few studies focus on a comparison of lipid force fields for different observables. A recent review 23 summarizes the latest progress, challenges and perspectives of such a comparison. In 2008, lipid structural properties and lipid self-diffusion of a fully hydrated DOPC lipid bilayer were compared for the CHARMM27 force field, 24 the Berger force field, 25 and a therein parameterized generalized Amber force field (GAFF). 26 Later on, Poger and Mark compared the structure and dynamics of lipids using three closely related united-atom lipid force fields (GROMOS54a7, 27 GROMOS53a6-Kukol, 28 and the widely applied variant of the Berger force field 29 ). The first comparison between united and all-atom force fields for both saturated and unsaturated phosphocholine lipids was presented by Piggot et al. in the same year: 30 The authors compared CHARMM36, 31 Berger, 25 GROMOS43a1-S3, 32 GROMOS53a6, 33 and GROMOS53a6-Kukol 28 and the influence of simulation parameters on membrane properties. This seminal work frequently serves as a guide in the choice of proper simulation parameters like the simulation time step, Coulomb and Lennard-Jones cut-offs, water model, etc. for a given force field. A recent broad initiative to compare membrane properties from different force fields (13 in total) is the NMR lipid blog platform (nmrlipids.blogspot.com) launched by S. Ollila. 34,35 The current work focused mainly on the choline and glycerol order parameters of phosphocholine lipids as well as the effect of

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

hydration and cholesterol concentration on those order parameters. Here, we compared four recent atomistic force fields suitable for simulations of biomembranes (i.e. containing parameters for phospholipids and sterols and compatible with parameters for proteins and small molecules), namely the CHARMM36, 31 GROMOS54a7, 27 Lipid14, 36 and Slipids 37 force fields. The latter two can be combined with proteins and small molecules described by force fields from the Amber family, e.g. the most recent Amber 14. 38 All simulations were performed consistently in GROMACS version 5.0.4 39 using force field specific parameter files including common integrators and other parameters assuring the comparability of the observed data and avoiding program version or integrator based biases. This comparative study focuses on both structural and dynamic properties of pure phosphocholine (namely DMPC and POPC) and phosphoethanolamine (POPE) bilayers at different temperatures. Structural properties as well as lipid diffusion, the water flux across the bilayer, lipid hydration, and lipid protrusions were assessed. The latter properties play an important role in membrane fusion 40–42 and may be influenced e.g. by the addition of antimicrobial or fusion peptides. 43,44 Moreover, the phase transition temperature of the gelto-liquid-crystalline phase transition, as well as the ability of a given force field to form a gel phase are compared for a DPPC bilayer using different force fields. Finally, the computational efficiency of the individual force fields is analyzed.

Methods Simulated systems Pure phospholipid bilayers consisting of either fully saturated 1,2-dimyristoyl-sn-glycero-3phosphocholine (DMPC), monounsaturated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) or 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) lipids were solvated by 50 w.w.% water each. The number of lipids was chosen such that the size of 4


Page 4 of 44

Page 5 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a rectangular membrane patch is at least 10 nm. Namely, in each system 384 lipids were present, i.e. 192 lipids per layer. Detailed information about the compositions of the individual systems is listed in 1. Table 1: Summary of simulation systems: number of water molecules, lipid:water ratios and final sizes of the simulation box. Each system contained 384 lipids. a The size of the box in x and y directions (corresponding to a membrane plane) equals. Property lipid:water [w.w %] lipid:water [mol %] number of waters final box-X size [nm]a final box-Z size [nm]

DMPC 300 K 50:50 1:34 13,000 10.07 7.73

DMPC 320 K 50:50 1:34 13,000 10.99 6.72

POPC 300 K 50:50 1:38 14,572 11.17 7.23

POPC 320 K 50:50 1:38 14,572 11.35 7.12

POPE 308 K 50:50 1:36 13,768 10.17 8.26

System preparation The simulation systems (with exception of the systems for the estimation of melting temperature) were prepared following our previously established procedure. 45 In short, the membranes were generated at coarse-grained (CG) resolution by the tool insane 46 and energy minimized by steepest-descent algorithm for 500 steps. Following a short equilibration, 100 ns CG MD simulation for each lipid type and temperature was performed to relax the membranes. The resulting CG membranes were converted to different atomistic force field descriptions by the tool backward 47 and were simulated atomistically for 200 ns at different temperatures (see 1). The phase behaviour of bilayer patches was analyzed by heating up a lipid bilayer in gel phase. A system composed of 288 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) lipids and surrounded by explicit water (1:40 lipid:water molar ratio) was prepared within the CHARMM36 force field and subsequently converted to the other three studied force fields. For every force field, the conversion was followed by a 100 ns equilibration simulation at 300 K in order to obtain lipid bilayer structures in gel phase. For the GROMOS54a7 force field (and to a lower extent also for the Slipids force field) some disruptions of the 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 6 of 44

Page 7 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Summary of simulation conditions for individual force fields. Following abbreviations were used: GRF, generalized reaction field; PME, particle-mesh Ewald; PR, Parrinello-Rahman thermostat; vdW, van der Waals. a The Verlet cutoff scheme is equal to the old switch scheme for van der Waals interactions. b The COM (center of mass) motion was removed linearly for the whole system. c LINCS 59 were applied. d CHARMM TIP3P water model with Lennard-Jones interaction sites on hydrogens. Parameter time step [fs] RCoulomb [nm] Coulomb method RvdW [nm] vdW method dispersion correction COM removalb [steps] neighborlist search [steps] barostat τp [ps] thermostat τT [ps] constraintsc water model

GROMOS54a7 27 2 1.4 GRF 49 (ǫRF = 61) 1.4 Verlet cut-offa 48 no 100 10 Berendsen 52 0.5 Berendsen 52 0.1 h-bonds SPC 56

CHARMM36 31 2 1.2 PME 50 0.8-1.2 Verlet cut-offa 48 no 100 10 PR 53 5 Nosé-Hoover 54,55 0.5 h-bonds TIPS3Pd 57

Lipid14 36 2 1.0 PME 50 1.0 Verlet cut-offa 48 EnerPres 51 100 10 Berendsen 52 1 Nosé-Hoover 54,55 0.5 h-bonds TIP3P 58

Slipids 37 2 1.5 PME 50 1.4-1.5 Verlet cut-offa 48 EnerPres 51 100 10 PR 53 10 Nosé-Hoover 54,55 0.5 all-bonds TIP3P 58

Analysis Standard GROMACS analysis tools and in house written scripts were applied. The initial 60 ns of the production simulations were excluded from the analysis for equilibration purposes. The remaining 140 ns (60 – 200 ns) were used for analysis either as whole or split into 20 ns intervals that were subsequently used for block averaging. Area per lipid For bilayers of 10 nm size, curving effects may significantly affect the estimation of the area per lipid. In order to evaluate this effect, a grid was layed on the phosphorus atoms separately for the upper and lower layer. 60 The area of this grid for each layer divided by the number of lipids in each layer was compared to the area per lipid as obtained from the lateral box area divided by the number of lipids in each layer.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Volume per lipid The volume per lipid (VL ) was determined as

VL =

Vbox − nW × VW . nL

Vbox is the volume of the simulation box, nW (nL ) is the number of water molecules (lipids) and VW is the volume of each water molecule obtained from additional lipid-free MD simulations for each water model, temperature, and the specific force field parameters. The obtained water volumes are listed in the Supplementary Table S1. Isothermal area compressibility modulus The isothermal area compressibility modulus (KA ) was determined 61 from

KA =

2kB T hAL i . nL (sAL )2

kB is the Boltzmann constant, T the simulation temperature, hAL i the average area per lipid, and sAL the standard error of AL . Bilayer thickness The thickness of the bilayer was determined on a grid using the tool g lomepro 62 as the phosphorus-phosphorus distance between the layers. Electron density profiles and scattering form factors The electron density profiles (EDP) along the membrane normal were calculated using a 0.2 ˚ A spacing and compared to the scattering density profile (SDP) model fit of the experimental data. 63 The EDP curves were obtained using

EDP (z) =

Z

X 0

Z

Y

ρ(x, y, z) dxdy . 0

8


Page 8 of 44

Page 9 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ρ denotes the 3D electron density. The experimental data were obtained from Norbert Kuˇcerka (personal communication). The X-ray and neutron scattering form factors, which offer direct (fit-free) comparison between experiment and simulations, 64 were calculated using the SIMtoEXP 65,66 program and compared to experimental data therein. Order parameters For the calculation of the order parameters in all-atom systems an in-house altered version of the GROMACS tool g order has been utilized. Here, the order parameter calculation is based on actual angles of the C-H bonds to the membrane normal. For GROMOS (united-atom force field) the standard calculation based on C-C bonds was applied. Order parameters are defined as 1 SCD = h3 cos2 θ − 1i . 2 θ is the angle between the membrane normal and the vector connecting C to its deuterium atom(s), and hi represents the ensemble average. It has been noted before 30 that the standard version of g order can’t be used to evaluate the order parameters in the vicinity of double bonds. A small How to guide describing the way how to deal with g order and double bonds is included in the Supplementary Information. Lipid headgroup orientation First, the average P-N vector orientation relative to the bilayer normal was calculated. For comparison to experiment, the strategy of Prakash and Sankararamakrishnan 67 was followed and the

31

P-13 C dipolar coupling values for headgroup carbons calculated as

∆ν = 12236.5h(3 cos2 θ − 1)r−3 /2i .

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

r is the instantaneous distance between the carbon and phosphorus atoms in ˚ A, θ is the instantaneous angle between the vector of these two atoms relative to the bilayer normal, and the hi brackets denote averaging over time and over lipids. These values can be directly compared to the experimental observables of Sanders. 68 However, data are available for DMPC at 313 K only. Lipid phosphorus hydration Lipid hydration was evaluated as the average number of water oxygens coordinated in a first solvation shell around the phoshorus atom of a given lipid. The radius of the first solvation shell was determined for each simulation separately as the first minimum of the radial distribution function of water oxygen atoms around lipid phosphorus atoms. Lipid tail protrusion Lipid protrusion events were counted if any acyl chain carbon atom was found more outside the bilayer than the phosphorus atom of the given lipid for at least two following frames (corresponding to 20 ps). The protrusion probability is given as the number of protrusion events per lipid in 1 µs. The lipid protrusion efficiency is obtained by multiplying the lipid protrusion probability and the median of the lipid protrusion time. Diffusion coefficients The lateral diffusion coefficients of the different lipid types were computed from the Einstein equation, using the slope of the lipid mean square displacement (MSD) versus time, 1 d h[r (t + t0 ) − r(t0 )]2 it0 . t→∞ 2n dt

Dl = lim

n is the number of translational degrees of freedom (i.e. for lateral diffusion n = 2) and r is the position of the lipid’s center of mass (COM). The MSD of the lipids’ centers of mass was

10


Page 10 of 44

Page 11 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculated separately for each monolayer and averaged over all molecules and over time. In order to obtain better statistics, the calculation was restarted via the initial time t0 every 500 ps. A linear fit between 10 – 50 ns of the MSD was used for the estimation of the diffusion coefficient. Several simulation studies 29,60,69,70 pointed out that a COM motion of monolayers may result in an overestimation of the diffusion coefficient; therefore this contribution was removed during the MSD calculation. Water flux across the membrane Based on the number of water molecules that passed the bilayer the water flux F was estimated as F =

N VW . hAM it

N is the number of water molecules that passed the bilayer, VW the volume of one water molecule, hAM i the average membrane area, and t the analyzed simulation time (140 ns).

Results In this comparative study, DMPC, POPC, and POPE bilayers were simulated using the CHARMM36, 31 GROMOS54a7, 27 Lipid14, 36 and Slipids 37 force fields.

Membrane structural properties The membrane properties area per lipid, the lipid volume, the isothermal area compressibility modulus, and the bilayer thickness analyzed from the bilayer simulations within the studied force fields are summarized in 3 and 4 and compared to literature experimental values. Experimental values like the area per lipid (APL) and bilayer thickness depend very much on the type of experiment and on the model used to extract these properties from the raw data. Accordingly, literature values show a large spread. Within this span, all force 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 44

Table 3: Basic lipid bilayer structural properties: area per lipid (APL) and volume per lipid (VPL). For cases when the experimental properties were determined at different temperatures, but not exactly at the temperature studied here, a linear intrapolation of the data was performed. The corresponding Rsquared values of the linear fit are a−e 0.99. Results for bilayers in a wrong phase are given in brackets (see text). System

Force Field

POPC 300K

CHARMM36 GROMOS54a7 Slipids Lipid14 CHARMM36 GROMOS54a7 Slipids Lipid14 (CHARMM36) GROMOS54a7 Slipids (Lipid14) CHARMM36 GROMOS54a7 Slipids Lipid14 CHARMM36 (GROMOS54a7) Slipids Lipid14

POPC 320K

DMPC 300K

DMPC 320K

POPE 308K

Curved APL (ACL ,˚ A2 ) 66.6 ± 0.4 63.2 ± 0.2 66.7 ± 0.4 67.9 ± 0.4 68.8 ± 0.2 64.9 ± 0.3 69.6 ± 0.4 70.9 ± 0.2 (57.7 ± 0.6) 60.8 ± 0.4 62.8 ± 0.4 (57.9 ± 0.1) 64.5 ± 0.3 62.9 ± 0.3 65.7 ± 0.4 66.8 ± 0.3 57.7 ± 0.5 (51.6 ± 0.4) 58.5 ± 0.3 56.2 ± 0.3

APL (AL ,˚ A2 ) 63.4 ± 0.4 60.5 ± 0.1 63.7 ± 0.3 64.6 ± 0.1 65.4 ± 0.2 62.2 ± 0.3 66.2 ± 0.4 67.3 ± 0.2 (52.5 ± 0.5) 57.9 ± 0.4 59.8 ± 0.4 (53.7 ± 1.0) 61.2 ± 0.2 60.2 ± 0.2 62.6 ± 0.3 63.5 ± 0.2 55.5 ± 0.4 (50.2 ± 0.3) 56.4 ± 0.3 54.4 ± 0.3

Exp APL (AL ,˚ A2 ) 63.8(300K) 64,a 63.0(297K) 71 68.3(303K) 72 66.0(310K) 73 66.5(320K) 64,a 66.0(310K) 73 59.2(300K) 64,b 59.4(300K) 74,c 59.6(303K) 75 60.6(303K) 72 63(320K) 64,b 64.3(320K) 74,c

58.0(308K) 76 59.7(308K) 77

VPL (VL ,˚ A3 ) 1225.1 ± 0.3 1224.6 ± 0.2 1208.0 ± 0.3 1194.3 ± 0.4 1245.4 ± 0.1 1238.4 ± 0.3 1226.7 ± 0.2 1214.2 ± 0.1 (1031.4 ± 4.6) 1071.6 ± 0.5 1052.6 ± 0.3 (1017.6 ± 4.8) 1085.6 ± 0.3 1084.9 ± 0.3 1069.7 ± 0.3 1061.6 ± 0.2 1169.4 ± 0.6 (1149.7 ± 1.1) 1153.3 ± 0.3 1129.0 ± 0.7

Exp VPL (VL ,˚ A3 ) 1256(303K) 72 1254(300K) 64,d

1272.27(320K) 64,d

1101(303K) 72,75 1096.5(300K) 64,e

1112.5(320K) 64,e

1175(308K) 77 1180(308K) 77

fields except for GROMOS give values close to the experimentally deduced area per lipid and bilayer thickness. However, a more direct comparison to experiment is required and will be provided in the next Section. Interestingly, for PC lipids the standard (noncurved) area per lipid is on average closer to experimental values while for PE lipids the curved area per lipid agrees better with experiments. A comparison of curved and standard areas per lipid suggests a lower bending modulus for PC lipids than for PE lipids; taking curvature into account the area per lipid increased by 5% and 3.7% for PC and PE lipids, respectively.

Significant deviations of the APL were observed for DMPC at 300 K for both the CHARMM36 and the Lipid14 force fields, as well as for POPE at 308 K for the GROMOS54a7 force field. These membranes underwent a phase transition to the gel or the ripple phase. As shown 12


Page 13 of 44




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 44

experimental values obtained from different approaches vary in the range of about 0.6%. 75 The volume per lipid was underestimated for all lipids by all force fields. Overall, the agreement with experiment is best for CHARMM36 (2.3% deviation for PC lipids and 0.7% for PE), closely followed by GROMOS54a7 (2.4% underestimation). The lipid volumes obtained by Slipids underestimate the experimental values by 3.8% (PC) and 2.1% (PE), and for Lipid14 by 4.6% (PC) and 4.1% (PE). Table 4: Isothermal area compressibility modulus (IACM) and bilayer thickness. R-squared values for linear intrapolation of experimental data are d 0.97 and e 0.99. Results for bilayers in a wrong phase are given in brackets. Systems

Force Field

POPC 300K


POPC 320K

DMPC 300K

DMPC 320K

POPE 308K

IACM [KA ,mNm−1 ] 352.8 ± 99.7 861.6 ± 271.9 395.9 ± 99.6 480.1 ± 172.0 280.7 ± 69.2 791.5 ± 172.1 271.1 ± 78.1 369.0 ± 115.8 (512.5 ± 140.3) 1013.6 ± 236.5 346.0 ± 157.8 (708.8 ± 281.2) 287.6 ± 58.5 840.0 ± 161.2 314.8 ± 49.0 351.0 ± 111.3 305.2 ± 83.4 (2069.7 ± 449.8) 361.4 ± 79.5 495.2 ± 110.4

Exp IACM [KA ,mNm−1 ] 180-330(298K) 81

234±23(294K) 145(300K) 83

84

233(303K) 83

Thickness [nm] 3.93 ± 0.02 3.90 ± 0.01 3.78 ± 0.02 3.83 ± 0.02 3.88 ± 0.01 3.84 ± 0.02 3.71 ± 0.02 3.76 ± 0.01 (3.88 ± 0.04) 3.54 ± 0.02 3.48 ± 0.02 (3.81 ± 0.06) 3.58 ± 0.01 3.45 ± 0.01 3.40 ± 0.01 3.47 ± 0.01 4.33 ± 0.03 (4.46 ± 0.02) 4.11 ± 0.02 4.32 ± 0.02

Exp Thickness [nm] 3.93(300K) 64,d

3.82(320K) 64,d

3.57 82 (300K) 83 3.70(300K) 64,e 3.60(303K) 75 3.53(303K) 72 3.53(320K) 64,e

3.95(308K) 77 4.05(308K) 76

In agreement with the reduced areas per lipid, the POPE bilayer thickness was slightly overestimated by all investigated force fields. The POPE thickness is best reproduced by Slipids (2.8% deviation to experiment) and significantly overestimated by GROMOS (by 11.5%), which is, however, not surprising due to the observed gel phase formation. CHARMM36 and Lipid14 both overestimate the POPE thickness by about 8%. The PC thickness (both POPC and DMPC, both temperatures, ripple phase system excluded from

14


Page 15 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


comparison) is best reproduced by CHARMM36 (1% deviation), followed by GROMOS54a7 (1.5%) and Lipid14 (1.9%) while Slipids understimates the PC thickness by 3.7%. Although the span of the experimental values for the isothermal area compressibility is large, all force fields overestimated this observable. The isothermal area compressibility modulus obtained for the GROMOS54a7 united-atom force field amounts to up to 200% of the values obtained by the all-atom force fields. A similar effect was observed by Piggot et al. 30 for a comparison of different GROMOS force fields to CHARMM36. Interestingly, the united-atom Berger force field yielded results more similar to CHARMM36. 30 Thus, the increased compressibility doesn’t result from the united-atom nature itself. Since both the Berger and the GROMOS54a7 force field carry the charges derived by Chiu et al., 85 the observed differences in the compressibility are likely coupled to differences in the van der Waals parameters.

Electron densities and scattering form factors Electron density profiles along the bilayer normal contain more detailed structural information than the bilayer thickness and the area per lipid. Unfortunately, the electron density profiles can not be obtained directly from experiment but require a fit (e.g. the SDP model 63 ). Therefore, the electron and neutron density profiles (EDP and NDP, respectively) obtained from simulations were transformed into the Fourier space and compared directly with the experimental X-ray and neutron scattering form factors. The physical meaning of form factors is, however, difficult to read. 3 summarizes both all electron density profiles and X-ray and neutron scattering structure form factors. Overall, a good agreement was obtained between simulated and experimental electron density profiles for systems staying in the liquid-crystalline membrane phase. Interestingly, the SDP model suggests a decreased electron density in the middle of the DMPC membrane at an increased temperature (320 K), which points to decreased interactions between the individual layers and probably even to a rise of the so-called free volume in the center of 15




Page 16 of 44

Page 17 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the bilayer. This increased drop in the center of the EDP at temperatures well above the melting temperature of the bilayer may, however, as well result from the method used to deduce the EDP from the underlying scattering data, 86 that uses gel phase structure as a reference. This reduction in the EDP was not reproduced by any of the studied force fields. A direct comparison of the form factors (see below) is therefore required. Interestingly, from the EDPs it can be seen that water (dashed lines) enters least deep into the lipid bilayer for the GROMOS54a7 simulations, which also correlates with the reduced solvation of the phosphate group and the decreased water flux across the bilayer in this force field (see below, 5). The agreement between the all-atom and the experimental neutron scattering data is best for CHARMM36 and Slipids force fields. Lipid14 slightly underestimates and GROMOS54a7 overestimates the neutron scattering form factors for PC lipids. Moreover, the neutron scattering form factors of all force fields decay (slightly) too fast for POPE. The X-ray scattering form factors were not perfectly reproduced by any force field. The best agreement for both the lobes’ heights and positions was observed for the Slipids force field with the exception of POPE, where the height of the first lobe is underestimated. The Lipid14 force field performed always well as for the height and position of the first lobe if the bilayer was in liquid-crystalline phase. With the exception of the POPE simulation, the reproduction of the second lobe is also satisfactory. The CHARMM36 force field shows a small shift of the maximum of the first lobe to smaller q values, consequently also the other lobes are slightly shifted in the same direction. The GROMOS54a7 performed comparable to the other force fields for the first and second lobe. Interestingly, the third lobe was more pronounced within the GROMOS54a7 force field. In general, all force fields showed significant deviations for POPC at 320 K as well as for POPE at 308 K.

17




Page 18 of 44

Page 19 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Acyl chain order parameters for POPC bilayers compared to experimental values. 9,90,91

Lipid acyl chain protrusion Another interesting property of lipid acyl chains is their protrusion probability from the hydrophobic membrane core to the polar lipid headgroup region. Lipid acyl chain protrusion was observed both in experiment 41 and in MD simulations before. 43,44,92 It was proposed to be connected to phenomenons like membrane fusion 40 as protrusion was shown to be enhanced in the presence of fusion peptides. 43,44 Protruding lipids are thought to connect the two opposing membranes 93 and thus lower the energy barrier of the transition state. 42 5 summarizes the lipid chain protrusion properties, namely, the lipid protrusion probability (LPP), the median of the lipid protrusion time (MPT), and the lipid protrusion efficiency (LPE). Extreme protrusion patterns were observed for lipid bilayers in the gel phase (GROMOS54a7 POPE, virtually no lipid acyl chain protrusion) and in the ripple phase of DMPC at 300 K using the Lipid14 force field. The latter system showed an enhanced protrusion efficiency. Additionally, lipid protrusion was seen to strongly depend on the force field: The protrusion probability is overall lowest within the GROMOS54a7 force field, connected probably to its drastically reduced diffusion (see below). 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Lipids using CHARMM36 parameters protrude more often, followed by lipids in Slipids. The Lipid14 force field resulted in a significantly enhanced protrusion rate as compared to the other force fields. The average protrusion time was generally below nanoseconds. Only the GROMOS54a7 force field led to significantly inreased protrusion times for all lipid types studied. The protrusion efficiency (estimated by multiplication of lipid protrusion probability and median lipid protrusion time) was 2- to 10-fold larger for the Lipid14 force field as compared to the other studied force fields. Thus, the Lipid14 force field may show enhanced fusogenic activity as compared to the other force fields. Overall, protrusion is decreased for PE lipids as compared to PC lipids and the number of lipid chain protrusions increases with inreasing temperature. Interestingly, a decreased protrusion efficiency was observed for POPE described by the CHARMM36 force field. This decreased protrusion efficiency is a result of both a lower protrusion probability and extremely short protrusion times of the protrusion events. For layer specific extended protrusion results including also the number of protrusion events and maximal protrusion times see Figure S4 and Table S2 in the Supplementary Information.

Lipid headgroup hydration As a measure of lipid hydration, the average numbers of coordinated water oxygens around the lipid phosphorus atoms as well as the radius of the first solvation shell (water oxygen atoms around phosphorus atoms) are summarized in 5. The largest solvations were observed within the CHARMM36 and the Lipid14 force fields, while the GROMOS54a7 force field showed a significantly decreased solvation with 1-2 water molecules less within the first solvation shell. As to be expected, the water coordination number correlates with the area per lipid for the different force fields. In contrast to observations of Ohta-Iino et al. 94 who studied membrane dehydration and lipid protrusion between two tightly spaced membranes, no dependence between the water 20


Page 20 of 44

Page 21 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5: Lipid protrusion characteristics, namely lipid protrusion probability (LPP), median lipid protrusion time (MPT) with the interquartile range (IQR), and lipid protrusion efficiency (LPE), defined as the product of LPP and MPT. Moreover, lipid hydration calculated as a water coordination number (WCN) to a phosphorus atom in its first solvation shell (FSS) and its radius are included. Values for membranes in a wrong phase are given in brackets. System

Force Field

POPC 300K


POPC 320K

DMPC 300K

DMPC 320K

POPE 308K

LPP [µs−1 ] 1.67 0.93 2.55 9.60 2.94 1.82 5.84 15.81 (1.26) 1.02 2.31 (9.95) 2.81 1.92 4.97 13.56 0.35 (0.04) 1.17 7.50

MPT [ns] 0.39 0.64 0.29 0.25 0.19 0.23 0.14 0.18 (0.45) 0.81 0.28 (0.44) 0.13 0.36 0.24 0.17 0.05 (0.68) 0.31 0.27

IQR [ns] 1.01 1.74 1.03 0.67 0.42 1.27 0.49 0.47 (1.87) 1.72 0.85 (1.27) 0.41 1.20 0.55 0.52 0.46 (0.66) 0.74 0.74

LPE [Number/lipid×10−3 ] 0.65 0.60 0.74 2.40 0.56 0.42 0.82 2.85 (0.57) 0.83 0.65 (4.38) 0.37 0.69 1.19 2.31 0.02 (0.03) 0.36 2.02

WCN [Number] 6.20±0.02 4.27±0.02 5.26±0.02 5.93±0.01 6.21±0.01 4.36±0.02 5.34±0.02 5.95±0.01 (5.88±0.03) 4.28±0.03 5.23±0.02 (5.65±0.04) 6.07±0.01 4.31±0.01 5.22±0.01 5.86±0.01 6.09±0.02 (4.57±0.03) 4.79±0.05 5.27±0.03

FSS Radius [nm] 0.446 0.456 0.440 0.434 0.450 0.460 0.444 0.436 (0.448) 0.456 0.442 (0.434) 0.450 0.460 0.442 0.436 0.450 (0.458) 0.442 0.432

coordination number and the lipid protrusion efficiency was obtained in single membrane simulations performed here.

Lipid headgroup orientation A detailed comparison of order parameters for carbon atoms in lipid headgroups was performed elsewhere. 34 Therefore, we focused here on the

31

P-13 C dipolar coupling (6) and the

orientation of the P-N vectors relative to the membrane normal (6). In all force fields, the replacement of the bulky choline group by an amine group in PE lipids allowed for larger P-N tilt angles. For the united-atom GROMOS54a7 force field the average P-N vector was almost parallel

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 44

Table 6: Average P-N vector tilt angle relative to the membrane normal. System DMPC 300K DMPC 320K POPC 300K POPC 320K POPE 308K

CHARMM36 (71.7) 70.9 69.1 70.1 77.5

GROMOS54a7 83.3 83.7 82.4 82.7 (92.6)

Lipid14 (69.6) 69.8 69.3 70.0 76.4

Slipids 68.3 69.3 68.3 69.1 76.1

Exp 61(313 K) 19 62(308 K) 15

to the bilayer plane for both PC and PE headgroups. The distribution of the P-N vectors in the gel-phase POPE bilayer (Figure S5) revealed that about 60% of all POPE headgroups are tilted by more than 90%, i.e. the amine groups are more deeply inserted in the membrane than the phosphate groups. For the all-atom force fields, the average P-N vector orientation of the PC headgroups relative to the membrane normal varied only by a few degrees (minimum 68.3◦ , maximum 71.7◦ ), independent of the lipid acyl chains and, interestingly, also of the lipid phase. This indifference of the headgroup orientation to the lipid chain is in agreement with experimental observations, which predicted the average P-N vector to be 61◦ and 62◦ for both DMPC at 313 K 19 and POPC at 308 K, 15 respectively. Compared to these experimental estimates, the orientation of the PC headgroup in all studied force fields was by 7-8◦ too strongly inclined towards the lipid bilayer plane. However, experiments do not allow to directly access the orientation of the P-N vector. The above values were calculated from a model fulfilling experimental observables. A more detailed comparison between experiment and simulations is achieved by comparison of the dipolar splittings between phosphorus and headgroup carbon atoms. The experimental data available for DMPC at 313 K 68 are shown in 6 together with results from the MD simulations at 320 K. Only small differences between the individual force fields for the predicted 31 P-13 C dipolar coupling were observed. The general pattern roughly agrees with the experimental one. A direct comparison of the absolute values is, however, hampered due to the difference in temperature of 7 K between experiment and simulation. The largest deviations from experiment were seen for carbons close to the phosphate group, namely for C11, C1 and

22


Page 23 of 44





Page 24 of 44

Page 25 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


coefficient is significantly lowered. Overall, the Slipids and CHARMM36 force fields showed the best agreement with experiment for both different lipids and different temperatures. A comparison to values reported by Piggot et al. 30 shows that the diffusion coefficient for POPC using CHARMM36 (at 298 K in case of Piggot et al., here 300 K) is reduced by about 50% here. This discrepancy is probably related to the increased number of lipids in this study, an effect ascribed to a reduction in the lipid correlated motions across periodic boundaries in larger system. 30 In all simulations using GROMOS54a7 the lipid diffusion was systematically underestimated. As observed earlier by Poger et al., 103 lipid diffusion coefficients are systematically lower in simulations using a reaction-field (RF) correction as compared to PME summation. An additional simulation using GROMOS54a7 in combination with PME summation for DMPC at 320 K (shown as a red cross in 7) resulted in a diffusion coefficient that is slightly increased as compared to the RF simulation. However, the absolute value is still 5-8 times below those obtained by all-atom force fields or reported from experiment. The same phenomenon was also observed by Piggot et al. 30 and was remedied by reducing the cut-off for nonbonded interactions to 1 nm. However, this solution is not compatible with the usage of a 1.4 nm cutoff in the GROMOS94 force field family. 103

Melting temperature The membrane phase transition temperatures in the investigated force fields was studied from heating simulations starting from a membrane in its gel phase. Because of the tendency of DMPC bilayers to form a ripple phase close to the melting temperature and due to the low melting temperature of POPC (271 K), the analysis was restricted to a DPPC bilayer as a test system. Snapshots of the DPPC bilayer in gel states for the different force fields are shown in 8 together with the area per lipid as a function of temperature from heating simulations starting at 300 K at a heating rate of 1 K/ns. For lower heating rates a decreased phase transition 25




Page 26 of 44

Page 27 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 8: Melting characteristics. The steepness reflects the speed of evolution in area per lipid at Tm , the transition width is reflected by the full width at half maximum (FWHM). Force field GROMOS54a7 CHARMM36 Lipid14 Slipids

Tm [K] 335.06±0.57 341.21±0.50 343.25±0.32 329.08±0.31

Steepness [ nm2 /K ] 0.0067±0.0008 0.0252±0.0015 0.0317±0.0024 0.0186±0.0021

FWHM [K] 14.45±3.03 4.97±0.29 5.08±0.39 6.18±0.73

was still partially disrupted. Accordingly, the melting curves for GROMOS54a7 started at a larger area per lipid. The Slipids gel phase DPPC bilayer exhibited a slightly larger area per lipid as compared to CHARMM36 ad Lipid14. In agreement with the observed phase transition of DMPC at 300 K using the CHARMM36 or Lipid14 force field to a ripple phase (see 2), the melting temperature of DPPC using these force fields is by 12-14◦ increased as compared to the Slipids system (329 K). The average melting temperature for GROMOS54a7 (335 K) was slightly higher than for Slipids. In order to test the influence of the larger area per lipid on the melting behavior, the Slipids DPPC system at 300 K was converted to the Lipid14 force field and heated (blue line in the Lipid14 plot). The main phase transition appears to be unaffected by this moderate change in starting configuration. The decreased melting temperature for Slipids is therefore due to the force field and not due to the larger area per lipid of the starting structure. The sharpest phase transition (steepest rise of APL upon phase transition) was observed for the Lipid14 force field (see 8). The GROMOS54a7 force field showed a very broad transition regime between the gel phase and the liquid crystalline phase.

Performance Different force fields may show a different computational efficiency based on different degrees of fine-graining (united-atom vs. all-atom) and different simulation parameters employed. The relative performance of the studied force fields is collected in 9. The long-range electrostatics (PME) and the cut-off radius for van der Waals and (short27




Page 28 of 44

Page 29 of 44

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


40 threads were run in parallel on one node). The performance loss is moderate for all force fields for the number of nodes investigated (2-6, corresponding to 40-120 compute cores).

Conclusions All studied force fields, namely the all-atom CHARMM36, Lipid14, and Slipids, as well as the united-atom GROMOS54a7 force field, were found to satisfactorily describe the overall structural characteristics of phospholipid bilayers. Also the lipid diffusion and water permeability were well described by the three all-atom force fields. However, also deficiencies in reproducing different observables are still seen for all lipid force fields. The GROMOS54a7 force field was e.g. found to underestimate the lipid self-diffusion and showed some weakness in describing the melting characteristics of membranes and is therefore less suitable for studies where dynamic lipid behavior is expected to play an important role. A summary of the performance of the different lipid force fields is graphically illustrated in 10. In detail, Slipids is very well suited for membrane studies closely above the melting temperature. On the other hand, Lipid14 and CHARMM36 showed well ordered membrane structures for the gel phase and a highly cooperative transition during melting. The acyl chain order and X-ray and neutron scattering form factors were best reproduced by Slipids and Lipid14. CHARMM36 yields improved values for the lipid volume and bilayer thickness as compared to the other all-atom force fields and describes most realistically the lipid diffusion. The computationally most efficient all-atom force field is Lipid14. We would like to note, that we did not observe any drawbacks or artifacts resulting from the usage of the short radii for nonbonded interactions (1.0 nm) using Lipid14. On the contrary, the small cutoff increased the computational efficiency to a value comparable to the unitedatom GROMOS54a7 force field. The combination of this computational efficiency and the significantly enhanced lipid acyl chain protrusion probability of Lipid14 renders this force field the most efficient choice for membrane fusion studies at atomistic resolution.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

"JFSRR29

TSURUVD6$5

!"

!#

!"

0123

4153

;&,$?&@AB/CC 8143 F@

A Critical Comparison of Biomembrane Force Fields: Structure and

Recommend Documents