Lipid Bilayer Simulations Force fields Simulation and Analysis

Lipid Bilayer Simulations: Force fields, Simulation and Analysis Jeffery B. Klauda Model Yeast Membrane Chemical Structure of Lipids

Lipids w Complex biomolecules • Contain a fatty acid chains and head group w Classified into 8 categories 1 Modified (Fig. 1)1 1 Fahy Fatty acyls Glycerolipids Glycerolphospholipids Sphingolipids Sterol Lipids Phenol lipids Saccharolipids Polyketides et al. J. Lipid. Res. 46: 839 (2005).

Glycerophospholipids w Some Common Subclasses of GP lipids 1 Phosphocholines Phosphonocholines Phosphoethanolamines Phosphonoethanolamines Phosphoserines Phosphoglycerols Phosphoglycerophosphates Phosphoinsitolmonophosphates 1 Fahy et al. J. Lipid. Res. 46: 839 (2005). Phosphates (Modified Fig. 4)1

Membranes in Single Cell Organisms Periplasm Cytoplasmic Membrane Channel Proteins Lipid/Cholesterol Bilayer Membrane Proteins Cytoplasm E. coli • Plasma membrane 1 contains many constituents • Membranes are located throughout the cell interior Cell Membranes 2 1 Fig. 2 Fig. 1 b from Engelman, D. M. Nature. 438: 578 (2005). 1 a from Mc. Mahon, H. T. et al. Nature. 438: 590 (2005).

Diversity of Lipid Types in Organisms w Yeast (Saccharomyces cerevisiae)1 • Mixture of fully saturated and unsaturated chains • Mixture of charged and zwitterionic head groups and typically 10 -20% sterols • Compositions depend on strain of yeast w Chlamydia (chlamydia trachomatis)2 • Exists in two forms reticulate body (metabolically active) and elementary body (infectious) • Bacterial membranes contain various chain types including branched chains (10 -20%) • Primarily zwitterionic and phosphatidylinositol head groups • 20 -30% sterols w E. coli (Escherichia coli)3 • Mixture of fully saturated and unsaturated chains • Fatty acid chains can contain cyclic moieties (cyclopropane) • Zwitterionic (~80% PE) and anionic (~20 %PG) head groups • Limited to no sterols 1 Daum 2 Wylie, J. L. et al. J. Bact. 179: 7233 (1997). G. et al. Yeast. 15: 601 (1999). A. & G. Larsson. Microbial Cell Factories. 3: 9 (2004). 3 Shokri,

Membrane Composition within Cell w Distribution of phospholipids (PL) vs. sterols 1 · Mammals in dark blue and yeast in light blue · Plasma membrane (PM) contains a significant amount of sterol (largest of all organelles) · Mammalian PM contain more sterol than yeast · Endoplasmic reticulum (ER) manufactures sterol, but levels are low · Large diversity of phospholipids between mammals and yeast and within a cell 1 van Meer, G. et al. Nature Rev. Mol. Cell. Bio. 9: 112 (2008).

Force Fields w Biomolecular Force Field (CHARMM) · Many terms to describe intra- and intermolecular interactions w All-atom Lipid Force Fields • CHARMM Family: CHARMM 27 r and CHARMM 36 (C 27 r 1 and C 362) • AMBER Family: GAFFlipid 3 and Lipid 144 • Stockholm Lipids (Amber-compatible): Slipid 5 1 Klauda, 2 Klauda, 3 Dickson J. B. et al. JPCB. 109: 5300 (2005). et al. Soft Matter. 8: 9617 (2012). 5 Jämbeck & Lyubartsev. JPCB. 116: 3164 (2012). 4 Dickson J. B. et al. JPCB. 114: 7830 (2010). et al. J. Chem. Theory Comput. 10: 865 (2014).

AMBER Lipids w Summary of Lipid 14 FF 1 Results (NPT Ensemble) Surface Area/lipid [Å2] MD DPPC DMPC DLPC DOPC 62. 0± 0. 3 59. 7± 0. 7 63. 0± 0. 2 69. 0± 0. 3 POPC POPE 65. 6± 0. 5 55. 5± 0. 2 Exp 63. 0± 1. 0 60. 6± 0. 5 63. 2± 0. 5 67. 4± 1. 0 68. 3± 1. 5 ~60 · Generally good agreement with experiment (slight tendency to underestimate) Deuterium Order Parameters · Overall excellent agreement with NMR SCDs · POPE SCDs of the saturated chain are somewhat high, which may indicate that the SA/lipid is too low · Decent splitting for the C 2 position (Fig. 71) · Unclear if the head group order parameters are in agreement with experiment • Procedure follows typical rules for AMBER FF optimization (RESP charges in gas phase) • Should only be used with the AMBER family of FF 1 Dickson et al. J. Chem. Theory Comput. 10: 865 (2014).

Stockholm Lipids (Slipids) w Summary of Slipids 1 -3 Results (NPT Ensemble) Surface Area/lipid [Å2] MD DPPC DMPC DLPC DOPC 62. 6± 0. 5 60. 8± 0. 5 62. 4± 0. 4 68. 0± 0. 5 POPC POPE 64. 6± 0. 4 56. 3± 0. 4 Exp 63. 0± 1. 0 60. 6± 0. 5 63. 2± 0. 5 67. 4± 1. 0 68. 3± 1. 5 ~60 · Generally good agreement with experiment (slight tendency to underestimate) Deuterium Order Parameters (Fig. 51) (Fig. 22) · Overall excellent agreement with NMR SCDs · Better POPE SCDs compared to Lipid 14 · Decent splitting for the C 2 position · Unclear if the head group order parameters are in agreement with experiment • Procedure similar to AMBER FF optimization (RESP charges in gas phase) • Extensions to PS, PG and SM lipids 3 1 Jämbeck 3 Jämbeck & Lyubartsev. JPCB. 116: 3164 (2012). 2 Jämbeck & Lyubartsev. J. Chem. Theory Comput. 8: 2938 (2012). & Lyubartsev. J. Chem. Theory Comput. 9: 774 (2013).

Force Fields Continued w Biomolecular Force Field (CHARMM) · Many terms to describe intra- and intermolecular interactions w United Atom/Coarse-grained Lipid Force Fields • United atom: C 27 -UA(acyl)1, C 36 -UA 2 and GROMOS 3 • Coarse-grained: MARTINI 4 and Shinoda/De. Vane/Klein 5 1 Henin, 2 Lee, Tran, Allsopp, Lim, Henin & Klauda. JPCB 118: 547 (2014). Shinoda & Klein. JPCB. 112: 7008 (2008). 4 Marrink et al. JPCB. 111: 7812 (2007). , O. et al. BJ. 72: 2002 (1997). 5 Shinoda, De. Vane, & Klein. JPCB. 114: 6836 (2010). 4 Berger

Issues with the CHARMM 27 r FF w Surface Tension • To maintain good agreement with density profiles and SCD, NPAT simulations at the experimental area are needed for MD simulations with C 27 r • Finite size effects may result in a non-zero surface tension, 1 but C 27 r values are too high 2 Surface Tension in dyn/cm LR LJ No LR-LJ Exp. Estimate DPPC bilayer (64 Å2/lipid, 323 K) 19. 7 16. 8 ~0 -5 DMPC bilayer (60. 7 Å2/lipid, 303 K) 19. 8 -- ~0 -5 w Freezing or Phase Change with NPT · Freezing of aliphatic chains at T > Tb · Issue with lipids that have 1 -2 fully saturated chains · Problematic when surface areas are not available for lipids and their mixtures 1 Klauda, 2 Klauda, J. B. et al. BJ. 90: 2796 (2006). J. B. et al. JPCB. 111: 4393 (2007).

Modification of CHARMM Charges w Charge/LJ Modification • Looked at small molecules and DPPC bilayer charges using semi-empirical AM 1 · Increase in polarization occurred going from the gas phase to realistic bilayer · Therefore, increasing the lipid charges in the glycerol region is justified • Adjusted charges/LJ Dipole moment of methylacetate (debye) Dipole QM C 27 r C 36 X/Y Ratio 1. 48 -7. 83 1. 52 Total 1. 65 2. 40 1. 52 · Adjustments are supported by AM 1 on the bilayer, small molecule dipoles and watermolecule interactions. · Small adjustments on the carbonyl carbon LJ parameters with the ester oxygen taken from previous optimizations 1 1 Vorobyov, I, et al. J. Chem. Theory and Comp. 3: 1120 (2007).

Dihedral Modifications Fits to QM of small molecules QM of bilayers (Alex Mac. Kerell) w Small Molecule Models of DPPC 1 Klauda et al. JPCB. 114: 7830 (2010).

Dihedral Modifications: CHARMM 36 w Glycerol FF Adjustments • Adjust the g 1 torsion q 2 q 4 b 1 g 1 MP 2/6 -31 g(d): 648 Energy Points • Issues fitting the q 4 and b 1 torsions kcal/mol 1 Klauda et al. JPCB. 114: 7830 (2010).

Empirical Fits of Torsions (C 36) w DPPC SCD Targets • MD simulations of the DPPC bilayer with an intermediate FF were used to empirically fit q 2, q 4, and b 1 torsions. • Populations of trans and gauche conformations of these torsions were optimized G+ T G- q 2 18% 36% 45% q 4 66% 3% 31% b 1 56% 43% 1% · The torsional potential was adjusted to bound the PMFs based on these fits and the optimal set was chosen. 1 Klauda et al. JPCB. 114: 7830 (2010).

Empirical Fits of Torsions (C 36) w Torsional surface scans from 20 ns MD simulations 1 Klauda et al. JPCB. 114: 7830 (2010).

DPPC Bilayer and C 36 w Deuterium Order Parameters (SCD): NPAT/NPT 1 vs. Experiment 2 NPAT A=64Å2 NPT · Excellent agreement with experiment and fairly independent of the ensemble. 1 Klauda, J. B. et al. JPCB. 114: 7830 (2010). 2 Seelig, A. & J. Seelig. Biochem. 13: 4839 (1974).

DPPC Bilayer w Density Profiles & Form Factors Compared to Experiment 1 Aexp=63± 1Å2 · Good agreement with the experimental form factors, F(q) · The methyl & methylene density is improved · NPT captures the overall and component densities correctly 1 Kučerka, N. et al. BJ. 95: 2356 (2008).

CHARMM 36 Lipids w Initial Parameterization with PC & PE lipids Surface Area/lipid [Å2] DPPCa DMPCb DLPCb MD 62. 9± 0. 3 60. 8± 0. 2 Exp 63. 0± 1. 0 60. 6± 0. 5 DOPCb POPEc 64. 4± 0. 3 69. 0± 0. 3 64. 7± 0. 2 59. 2± 0. 3 63. 2± 0. 5 67. 4± 1. 0 68. 3± 1. 5 ~60 w Additional Lipids • Lipids with polyunsaturated chains 2 DAPC • Branched and cyclic-containing chains (important for certain bacteria)3, 4 • Sterols (cholesterol, oxysterols, ergosterol)5 • Various published parameters: PS, PG, PA and Cardiolipin • Other lipid parameters on the way: PI, PIP, SM, and CER 1 Klauda 2 Klauda 3 Lim et al. JPCB. 114: 7830 (2010). & Klauda. BBA: Biomemb. 1808: 323 (2011). 5 Lim et al. JPCB. 116: 203 (2012). 4 Pandit et al. JPCB. 116: 9424 (2012). & Klauda. BBA: Biomemb. 1818: 1818 (2012). a 323 K b 303 K c 310 K

CHARMM-GUI w CHARMM-GUI. org – Membrane Builder 1, 2 Dr. Im (KU) • Allows for easy building of lipid membranes • Select from 140+ lipids and any mixture from these lipids • Builds membranes and provides rigorously tested equilibration inputs for CHARMM and NAMD simulations • Membrane proteins can be easily incorporated into the bilayer • Freely available to any researcher 1 Jo, 2 Jo, Kim, Iyer & Im. J. Comput. Chem. 29: 1859 (2008). Lim, Klauda & Im. Biophys. J. 97: 50 (2009).

CHARMM-GUI w CHARMM-GUI. org – Membrane Builder 1, 2 • Can easily build heterogeneous bilayers • Specify water hydration in three ways (defaults are safe for fully hydrated bilayers) • Can choose ratio or number of lipids for each leaflet • Reported surface area per lipid is based on simulations with a pure membrane • Further steps ask for ion concentration, ring penetration checks, ensemble and temperature • At the end (or during the process) you can download the files in. tgz format (all files needed to simulate bilayer) 1 Jo, 2 Jo, Kim, Iyer & Im. J. Comput. Chem. 29: 1859 (2008). Lim, Klauda & Im. Biophys. J. 97: 50 (2009).

CHARMM-GUI w Output Initial Structure of Bilayer • Water is initially away from bilayer (will quickly fill in the vacuum space). • Lipid head groups are aligned to a specific z-position based on prefered location in the bilayer • Chains can tangle and careful equilibration is required w Restraints During Equilibration • Water away from hydrophobic core • Head group and tails to appropriate regions • Double bonds in their respective cis or trans conformation • Ring conformations (chair & upright for PI lipids) 1 Jo, 2 Jo, Kim, Iyer & Im. J. Comput. Chem. 29: 1859 (2008). Lim, Klauda & Im. Biophys. J. 97: 50 (2009).

MD Simulations of Membranes w Caveats of CHARMM-GUI with membranes • Membrane surface area/volume · Primarily based on SA from pure lipid bilayers with C 36 force field at 303 K · Some lipids have high gel transition temperatures >303 K and values are based on higher temps · This can result in poor initial guess for mixed lipid systems, especially with sterols · If the SA is known or can be estimated a priori then this is preferred • Membrane equilibration · Although we have tested this extensively there might be some issues · Pay careful attention to your bilayer lipids · Make sure all bonds are maintained after equilibration, otherwise results will be off · Building the membrane may cause chain overlap · Internal checks for ring penetration by chain (chain through cholesterol or amino acid rings) · If these exist, then you need to rebuild the system!

Simulation Snapshot ERG, YOPS, DYPC and water

Multilayer System/Periodic Boundary Conditions

ST-Analyzer w Web-based Interface for Simulation Trajectory Analysis 1 Dr. Im (KU) • Allows for easy collection of data on membranes and proteins • Can be setup to on a workstation or a cluster environment with batch submission of analysis 1 Jeong et al. J. Comput. Chem. 35: 957 (2014).

Membrane Area per Lipid w Equilibrated? Thermal Equilibration NPT-Production • Things to consider with membrane equilibration · Possible transient stability in volume/surface area · Must run for long periods of time: 10 -30 ns for simple single lipids and 50 -300 ns (or longer) for complex mixtures (General rules of thumb without phase transitions) · Current run suggest longer times (beyond 20 ns) are needed

Membrane Area per Lipid: Examples • Equilibration is slower during changes in phase (La to gel- like phase) • 100 ns or greater can be required to obtain a fully equilibrated bilayer even for single lipids DPPC at 200 ns (303 K) z

Lipid Bilayer Structure: Simulation w Molecular Dynamics • Simulations can easily obtain density profiles Bulk Water Headgroup Dz=0. 1Å Count number electrons/bin and average 1 Jo, 2 Jo, Kim, Iyer & Im. J. Comput. Chem. 29: 1859 (2008). Lim, Klauda & Im. Biophys. J. 97: 50 (2009). SM=Structural Model Hydrophobic/Chain

Lipid Bilayer Structure: Experiment w Form Factors F(q) • F(q) is transformed into real space to get structural properties F(q) Structural Models EDP, A Fourier Reconstruction Only Total EDP & Fourier Wiggles HB Fit to Exp. F(q) for the DMPC Bilayer 1 1 Kučerka, N. et al. Biophys. J. 88: 2626 (2005).

Development of H 2 Structural Model w Density Profile • Component electron density used to guide model development Asim=60. 7 Å2 Black & Blue: Simulation Red: H 2 fit to density w New Hybrid Model (H 2)1 • Consists of five physical components 1 Klauda, J. B. et al. Biophys. J. 90: 2796 (2006). BC=water+choline CG=carbonyl-glycerol chol = choline

Comparing to X-ray/Neutron Scattering w Model Free Comparison • Form Factors (symmetric bilayers, where D-repeat spacing) fa(q): atomic form factors (depend on q for X-ray (not neutron data)) na(z): atomic number distribution (density of atoms of each type) r. W: scattering density of water (solvent) • Method to use and program · Calculate atomic densities (na(z)) (in CHARMM or ST-Analyzer) and use SIMto. EXP program 1 · Load in atomic density to SIMto. EXP program to get F(q) 1 Kučerka et al. J. Membr. Biol. 235: 43 (2010).

Examples for C 36 POPC Form Factors & Density Profiles 1 · Excellent agreement between experiment and MD simulation form factors. · Can easily obtain density profiles of groups within the bilayer 1 Zhuang, Makover, Im & Klauda. BBA-Biomemb. Submitted (2014).

Overview of Lipid Dynamics and Internal Structure w Range of Lipid Motions 0 Wobbling in a Cone Model 1 Bond Vibrations 1 ps Hydrogen Bonds 1 ns Internal Isomerization (C-H, P-H, etc. ) and Wobbling 3 10 ns Isomerization Axial Rotation Lipid Axial Rotation 3 Lateral Diffusion 2 1 ms Vesicle Rotation Wobbling w Internal Structure • Orientation of bonds • Angle of bond vectors with respect to bilayer normal w Methods to obtain these Quantities • NMR • Molecular dynamics 1 Pastor, 3 Klauda R. W. et al. Accounts. Chem. Res. 35: 438 (2002). et al. Biophys. J. 94: 3074 (2008). 2 Klauda et al. J. Chem. Phys. 125: 144710 (2006).

NMR Background w Nuclear Magnetic Resonance (NMR) • Magnetic nuclei (13 C/31 P) respond to an oscillating magnetic field • Spin-lattice relaxation rates (R 1) • Dipolar term: nuclear spin interaction between neighbors Spectral Density Reorientational Correlation Function 2 nd Order Legendre Polynomial Unit vector between P and its neighboring H

NMR Background w Nuclear Magnetic Resonance (NMR) • Magnetic nuclei (13 C/31 P) respond to an oscillating magnetic field • Spin-lattice relaxation rates (R 1) • Chemical Shift Anisotropy: on nucleus · Based on sold-state measurements on lipids 1 · Major principal axis 1 (s 33) is used to obtain the spectral density · The asymmetry in principal axis is accounted for by h • Field dependence · Dipolar contribution is important at low field · CSA is important at high field 1 Herzfel, J. et al. Biochem. 17: 2711 (1978).

Deuterium NMR w Deuterium NMR • Order parameters · qi is the angle of a C-D vector with the bilayer normal (usually the z axis) • Internal structure of lipids w How obtain this via MD Simulations • Calculate the C-H angle (MD simulations without deuterium) • Do this for every carbon • Simple trig calculation

Deuterium NMR: Examples SCD’s for POPC and DLPC 1 · Higher values indicate more order (lower disorder) · Double bond adds a kink to the chain and more disorder · SCDs depend on temperature and agree fairly well with experimental data 1 Zhuang, Makover, Im & Klauda. BBA-Biomemb. Submitted (2014).

Summary • There are many lipid types that can exist in biology and each has it own function to the cell • Lipid diversity in biology can vary between different head groups to chain types • Lipids from in vivo membranes are diverse between organisms and organelles with a single organism • There are several options for lipid force fields to run MD simulations (all-atom, united-atom and coarse-grained) • CHARMM 36 lipid force field has been parameterized and agrees well with a multitude of experiments (dynamical and structural) for all regions of the lipid • CHARMM-GUI allows for easy building of simple and complex membranes with/without proteins • ST-Analyzer allows for easy access and analysis of simulation trajectories from many different simulation program platforms • A key test for bilayer equilibration is the surface area per lipid • Structural (SCD and density profiles) and dynamical properties (diffusion and relaxation rates) can easily be obtained with proper analysis of MD simulations