Molecular Dynamics Simulations of Ceramide and Ceramide-Phosphatidylcholine Bilayers.

Subscriber access provided by Eastern Michigan University | Bruce T. Halle Library

Article

Molecular Dynamics Simulations of Ceramide and Ceramide-Phosphatidylcholine Bilayers Eric Wang, and Jeffery B. Klauda J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b08967 • Publication Date (Web): 10 Oct 2017 Downloaded from http://pubs.acs.org on October 16, 2017

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 38

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

Molecular Dynamics Simulations of Ceramide and Ceramide-Phosphatidylcholine Bilayers Eric Wang1 and Jeffery B. Klauda1,2* 1

Department of Chemical and Biomolecular Engineering, 2Biophysics Graduate Program University of Maryland, College Park, MD 20742, USA *Corresponding Author: [email protected] Abstract Recent studies in lipid raft formation and stratum corneum permeability have focused on

the role of ceramides (CER). In this study, we use the all-atom CHARMM36 (C36) force field to simulate

bilayers

using

N-palmitoyl-sphingosine

(CER16)

or

α-hydroxy-N-stearoyl

phytosphingosine (CER[AP]) in 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) or 1palmitoyl-2-oleoylphosphatidylcholine (POPC) which serve as general membrane models. Conditions are replicated from experimental studies for comparison purposes, and concentration (XCER) is varied to probe the effect of CER on these systems. Comparison with experiment based on deuterium order parameters and bilayer thickness demonstrates good agreement which supports further use of the C36 force field. CER concentration is shown to have a profound effect on nearly all membrane properties including surface area per lipid, chain order and tilt, area compressibility moduli, bilayer thickness, hydrogen bonding, and lipid clustering. Hydrogen bonding in particular can significantly affect other membrane properties and can even encourage transition into the gel phase. Despite CER’s tendency to condense the membrane, an expansion of CER lipids with increasing XCER is possible depending on how the balance between various hydrogen bond pairs and lipid clustering is perturbed. Based on gel phase transitions, support is given for phytosphingosine’s role as a hydrogen bond bridge between sphingosine ordered domains in the stratum corneum.

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

1.

Page 2 of 38

Introduction Sphingolipids are membrane lipids specific to eukaryotes with cellular roles including

stress resistance1, metabolism, and lipid raft formation2. Lipid rafts, ordered microdomains within the cell membrane enriched in sphingomyelin (SM) and cholesterol (Chol), have been implicated in a wide variety of cellular processes including signal transduction3 and nerve terminal function4. Additionally, rafts have been suggested to contain significant amounts of ceramides (CER) during signal transduction processes with functional CER-enriched macrodomains forming5-6. CER lipids make up the backbone of sphingolipids and are produced through the hydrolysis of SM by the enzyme sphingomyelinase7 which removes the phosphocholine headgroup. CER functions in cell signaling, acting as an important second messenger that affects membrane properties such as lipid phases and Chol concentration1, and the treatment of cancer cells with apoptotic agents increases CER production which has led to questions of its function as a tumor suppressor lipid8. Of particular interest is the high proportion of CER lipids in skin membranes, especially in the outer layer called the skin’s stratum corneum (SC)9. The SC is composed of corneocytes, dead cells containing keratin and water, within an intercellular lipid matrix10. The majority of skin penetration occurs through the lipid matrix and is closely related to the concentration and types of CER present. The involved CER types can be categorized into three groups based on the presence of either a sphingosine, phytosphingosine, or 6-hydroxysphingosine headgroup11. To properly understand the effect of CER on membranes, binary mixtures of CER and phosphatidylcholine (PC) lipids are often studied as a model since pure PC membranes function well as general membrane models. Experimental techniques such as NMR and neutron/X-ray scattering are essential to the accurate study of CER-containing membranes, but these are resource demanding and do not provide access to atomic-level detail. A complementary approach is the use of molecular dynamics (MD) simulations which has enjoyed a strong increase in accessibility in recent years as computational power has rapidly expanded12-14. Molecular dynamics studies of pure CER bilayers have noticed the effect of head group size15, force field parameters16, and hydroxylation17 on membrane properties. Small head group size is extremely influential and produces weaker peaks in electron density profiles compared to SM whereas varying force field parameters demonstrate consistent results with some discrepancies in hydrogen bonding. 2 ACS Paragon Plus Environment

Page 3 of 38

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 between sphingosine and phytosphingosine head groups, differing in the presence of a double bond and hydroxylation, revealed the formation of a distinct hydrogen bonding network in phytosphingosine CER. These results, particularly on the effect of head group type, raise questions about properties in binary mixtures. However, although CER/PC mixtures have been extensively studied by experimental methods, few MD simulations have been performed, and these are often run for relatively short timescales (≤ 100 ns)5,

18

. While pure PC systems

equilibrate relatively quickly, binary mixtures potentially relax much more slowly, especially those involving CER, which has a significantly different structure from PC lipids with its missing phosphocholine headgroup. These studies also do not study the effect of different CER hydroxylations on the membrane. Therefore, current research is lacking in studies that have run long enough to ensure equilibration and that probe the effect of CER structure. The purposes of this study are: (1) to understand how CER structure and concentration in PC bilayers influences membrane properties with a focus on barrier function in the SC; (2) to verify the accuracy of the CHARMM36 (C36) force field by comparing against previous experimental and computational results. The bilayer mixtures studied in this work are mixtures of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) with either N-palmitoyl-sphingosine (CER16) or α-hydroxy-N-stearoyl phytosphingosine (CER[AP]) and mixtures of 1-palmitoyl-2oleoylphosphatidylcholine (POPC) with CER16 at varying CER concentrations (XCER) and temperatures. MD simulations of single-component lipid bilayers of CER16 and CER[AP] are performed to compare with the binary mixtures with PC lipids. The simulation conditions were chosen to match compositions from previous experimental studies5, 8, 19 and various membrane properties were calculated. The results generally show good agreement with past results and elucidate how CER affects PC membranes.

2.

Methods

2.1. System setup and simulation protocol Previous experimental and computational studies have produced data for CER/PC mixtures which vary in bilayer composition and experimental conditions5,

8, 19

. The simulated

systems in our work were built to mimic these experimental systems and provide comparison with the C36 force field20-21, and additional systems were built with a constant XCER=0.60 to analyze the effect of CER concentration on bilayer properties. The specific lipid ratios and 3 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

Page 4 of 38

system conditions can be found in Table 1. All simulations contained 50 lipids per leaflet with 30 water molecules per lipid except for a single system containing 24 waters per lipid to probe the effect of lower hydration. Triplicate replicas of each system were generated in rectangular boxes and tetragonal cells (X = Y ≠ Z) using CHARMM-GUI Membrane Builder and subsequently equilibrated while slowly removing restraints using the standard Membrane Builder six-step process22-25. CER[AP] was not available in CHARMM-GUI Membrane Builder by default and was obtained by mutation of N-stearoylsphingosine. Production runs were run for 400-1500 ns such that all replicas had at least 150 ns of equilibrated data which was determined by examining plots of surface area per lipid (SA/lip) versus time. The final production time for each system can be found in Table 1, which totals 15.3 µs of simulation time for all replicas.


Page 5 of 38

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 1. a) Chemical structure of POPC, DMPC, CER16, and CER[AP]. Atom labels include terminal tail carbons, tilt vector carbons, RDF selection atoms, and component surface area representative atoms. b) Snapshot of DMPC and CER16 at XCER=0.20 at the end of the simulation. DMPC lipids are shown in cyan with colored headgroups, CER16 molecules are in yellow, and water is shown as purple dots. c) Snapshot of DMPC and CER16 at XCER=0.60 at the end of the simulation. Color scheme is same as in (b)



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 6 of 38

Table 1. System parameters and general bilayer properties. SA/lip is the average for all lipids (PC and CER). System

XCER

T (K)

# Water

Time (ns)

SA/lip (Å2)

KA (N/m)

DMPC/CER16 DMPC/CER16 DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16

0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60

310 310 310 310 303.15 303.15 303.15 330.15 330.15

3000 2400 3000 3000 3000 3000 3000 3000 3000

400 500 500 500 1500 400 500 400 400

55.3±0.2 54.9±0.1 43.2±0.2 43.6±0.3 46.3±0.4 43.9±0.4 45.3±0.6 65.2±0.1 54.2±0.1

0.26±0.01 0.22±0.01 1.45±0.15 1.50±0.10 1.26±0.19 1.83±0.05 2.46±0.16 0.25±0.01 0.27±0.03

MD simulations were carried with the CHARMM36 lipid force field20 with updates for sphingolipids21 together with the TIP3P water model26-27. The constant molecule number, pressure, and temperature (NPT) ensemble was used in NAMD28 with a constant pressure of 1 bar and with temperatures set as specified in Table 1. Constant pressure was maintained through the Nosé-Hoover Langevin-piston algorithm29-30 while a constant temperature was maintained through Langevin dynamics. The particle-mesh Ewald (PME) method31 was used to model longrange electrostatics with an interpolation order of 6 and direct space tolerance of 10-6, and electrostatics were calculated every step. The Lennard-Jones potential was used with a forceswitching function32 in the interval of 10 to 12 Å to describe van der Waals interactions. Bond lengths involving hydrogen were kept constant by the RATTLE algorithm33. The simulation time step was 2 fs and coordinates were saved every 5 ps. Visual Molecular Dynamics34 (VMD) was used to visualize bilayers and create snapshots, and gnuplot was used to create plots35. 2.2. Analysis Analysis was consistently performed using the last 150 ns. Based on SA/lip plots, all systems had reached equilibrium by the beginning of analysis. The data analyses presented in this work are overall SA/lip, deuterium order parameters (SCD), component SA/lip, area compressibility modulus (KA), gel phase tilt angle, electron density profiles (EDPs), bilayer thickness, hydrogen bonding (H bonding), two-dimensional radial distribution functions (2D RDFs), and lipid clustering.


Page 7 of 38

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


Most analyses were performed with CHARMM36-37, and averages with standard errors were calculated from triplicates. SCD was calculated by the formula,

= 〈 − 〉

(1)

where θ is the angle between the C-H vector and the bilayer normal. At each carbon, the average of all C-H vectors is calculated. Published quadrupolar splitting values were converted to order parameters by the formula,

∆ = | (cos ( ))|| |

(2)

where ∆vQ is the quadrupolar splitting, χQ is the quadrupolar coupling constant (167 kHz for

aliphatic C-2H bonds38), P2 is the second-order Legendre polynomial (3 − 1), and θ is the angle of the applied magnetic field with respect to the bilayer normal. The tilt angle distribution of the gel phase was calculated by determining the angle between the lipid tail vector (defined as the C1-C12 vector for DMPC, C2-C16 vector for CER[AP] and the sphingosine chain of CER16, and C2-C14 vector for the N-linked chain of CER16) and the bilayer normal. The phosphorusnitrogen vector was considered for P-N tilt analysis. The average tilt angle was calculated from the first moment of the distribution. The overall SA/lip was obtained by the area of the simulation box, which is the square of the edge length for a tetragonal cell, divided by the number of lipids per leaflet and is a reasonable model for the small systems presented in this work. Component SA/lip values were calculated by determining the X and Y coordinates of representative atoms for each lipid (C2, C21, and C31 for PC lipids, and C2S, C1F, and C4S for CER lipids as shown in Figure 1). Quickhull39 uses the coordinates of these atoms to create a tessellation in which each atom is represented by a polygon with calculable area. The sum of the polygonal area is the area of the simulation box, and summing areas of representative atoms gives the SA/lip for individual lipid components. The area compressibility modulus (KA) was calculated by the formula, "# $〈!〉 ! = %& ( 〈'〉

(3)

where kB is Boltzmann’s constant, T is the absolute temperature, 〈)〉 is the average SA/lip, N is the number of lipids per leaflet, and *〈!〉 is the variance of the area.

EDPs were determined by recentering the bilayer to place the interface at Z=0 Å, calculating densities for all atoms with a slab thickness of 0.2 Å, and combining atom densities 7 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

Page 8 of 38

into group and total densities. The overall bilayer thickness (DB), headgroup-to-headgroup distance (DHH), and hydrophobic distance (2DC) were calculated from EDPs. DB is defined as the distance between the half-maximums of the water EDP, DHH is the distance between the peaks of the total EDP, and 2DC is the distance between half-maximums of the hydrophobic EDP. For comparison with experimental neutron scattering length density (NSLD) data, SIMtoEXP40 was used to obtain an NSLD profile. From the profile, the distances to the primary and secondary peaks were averaged and the distance between the averages was calculated as the membrane thickness (dm). The thickness of the hydrophobic core (D) was calculated as the distance between half-maximums. Two-dimensional (2D) RDFs were calculated with a ∆r value of 0.1 Å for CER-CER, PC-PC, and PC-CER interactions. The phosphorous atom (P) of PC was considered for PC-PC RDFs and the carbon (C2) atom was considered for PC-CER RDFs. For both CER-CER and PCCER RDFs, the C2S atom of CER is used. Additionally, RDFs for the amide-carbonyl H bond used the N atom of CER and O22 atom of PC, and RDFs for the hydroxyl-carbonyl H bond used the O4 atom of CER[AP] and O22 atom of PC. The number of intramolecular and intermolecular H bonds (NHB) was calculated using the definition of an H bond as the distance between the donor and acceptor pairs being less than 2.4 Å and the donor-acceptor angle being greater than 150°. The results report the fraction of lipids with the considered H bond. Lipid clustering was analyzed using a python scikit-learn software41 with a density-based spatial clustering of applications with noise (DBSCAN) algorithm42. A cut-off distance of 5.0 Å between head groups was used for PC-CER clustering and a distance of 6.0 Å was used for PC-only clustering. Only cluster sizes of three lipids and more were considered. All reported p-values are based on t-test for equal means. 3.

Results

3.1. Overall Surface Area per Lipid (SA/lip) Overall SA/lip is a simple measure of lipid size in the membrane, describing the average lipid density without considering the natural differences in area among lipid types. Plots of the SA/lip as a function of time (Figure S1) generally show quick equilibration within 100 ns except for DMPC/CER[AP] at XCER=0.14 which relaxes much slower than other systems. Therefore, this system required much longer equilibration time and was ran until 1500 ns. The average SA/lip (Table 1) is largest for systems involving POPC and low CER concentrations with the 8 ACS Paragon Plus Environment

Page 9 of 38

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


largest value being 65.2 Å2 in POPC/CER16 at XCER=0.10. This effect of SA/lipid on CER concentration is expected since the absence of a large head group encourages tight packing in CER. This eclipses other contributions from H bonding as the addition of H-bond donors and loss of a double bond from CER16 to CER[AP] does not significantly affect SA/lip (p = 0.25). Likewise, the decrease of hydration in the DMPC/CER16 system does not produce any significant change in SA/lip. In contrast, the SA/lip of POPC/CER16 at XCER=0.60 is 25% greater than DMPC/CER16 at the same CER concentration (p=0.00003), which suggests that the unsaturated tail of POPC strongly opposes the tight packing encouraged by CER. Even between gel phase systems, these trends strongly correlate with ordering of lipid tails which will be explained in the next section. In this case, the gel phase is defined as a system having an overall SA/lip below 50 Å2 which agrees with tight packing that is characteristic of the gel phase. 3.2. Deuterium Order Parameters (SCD) Previous studies have compared the C36 force field43 to experimental data44-47 using pure bilayers of DMPC and POPC through comparison of order parameters and X-ray/neutron form factors, but comparison of PC and CER has so far only been provided for DMPC/CER16 systems through the GROMOS force field5, 18. The order parameters of this study are compared against both the GROMOS force field and NMR measurements as well as the NMR measurements of POPC/CER16 systems19, but the lack of site-specific deuteration in these studies often necessitates sorting values to be monotonically decreasing.

Figure 2. Comparison of SCD against experimental and GROMOS data for POPC lipids at XCER=0.10. Crosses represent the N-linked chain of CER16 while squares represent the sn-1



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 10 of 38

chain of POPC. Note that GROMOS data was run at 303K while C36 and NMR data were taken at 330K. Generally, SCD of the C36 force field are in good agreement with previous experimental and computational studies with some discrepancies. In DMPC/CER16 systems, the results from C36 simulations consistently predict lower order than NMR data (Figure S2), and since the NMR was carried out at a higher temperature of 325K, the true discrepancy is expected to be slightly larger. The gel-transition temperature of DMPC/CER16 mixtures is not precisely characterized because of CER’s tendency to broaden the transition temperature range48, but the tight packing of CER drastically increases the transition temperature of pure CER which is found to be around 353-363K49. CER therefore has been found to increase the melting temperature in PC/CER mixtures50 in a chain length dependent manner51. For DMPC/CER16 membranes at XCER=0.20, the transition occurs over the range 305-323K52 which overlaps the simulation of this work at 310K. Since the simulation occurs within the transition range, it is possible that the C36 force field inaccurately characterizes the phase transition which causes the observed decrease in order. This is supported by the strong agreement of POPC/CER16 membrane systems run at 330K, well above the melting temperature of 308K47, with experiment (Figure 2 and S3). C36 and GROMOS have differences in DMPC and POPC order parameters and large differences in CER16 order parameters as the C36 predicts similar order between DMPC and CER16 lipids while GROMOS predicts ~20% increase in CER16 (Figure S4). This is consistent with previous results demonstrating order parameter drop-off occurring at elevated temperatures for the GROMOS Berger force field17. The same study noted that the CHARMM force field better captured phase transitions, and the results of this work suggest that C36 better predicts the order of CER16 while less accurately predicting the order of DMPC. Direct comparison of POPC/CER16 might suggest that C36 is more accurate, but the higher order in GROMOS could also be explained as a result of the lower temperature. A tail comparison between all systems is given in Figure S5. The sn-2 tail of PC shows a strong dependency on CER concentration. For DMPC, the XCER=0.60 systems all display high order that plateau at the interior carbons which is characteristic of the gel phase. However, even systems in the gel phase show differences depending on XCER as DMPC/CER[AP] increases 20% from XCER=0.14 to XCER=0.60. The effect of lowering hydration is not statistically significant. In contrast, the unsaturation in POPC has a dramatic effect as the significant drop at the C9-C10


Page 11 of 38

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


carbons demonstrates. Even at carbons far away from the double bond, the order is lower than in saturated tails of DMPC. The sn-1 tail of PC tells the same story, but the N-linked and sphingosine tails of CER show the effect of adding hydroxyls. The N-linked chain of CER[AP] has much lower order at the C2 carbon where a hydroxyl is added whereas the sphingosine chain increases at the C4-C5 carbons which indicates the removal of the double bond is the dominating effect. The order parameters of these binary mixtures agrees with trends expected from pure PC membranes43 as order increases with XCER. Pure CER bilayers demonstrate similar or lower order than high XCER PC/CER mixtures which is likely a result of residual disordered lipids in pure CER systems. 3.3. Component Surface Areas (SA) and Area Compressibility Modulus (KA) The inability of overall SA/lip to differentiate between lipid types is resolved by analysis of component SA/lip (Table 2). For all systems, increasing XCER to 0.60 causes a clear decrease in PC SA/lip although the extent varies depending on the CER type. A 20% decrease is apparent in DMPC/CER16 systems whereas only a 4% decrease appears for DMPC/CER[AP] systems. The initial XCER varies from 0.20 for DMPC/CER16 to 0.14 for DMPC/CER[AP], so one would expect the discrepancy at identical concentrations to be even larger. An examination of the chain ordering reveals a larger difference in the PC ordering for DMPC/CER16 systems than DMPC/CER[AP] systems (60% vs 20%) due to the entrance into the gel phase (Figure S5). POPC/CER16, which never enters the gel phase and only moderately increases in chain order, shows a smaller decrease of 10% for a larger concentration difference which further indicates that phase transition drastically impacts the SA/lip through the chain ordering. Table 2. SA/lip (Å2) of specific lipid types in bilayer mixtures. Subscript LH refers to the lower hydration system. System

XCER

PC (Å2)

CER (Å2)

DMPC/CER16 DMPC/CER16LH DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16

0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60

56.6±0.2 56.2±0.1 44.9±0.2 -47.2±0.4 45.1±0.4 -66.1±0.1 58.5±0.2

50.2±0.1 49.8±0.5 42.1±0.2 43.6±0.3 40.5±0.3 43.2±0.5 45.3±0.6 57.1±0.2 51.3±0.1 11

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

Page 12 of 38

The SA/lip of CER generally follows similar patterns although to a lesser extent as the SA/lip generally decrease 10-15% as XCER increases to 0.60. However, an exception occurs in the DMPC/CER[AP] system in which SA/lip instead increases ~7%. The chain order of CER lipids increases in these systems so an appeal to ordering is unsatisfactory. Previous results found a similar increase for mixtures of SM and cholesterol (Chol) which were explained by Chol’s disruption of H-bond networks between SM lipids which allows the lipids to separate from each other while maintaining chain order53. CER is structurally similar to SM and induces similar membrane properties as Chol as a “Chol surrogate”18, so a similar appeal to H-bonding effects is reasonable here. For CER SA/lip, the most significant H bond is the CER-CER amide-carbonyl bond which universally increases with CER concentration (see section 3.6). Since this H bond involves multiple CER lipids, it is expected that this would have a larger impact on CER SA/lip than CER-PC bonds. The significance is especially notable in the POPC/CER16 system which demonstrates a decrease in CER SA/lip despite a decrease in overall inter H bonding. For CER16 in POPC, the 5-fold increase in CER-CER amide-carbonyl bonding outweighs the 0.05), so the hydroxylation of CER and double bond of POPC have no effect on the P-N oritentation and the main effect is the concentration of CER.


Page 17 of 38

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 4. Angle distributions of the P-N vector relative to the bilayer normal for a) DMPC/CER16, b) DMPC/CER[AP], and c) POPC/CER16 at varying CER concentrations. Pure DMPC was calculated from previous trajectories62 and is provided for comparison. 3.5. Thicknesses and Electron Density Profiles (EDPs) Further comparison with experiment8 can be achieved through the comparison of bilayer thicknesses calculated from NSLD profiles of the DMPC/CER[AP] system at XCER=0.14 (Figure 5). Profiles calculated from C36 data show the appearance of dual secondary peaks at around 1.8 nm and 2.2 nm. These peaks correspond to the polar moieties of PC and CER with the different distances of the phosphate and hydroxyl groups producing the dual peak pattern. These dual peaks are not found in experiment which may result from experimental Fourier transforms providing excess smoothness and consolidating both peaks into one. The secondary peaks reduce the suitability of Gaussian fitting to calculate the membrane thicknesses as was carried out experimentally, so the thickness dm for this work is calculated as the distance between peak 17 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

Page 18 of 38

averages and D as the distance between half-maximums, both of which are described in the Methods section. A comparison of the thicknesses yields reasonable agreement although C36 predicts a slightly thinner membrane for both overall and hydrophobic thickness. This likely results from the experiment being carried out at 60% humidity which can produce a slightly thicker membrane than at full hydration63 as was done in this work.

Figure 5. Comparison of membrane thickness (dm) and hydrophobic thickness (D) between experimental NSLD and C36 NSLD data for DMPC/CER[AP] at XCER=0.14. Thicknesses shown here are not calculated from EDP and are not the same as DB and 2DC. b) NSLD profile from C36 simulations used to calculated dm and D in (a) The thicknesses dm and D were only calculated for the DMPC/CER[AP] system to serve as comparison with experiment. The bilayer thicknesses DHH, DB, and 2DC are calculated from EDPs and serve as the main source of analysis (Table 4). Lowering hydration, DMPC/CER16LH, has minimal effect on all bilayer thicknesses, implying that this level of hydration (24 waters per lipids) is at or above the full hydration limit for this mixture. Slight increases are seen as might 18 ACS Paragon Plus Environment

Page 19 of 38

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


be expected from lower hydration, but these differences are statistically insignificant (p > 0.1). A clear increase in all membrane thicknesses is apparent for CER16 systems as XCER increases to 0.60, which is due to the greater ordering of the acyl chains. The hydrophobic thickness 2DC increases most significantly as it is the most impacted by chain ordering. This is especially pronounced for DMPC/CER16 in which 2DC increases by 20% whereas DHH increases 8% and DB increases 15%. Similarly, POPC/CER16 sees 2DC increase the most at 11% while DHH increases 9% and DB 6%. One might expect tilting of the gel phase to diminish membrane thickening, but the large concentrations of CER needed to enter the gel phase also produce relatively upright membranes so that chain ordering is dominant. CER also introduces an effect which opposes the increase of chain ordering which is the greater concentration of polar groups (hydroxyls) deeper within the membrane. In DMPC/CER[AP] systems, increasing the CER[AP] concentration affects DHH which decreases instead of increasing because of greater chain ordering as DB and 2DC do. Therefore, for high XCER, using DHH as a measure of bilayer thickness is not always consistent. Indeed, an increase in DHH is seen in for CER[AP] as XCER increases from 0.60 to 1.00 which is inconsistent with the expected decrease due to the lack of phosphocholine groups as seen in pure CER16. Instead of using DB and DHH, the depth of water penetration can be evaluated as the difference between DB and 2DC where smaller differences reflect greater penetration into the bilayer. Using this metric, it is clear that increasing CER concentration universally increases water penetration due to polar head groups being located deeper within the membrane. This might be counterintuitive given CER’s importance in barrier function of the SC, but the SC does not contain any phosphate moieties which would draw water away from CER head groups. For bilayers lacking phosphates, increasing CER concentration would instead decrease water permeability as the tight packing, high order, and increased hydrophobic thickness become dominant.



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 20 of 38

Table 4. Thicknesses (Å) of each system for different regions. DHH is the head-to-head distance, DB is the overall bilayer thickness, and 2DC is the hydrophobic thickness. Subscript LH refers to the lower hydration system. System

XCER

DHH

DB

2DC

DMPC/CER16 DMPC/CER16LH DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16

0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60

39.5±0.1 39.8±0.1 42.8±1.8 38.5±0.2 44.1±0.5 38.5±0.3 40.6±0.3 39.3±0.1 42.9±0.1

38.3±0.1 38.5±0.1 44.4±0.3 44.0±0.3 44.7±0.5 45.9±0.6 45.2±0.6 37.7±0.1 40.2±0.1

28.2±0.1 28.4±0.1 34.7±0.3 35.5±0.3 31.6±0.4 34.7±0.4 35.8±0.5 28.8±0.1 32.2±0.1

A strong dependency is shown with CER concentration both in the total and component EDPs. At higher CER concentrations, the total EDP becomes steeper in the center of the membrane and flatter in the 5-10 Å region which correlates with an increase in chain ordering diminishing orthogonal movement (Figure 6). The most notable effect is the appearance of dual peaks in gel phase systems which correspond to the different head group distances for PC and CER. These peaks are very distinct in the gel phase whereas fluid phases such as POPC/CER16 at XCER=0.60 instead demonstrate a flattening out of the head group region. Both the flattened and dual peak curves represent a blurring of the bilayer edge which can complicate determination of bilayer thickness as mentioned previously. Fluid PC/CER16 systems at low XCER show a shoulder in the head group region which transitions either into a dual peak or flattened curve depending on whether a phase transition occurs. DMPC/CER[AP] possesses dual peaks even at low XCER, but the phosphate peak is notably larger than the CER peak, but the amplitudes even out as XCER increases to 0.60. Pure CER systems generate smaller head group peaks due to the lack of phosphates which are shifted outward from CER head group peaks of PC/CER mixtures due to bilayer thickening, especially in the hydrophobic region.


Page 21 of 38

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 6. Total electron density profiles at various CER concentrations of a) DMPC/CER16 mixtures. The lower hydration systems are very similar to the full hydration systems and are omitted for clarity. b) DMPC/CER[AP] systems and c) POPC/CER16 systems Examination of component EDPs can reveal details not apparent from total EDPs (Figure S8-9). Perhaps the most notable effect is the drastic increase of the sphingosine EDP which is expected from increasing CER concentration. The broadness of the peak, calculated as the halfmaximum width, does not change despite increasing the concentration. For DMPC/CER16, this width is 6.2±0.1 Å for XCER=0.20, 5.9±0.1 Å for XCER=0.60, and 6.3±0.3 Å for pure CER. With the increase in sphingosine/phytosphingosine comes a decrease in PC head groups such as phosphate and choline. An increase in hydroxyl EDP is shown, especially for CER[AP] systems which have more hydroxyl groups present. All head groups move outward with XCER which is in agreement with general increases in bilayer thicknesses. In contrast, the water EDP shifts inward 21 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

Page 22 of 38

with a bend appearing at around 75% of the bulk density. This represents an increase of water penetration into the bilayer due to the higher sphingosine/phytosphingosine density and is reflected in the bilayer thicknesses. As pure CER does not see the same bend, this is likely due to the differing distances of polar groups in PC and CER. Another effect of CER is the increase of methyl density and decrease of methylene density in the bilayer center, but this effect diminishes at high concentration as pure CER shows little difference compared to XCER=0.60 systems. The decrease of methylene density is most significant and is the dominating factor that causes the dip in the total EDP. Changes in both densities represent a retardation of orthogonal movement which is highly correlated with the chain ordering and occurs in both the fluid and gel phases. Lastly, there are some notable differences in the EDP between POPC and DMPC, mainly resulting from unsaturation in POPC tails producing methine peaks. These peaks cause methylene peaks in POPC to become less smooth than in DMPC, but this effect reduces at high XCER because of lower methine density. 3.6. Hydrogen Bonding (H bonding) Pure PC bilayers do not have H bonding potential due to the lack of a bond donor in the PC headgroup, but sphingolipids such as CER have high bonding potential due to the presence of both donor and acceptor groups. CER[AP] is especially significant for having two additional hydroxyl groups that increase H bonding. Overall H bonding and individual pair probabilities for the most bonded pairs are reported (Table 5). Based on overall probabilities, inter H bonding is important for both PC and CER with both lipids generally having greater than 10% of lipids being involved in an inter H bond. For systems involving CER[AP], the probability can be greater than one suggesting that many lipids are involved in more than one H bond, which is not surprising given the large number of donor groups on CER[AP] and the presence of multiple acceptor groups in PC. CER often has greater H-bonding probability than PC, but this difference diminishes at high CER concentration where PC greatly increases H bonding while CER either marginally increases or decreases. The most stark contrast is in the DMPC/CER[AP] systems which transition from a CER:PC H-bonding ratio of 6.7 to 1.0. At identical XCER, DMPC tends to have more H bonding than POPC which is related to both the transition into the gel phase and the lower temperature that stabilize H-bond formation. Pure CER[AP] possesses higher H bonding compared to pure CER16 which agrees with both experimental64 and computational17 results. This increase is due to significant amounts of hydroxyl-hydroxyl bonding involving the 22 ACS Paragon Plus Environment

Page 23 of 38

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


O4 and O2 donors which are unique to CER[AP]. There is not one particular hydroxyl-hydroxyl bond that dominates, but the total contribution from several different pairs creates a strong effect (Table S3). The hydroxyl-hydroxyl bonds weaken in binary mixtures due to the lower concentration of CER and a preference for CER to bond with PC over other CER lipids. In contrast to inter H bonding, intra H bonding is negligible with very low probabilities that have large standard errors (Table 5). This is in contrast to SM bilayers which possess extremely high intra-H-bonding capabilities21. The absence of a phosphate group which is flexible enough to bend toward donor groups contributes to nearly all of the intra H bonding in SM. Table 5. Intermolecular and intramolecular hydrogen bond (defined as when the distance between donor and acceptor is less than 2.4 Å and the angle is greater than 150º) count for overall lipids in each system. Values are reported as the fraction of lipids with the considered hydrogen bond. Standard errors are represented in parentheses. Subscript LH refers to the lower hydration system. Type Inter

Lipid PC

CER

Intra

CER

System DMPC/CER16 DMPC/CER16LH DMPC/CER16 DMPC/CER[AP] DMPC/CER[AP] POPC/CER16 POPC/CER16 DMPC/CER16 DMPC/CER16LH DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16 DMPC/CER16 DMPC/CER16LH DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16

XCER 0.20 0.20 0.60 0.14 0.60 0.10 0.60 0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60 0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60

NHB 9.02±0.24 9.27±0.61 21.72±0.60 21.54±1.47 41.44±0.91 3.97±0.05 11.34±0.08 10.68±0.08 10.94±0.21 34.62±0.60 34.76±0.76 23.08±0.73 63.79±1.42 58.75±2.17 4.42±0.06 24.36±0.13 0.10±0.05 0.03±0.02 0.22±0.05 0.48±0.08 0.39±0.38 2.96±0.79 7.81±0.60 0.04±0.01 0.26±0.04

NHB/lip 0.113±0.003 0.116±0.008 0.543±0.015 0.250±0.017 1.036±0.023 0.044±0.001 0.283±0.002 0.534±0.004 0.547±0.011 0.577±0.010 0.348±0.008 1.648±0.052 1.063±0.024 0.587±0.022 0.442±0.006 0.406±0.002 0.005±0.003 0.002±0.001 0.004±0.001 0.005±0.001 0.028±0.027 0.049±0.013 0.078±0.006 0.004±0.001 0.004±0.001



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 24 of 38

Inter H bonding can be divided into CER-PC and CER-CER interactions and into the individual bond pairs, although hydroxyl-hydroxyl bonds will be excluded as these have been investigated in the previous paragraph. The most dominant bond pair is the CER-PC amidecarbonyl interaction (Table S4) which is expected since both pairs are located at similar distances along the bilayer normal. For CER, the probability of this bond ranges from 0.13-0.59 whereas the range is 0.03-0.39 for PC. The larger range for PC exists simply because low XCER systems naturally exclude many PC lipids from engaging in CER-PC H bonds. Accordingly, increasing XCER tends to decrease bonding for CER and increase bonding for PC. As PC cannot form H bonds with other PC lipids, this suggests that PC increases in PC-CER bonding which can be verified in Table S4. Greater association with CER diminishes crowding of PC choline groups which increases the P-N vector tilt as previously mentioned. Many CER-PC bonds are present in small amounts, generally having probabilities less than 0.10. Exceptions exist for systems involving CER[AP] which show bonding in the added hydroxyls at the 2nd and 4th carbon positions. The O4 hydroxyl-carbonyl bond is especially prominent with a probability of 0.362 at XCER=0.14. CER-CER inter H bonds show interesting trends as the amide-carbonyl bond increases with XCER but the hydroxyl-carbonyl bond does not change (Table S3). An anomaly exists for the DMPC/CER[AP] system at XCER=0.14 which shows no amide-carbonyl bonding. CER lipids tend to bond with PC lipids over other CER lipids even at high XCER, so these compete with CER bond acceptors and decrease CER-CER bonding. Since DMPC/CER[AP] exists in the gel phase, this reduces mobility so that rare and transient H bonds do not appear as they do for POPC/CER16 in the fluid phase. It is clear from our work that studies of CER systems need not consider intra H bonding for it is not significant (Table S5) and requires unfavorable conformational changes with CER’s small head group. 3.7. Lipid Clustering Clustering describes the local packing of lipids and is often driven by H bonding. The time evolution of lipid clustering could be used as a measure of convergence. Comparing with SA/lip over time (Figure S10), a negative correlation can be seen, but clustering does appear to reach a steady state after 100ns and is equilibrated within our analysis period. CER concentration, which has profound effects on H bonding as previously described, affects clustering of both PC and CER lipids (Table 6). Increasing CER concentration universally 24 ACS Paragon Plus Environment

Page 25 of 38

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


increases lipid clustering which correlates with a strong increase in inter H bonding (Figure 7). Increasing XCER to 0.60 generally increases the number of clusters 2- to 2.5-fold. Unlike other membrane properties, cluster number is not heavily impacted by phase transition as systems that enter the gel phase and systems that remain fluid increase clustering by similar amounts. Further increasing XCER to 1.00 causes similar increases in cluster number of around 2.5-fold from XCER=0.60. Clusters containing only CER lipids are considered with a 5.0 Å cutoff, and all systems show an average near 4 clusters with the difference between CER16 and CER[AP] being insignificant (p=0.17). As PC lipids have larger head groups than CER lipids, clustering analysis focusing exclusively on PC lipids with a larger cutoff distance of 6.0 Å is also considered. At low XCER, the clustering is greatest with a strong increase for DMPC/CER[AP] in the gel phase. At high XCER, this phase dependency disappears as all systems have between 3.0 and 3.5 clusters. Lipid ratios in clusters also is affected by CER concentration with a general increase in CER fraction and decrease in PC fraction (Table S6). This causes a tendency for the cluster fractions to approach the membrane fractions as PC is underrepresented in cluster relative to its concentration in the membrane at low CER concentration. This is likely a consequence of low H bonding in PC lipids which fundamentally lies in competition between PC lipids for CER’s Hbond donors at low CER concentration. This tendency for underrepresentation is not broken even at high CER concentrations and is fundamentally a consequence of PC’s lack of an H-bond donor group. Any H bonds involving PC will also involve CER, so CER will always be overrepresented in the cluster fraction. Furthermore, there is competition between PC lipids for H bonds at low XCER and an increase in CER-CER H bonding at high XCER, both of which drive fractional imbalance. Clustering can also explain the component SA/lip trends previously mentioned. As PC molecules tend to participate in clusters more at higher CER concentrations, the SA/lip will tend to decrease. CER is more complicated as a decreased tendency to participate in clusters can cause an increase in SA/lip as is the case in DMPC/CER[AP], but activity outside the clusters can dominate as chain ordering increases and encourages tight packing.



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 26 of 38

Table 6. Average number of total clusters (NT), CER-only clusters (NCER), and PC-only clusters (NPC) in bilayer mixtures. Cluster sizes of 3 lipids and greater were considered with a cutoff distance of 5.0 Å for NT, 5.0 Å for NCER, and 6.0 Å for NPC. NCER for low-XCER is not considered due to the rarity of clusters. Subscript LH refers to the lower hydration system. System DMPC/CER16 DMPC/CER16LH DMPC/CER16 CER16 DMPC/CER[AP] DMPC/CER[AP] CER[AP] POPC/CER16 POPC/CER16

XCER 0.20 0.20 0.60 1.00 0.14 0.60 1.00 0.10 0.60

NT 4.81±0.05 4.81±0.11 12.03±0.68 30.78±0.69 5.74±0.36 13.79±0.84 31.73±1.0 3.60±0.02 7.55±0.12

NCER --4.08±0.20 30.78±0.69 -3.63±0.18 31.73±1.0 -3.64±0.05

NPC 6.87±0.05 6.87±0.08 3.25±0.06 -15.13±0.83 3.42±0.08 -7.20±0.10 3.28±0.05


Page 27 of 38

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 7. Visual representation of lipid clustering for top leaflet in a) DMPC/CER[AP] at XCER=0.14, b) DMPC/CER[AP] at XCER=0.60, and c) pure CER[AP]. Red circles represent DMPC lipids while yellow circles represent CER[AP]. Clusters are taken from the final snapshot of the simulation rather than from a time average. 3.8. 2D Radial Distribution Functions (2D RDFs) 2D RDFs provide a more detailed measure of lateral ordering in the membrane that is specific to distances, although this is limited to atom pairs located approximately in the same plane. PC-PC interactions demonstrate strong dependency on CER concentration although there is a consistent shape made up of two major peaks with distances at 6.5 and 8.6 Å (Figure 8). These distances correspond to phosphate-phosphate conformations of PC lipids involved in 27 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

Page 28 of 38

clusters that are driven by H bonding with CER lipids (Figure 9). The amplitudes of these peaks are heavily dependent on CER concentration with the most notable case being DMPC/CER[AP] at XCER=0.14 demonstrating a large first peak that reaches a value of 1.6. Such a large peak is not found in the RDFs of pure DMPC65 which further suggests interactions with CER lipids promote local PC clustering. Clustering analysis supports that this first peak could be driven by DMPC clusters as DMPC cluster fractions increase at lower XCER which would increase the proportion of DMPC in close proximity which is represented by the first peak. As XCER increases, the first peak shrinks to become equal to the second peak which indicates no preference for either conformation, which is also true for CER16 systems. This is explained by decreasing cluster fractions so that PC lipids are less often in close proximity which eliminates preference for the close phosphate-phosphate conformation.

Figure 8. Comparison of PC-PC (phosphate-phosphate) (a-c) and CER-CER (C2S-C2S) (d-f) 2D RDFs at different CER concentrations for a,d) DMPC/CER16 systems, b,e) DMPC/CER[AP] systems, c,f) POPC/CER16 systems. Lighter shades of each color represent ±standard error. 28 ACS Paragon Plus Environment

Page 29 of 38

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


DMPC/CER16 at XCER=0.20 and lower hydration is nearly identical to the full hydration system and is omitted for clarity. CER-CER interactions indicate high clustering with a peak consistently located at 5.1 Å (Figure 8). Peak amplitude is higher than in PC-PC RDFs as peak amplitudes can exceed 2.0, but there is no distinct secondary peak. Increasing CER concentration has a tendency to decrease peak amplitude which is likely due to a decrease in CER-CER H bonding which disturbs CERCER clustering. Clustering analysis shows an increase in CER clustering which may seem contradictory, but that analysis does not distinguish between CER-CER and PC-CER interactions. POPC/CER16 at XCER=0.10 is noteworthy for consistently remaining below the bulk density for large distances up to 30 Å which may simply be a result of the low concentration. DMPC/CER[AP] at XCER=0.14 also does not quickly approach the bulk density and does not show the same increase in peak amplitude, but the standard interval is too large to make definite statements about the interactions occurring. The tendency for DMPC/CER[AP] systems to have larger standard errors and more fluctuations at far distances is related to a gel phase transition as lower mobility translates into many secondary peaks in the RDF. DMPC/CER16 at XCER=0.60, which exists in the gel phase as well, also shows signs of this behavior although to a lesser degree. PC-CER RDFs show similar activity to CER-CER RDFs although with lower peak amplitudes ranging from 1.1-1.9 and similar peak distances at 4.9 Å (Figure S11). CER concentration generally shows the same trend of decreasing peak amplitude, especially for POPC/CER16 where the peak reduces to near the bulk density. However, an exception exists for DMPC/CER16 systems where the peak amplitude significantly increases (p=0.007). This is surprising given that both PC-PC and PC-CER RDFs decrease with XCER but is likely a result of increased CER inter H bonding (Table 5) which decreases in other systems and may be fundamentally related to a fluid to gel transition that other systems do not go through.



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 30 of 38

Figure 9. Snapshots of CER-PC complexes in DMPC/CER[AP] at XCER=0.14 showing a) phosphate-phosphate distances corresponding to the two major peaks in PC-PC RDFs and b) CER-PC CER amide-carbonyl bond and C2-C2S distance. Analysis of some of the most dominant H-bond pairs was determined from H-bond analysis as the CER-PC amide-carbonyl bond and the O4 hydroxyl-carbonyl bond (Figures S1213). For amide-carbonyl bonds, a very sharp first peak which corresponds to the first bond at is located at 2.8 Å. Although this is greater than the H-bonding cutoff, the RDF differs because it uses the heavy atom distance while H bonding uses the hydrogen-acceptor distance. Each atom selection is more common in its respective analysis, and this accounts for the greater RDF distance. Increasing CER concentration increases peak height which corresponds which greater bonding from H-bonding analysis. POPC systems notably have smaller peaks which agrees with a lower total bond count and is due to both higher temperature and remaining in the fluid phase which tends to disrupt bonding. DMPC/CER[AP] systems possess a secondary peak at 4.5 Å which possibly points to the existence of DMPC clusters that produce the major first peak seen in phosphate-phosphate RDFs. The O4 hydroxyl-carbonyl bond is unique to the DMPC/CER[AP] systems and is much weaker than amide-carbonyl bonding with peak amplitudes not exceeding 1.2. Secondary peaks at 4.6 Å occur similarly to amide-carbonyl RDFs and likewise support the existence of occasional DMPC clusters, but differ in that high CER concentration causes a higher peak for hydroxyl-carbonyl bonds instead of a lower peak as in amide-carbonyl bonds.


Page 31 of 38

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


4.

Discussion and Conclusion Computational studies involving binary mixtures of CER and PC lipids have been

relatively lacking in comparison to experimental research. Pandit et al. ran GROMOS FF simulations18, with parameters from previous work66, of POPC/CER16 bilayers which found increased chain order. However, these values are generally around 20% larger than NMR and C36 results from this work. Given the quick equilibration seen in our simulations of these systems, this is not likely a timescale issue and is more likely a combination of a lower temperature and inherent bias of the force field derived from Chiu et al67. Dutagaci et al.5 performed simulations using the Berger force field68 which provided qualitative understandings of how CER induces order in DMPC membranes. However, quantitative comparison suggests mixed results and reveals that GROMOS Berger tends to overestimate order compared to experimental NMR. Accordingly, decrease of SA/lip with the introduction of CER is too large as GROMOS Berger predicts a 10 Å2 decrease while experiment suggests a 4 Å2 decrease. The C36 force field exhibits this problem to a much milder extent, generally more accurately predicting order parameters and exhibiting an area decrease of 6.3 Å2 from 61.243 to 54.9 Å2. H-bond analysis results also do not agree as Dutagaci et al. found that hydroxyl-carbonyl pairs are more prevalent than amide-carbonyl pairs which contradicts the C36 results of this work. Experimental results from Ekman et al.69 support C36 results as the methylation of the amide attenuates ordered phase formation while methylation of the hydroxyl has no discernable effect. In general, although both C36 and Berger force fields correctly predict qualitative trends in chain ordering, C36 more accurately predicts ordering and H bonding in PC/CER bilayers, which may hint at better modeling of lipids in diverse membranes. There has been a wide variety of experimental studies on CER and PC/CER mixtures with a large majority finding evidence of phase separation with CER lipids generally occupying gel phases while PC lipids prefer liquid crystalline phases. Massey70 proposed the formation of DPPC-CER complexes (likely at the micrometer scale) which may be seen at the nano-scale for the DMPC/CER[AP] system of this work (DPPC is similar in structure to DMPC). (Figure 9). Microdomain formation has been suggested to be driven by multiple factors with the presence of both H-bond donors and acceptors in CER being vital as these allow the formation of H-bond networks that cannot be formed with PC. Although PC-CER interactions are more favorable than CER-CER interactions in PC/CER mixtures, more H bonds form in pure CER16 than in 31 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

Page 32 of 38

PC/CER16 mixtures, which suggests that CER16 domains are favorable but can be hindered over short timescales by strong PC-CER16 interactions. Chain saturation is also essential as unsaturated chains tend to produce more miscible CER lipids and lower ordered phase separation, and it has been proposed by Holopainen et al.71 that different chain structures in CER could lead to different functions depending on their effect on physiological membranes such as long chain CER promoting leaflet coupling versus short chain CER serving as second messengers in apoptosis. Rerek et al. make a similar claim in pointing out potentially different functions of sphingosine and phytosphingosine head groups in CER64. Phytosphingosine has greater H bonding which is supported by the 2-fold increase found in this work (Table 5) and contributes to proper barrier function in the SC. However, the double bond in sphingosine allows orthorhombic chain packing which also contributes to proper barrier function, so both head groups can exist in the SC and possess different roles. Sphingosine, with its propensity for orthorhombic packing, tends to form distinct domains as found by the previously mentioned studies, and phytosphingosine serves to join those domains through strong H bonding although phytosphingosine does not itself form distinct domains. The systems of this work are too small to see the distinct ordered domains at the µm-level seen in experimental studies and would need to be much larger as Pandit et al. similarly did not see any major lateral organization for systems with twice as many lipids. However, these systems possibly offer views of different areas in segregated PC/CER mixtures by varying the CER concentration since CER is known to populate gel phase domains in greater proportions. Additionally, systems involving CER[AP] offer further support for the hypothesis of Rerek et al. as both low and high CER[AP] concentrations induce the gel phase, demonstrating that CER[AP] does not favorably form segregated domains and is likely relegated to a role in joining such domains together. In this study, the nuance of H bonding in CER and PC/CER bilayers and its impact on various membrane properties is elucidated. H bonding is a strong driving force for transition into the gel phase in mixtures, and the balance of various H bonds can determine if CER expands in surface area despite greater chain ordering. Simulations performed by Guo et al.17 with a modified CHARMM force field found that hydroxyl H bonds are dominant over amide H bonds which contradicts the results of this work and experimental results from Ekman et al. as previously mentioned. Likely, this results from not considering the carbonyl group in H-bond analysis which is the most prominent acceptor group. Tilt angle in C36 ranges from 16-18º 32 ACS Paragon Plus Environment

Page 33 of 38

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


which is lower than the modified CHARMM values of 24.3 º and closer to the Berger GROMOS values of 17 º. The bilayer thickness of this work is also slightly greater as CER16 has a value of 44.0 Å2 whereas Guo et al. report a value of 42.5 Å2, but this is likely in part due to a difference in thickness definition as Guo et al. used 1/e of the maximum rather than half to calculate DB. However, the overall trend of greater thickness in phytosphingosine CER over sphingosine CER is maintained. The systems studied in this work are pure CER bilayers as well as binary mixtures of PC and CER with varying lipid type and CER concentration. Comparison of deuterium order parameters with experimental data showed good agreement with a slight underestimation for DMPC lipids which is attributed to an inaccurate prediction of the phase transition, blurred by the introduction of CER. However, good agreement with experiment is shown based on bilayer thickness calculated from NSLD data, and a general improvement in order and H-bond agreement is demonstrated over GROMOS FF simulations of similar studies. Consistency is shown in SA/lip and thicknesses calculated from previous C36 simulations of pure PC bilayers which supports the use of the C36 force field for ternary mixtures of PC and CER and potentially even more diverse membranes, although these simulations may be limited by the capability of current supercomputing architecture. CER is found to have considerable effects on membrane properties, often mediated through H bonding which is essential for the domain formation that is prevalent in rafts and in the SC. 5.

Supporting Information Available Figures for SA/lip as a function of time, order parameter plots, lateral clustering of lipids,

snapshots of pure ceramide bilayers, comparison of neutron scatting length density, component electron density profiles, radial distributions are included. Tables of average tilt angle, P-N vector orientation, lipid hydrogen bonds, and clustering analysis are included in the supporting materials. 6.

Acknowledgements This research is supported by the NSF grant MCB-1149187 and by a grant to the

University of Maryland from the Howard Hughes Medical Institute through the Science Education Program. The high performance computational resource used for this research is 33 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

Page 34 of 38

Deepthought which is maintained by the Division of Information Technology at the University of Maryland. Production runs used the Extreme Science and Engineering Discovery Environment (XSEDE) allocation by grant number MCB-100139 on Stampede and Bridges and the computational resources of the NIH HPC Biowulf cluster provided by Dr. Richard Pastor. We especially thank Xiaohong Zhuang for assisting with clustering analysis and providing pure PC data. 7.

References

1. van Meer, G.; Voelker, D. R.; Feigenson, G. W., Membrane Lipids: Where They Are and How They Behave. Nature Reviews Molecular Cell Biology 2008, 9 (2), 112-124. 2. Hla, T.; Dannenberg, A. J., Sphingolipid Signaling in Metabolic Disorders. Cell Metabolism 2012, 16 (4), 420-434. 3. Yasuda, T.; Kinoshita, M.; Murata, M.; Matsumori, N., Detailed Comparison of Deuterium Quadrupole Profiles between Sphingomyelin and Phosphatidylcholine Bilayers. Biophysical Journal 2014, 106 (3), 631-638. 4. Saghy, E.; Szoke, E.; Payrits, M.; Helyes, Z.; Borzsei, R.; Erostyak, J.; Janosi, T. Z.; Setalo, G.; Szolcsanyi, J., Evidence for the Role of Lipid Rafts and Sphingomyelin in Ca2+-Gating of Transient Receptor Potential Channels in Trigeminal Sensory Neurons and Peripheral Nerve Terminals. Pharmacological Research 2015, 100, 101-116. 5. Dutagaci, B.; Becker-Baldus, J.; Faraldo-Gomez, J. D.; Glaubitz, C., Ceramide-Lipid Interactions Studied by MD Simulations and Solid-State NMR. Biochimica Et Biophysica Acta-Biomembranes 2014, 1838 (10), 2511-2519. 6. Gulbins, E.; Kolesnick, R., Raft Ceramide in Molecular Medicine. Oncogene 2003, 22 (45), 7070-7077. 7. Gulbins, E.; Dreschers, S.; Wilker, B.; Grassme, H., Ceramide, Membrane Rafts and Infections. Journal of Molecular Medicine-Jmm 2004, 82 (6), 357-363. 8. Kiselev, M. A.; Zemlyanaya, E. V.; Ryabova, N. Y.; Hauss, T.; Almasy, L.; Funari, S. S.; Zbytovska, J.; Lombardo, D., Influence of Ceramide on the Internal Structure and Hydration of the Phospholipid Bilayer Studied by Neutron and X-ray Scattering. Applied Physics a-Materials Science & Processing 2014, 116 (1), 319-325. 9. Bouwstra, J. A.; Gooris, G. S.; Dubbelaar, F. E. R.; Weerheim, A. M.; Ijzerman, A. P.; Ponec, M., Role of Ceramide 1 in the Molecular Organization of the Stratum Corneum Lipids. Journal of Lipid Research 1998, 39 (1), 186-196. 10. White, S. H.; Mirejovsky, D.; King, G. I., Structure of Lamellar Lipid Domains and Corneocyte Envelopes of Murine Stratum-Corneum - An X-ray Diffraction Study. Biochemistry 1988, 27 (10), 37253732. 11. Bouwstra, J. A.; Ponec, M., The Skin Barrier in Healthy and Diseased State. Biochimica Et Biophysica Acta-Biomembranes 2006, 1758 (12), 2080-2095. 12. Khakbaz, P.; Klauda, J. B., Probing the Importance of Lipid Diversity in Cell Membranes via Molecular Simulation. Chem. Phys. Lipids 2015, 192, 12-22. 13. Lyubartsev, A. P.; Rabinovich, A. L., Force Field Development for Lipid Membrane Simulations. Biochimica Et Biophysica Acta-Biomembranes 2016, 1858 (10), 2483-2497. 14. Khakbaz, P.; Monje-Galvan, V.; Zhuang, X.; Klauda, J. B., Modeling Lipid Membranes. In Biogenesis of Fatty Acids, Lipids and Membranes, Geiger, O., Ed. Springer International Publishing: Cham, 2017; pp 1-19.


Page 35 of 38

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


15. Pandit, S. A.; Scott, H. L., Molecular-Dynamics Simulation of a Ceramide Bilayer. Journal of Chemical Physics 2006, 124 (1), 7. 16. Papadimitriou, N. I.; Kainourgiakis, M. E.; Karozis, S. N.; Charalambopoulou, G. C., Studying the Structure of Single-Component Ceramide Bilayers with Molecular Dynamics Simulations using Different Force Fields. Molecular Simulation 2015, 41 (13), 1122-1136. 17. Guo, S.; Moore, T. C.; Iacovella, C. R.; Strickland, L. A.; McCabe, C., Simulation Study of the Structure and Phase Behavior of Ceramide Bilayers and the Role of Lipid Headgroup Chemistry. Journal of Chemical Theory and Computation 2013, 9 (11), 5116-5126. 18. Pandit, S. A.; Chiu, S. W.; Jakobsson, E.; Grama, A.; Scott, H. L., Cholesterol Surrogates: A Comparison of Cholesterol and 16 : 0 Ceramide in POPC Bilayers. Biophysical Journal 2007, 92 (3), 920-927. 19. Hsueh, Y. W.; Giles, R.; Kitson, N.; Thewalt, J., The Effect of Ceramide on Phosphatidylcholine Membranes: A Deuterium NMR Study. Biophysical Journal 2002, 82 (6), 3089-3095. 20. Klauda, J. B.; Venable, R. M.; Freites, J. A.; O'Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W., Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. Journal of Physical Chemistry B 2010, 114 (23), 78307843. 21. Venable, R. M.; Sodt, A. J.; Rogaski, B.; Rui, H.; Hatcher, E.; MacKerell, A. D.; Pastor, R. W.; Klauda, J. B., CHARMM All-Atom Additive Force Field for Sphingomyelin: Elucidation of Hydrogen Bonding and of Positive Curvature. Biophysical Journal 2014, 107 (1), 134-145. 22. Jo, S.; Kim, T.; Iyer, V. G.; Im, W., CHARMM-GUI: a Web-based Graphical User Interface for CHARMM. J Comput Chem 2008, 29 (11), 1859-65. 23. Jo, S.; Kim, T.; Im, W., Automated Builder and Database of Protein/Membrane Complexes for Molecular Dynamics Simulations. Plos One 2007, 2 (9), 9. 24. Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; Davila-Contreras, E. M.; Qi, Y. F.; Lee, J. M.; Monje-Galvan, V.; Venable, R. M., et al., CHARMM-GUI Membrane Builder Toward Realistic Biological Membrane Simulations. Journal of Computational Chemistry 2014, 35 (27), 1997-2004. 25. Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W., CHARMM-GUI Membrane Builder for Mixed Bilayers and Its Application to Yeast Membranes. Biophysical Journal 2009, 97 (1), 50-58. 26. Durell, S. R.; Brooks, B. R.; Bennaim, A., Solvent-Induced Forces between Two Hydrophilic Groups. J. Phys. Chem. 1994, 98 (8), 2198-2202. 27. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926-935. 28. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K., Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781-1802. 29. Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R., Constant-Pressure Molecular-Dynamics Simulation - the Langevin Piston Method. Journal of Chemical Physics 1995, 103 (11), 4613-4621. 30. Martyna, G. J.; Tobias, D. J.; Klein, M. L., Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101 (5), 4177-4189. 31. Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald - an NLog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089-10092. 32. Steinbach, P. J.; Brooks, B. R., New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. J. Comput. Chem. 1994, 15 (7), 667-683. 33. Andersen, H. C., Rattle - a Velocity Version of the Shake Algorithm for Molecular-Dynamics Calculations. Journal of Computational Physics 1983, 52 (1), 24-34. 34. Humphrey, W.; Dalke, A.; Schulten, K., VMD: Visual Molecular Dynamics. Journal of Molecular Graphics & Modelling 1996, 14 (1), 33-38. 35. Racine, J., gnuplot 4.0: A Portable Interactive Plotting Utility. Journal of Applied Econometrics 2006, 21 (1), 133-141.



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 36 of 38

36. Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S., et al., CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30 (10), 1545-1614. 37. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., CHARMM - a Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4 (2), 187-217. 38. Dufourc, E. J.; Parish, E. J.; Chitrakorn, S.; Smith, I. C. P., Structural and Dynamical Detials of Cholesterol Lipid Interaction as Revealed by Deuterium NMR. Biochemistry 1984, 23 (25), 6062-6071. 39. Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H., The Quickhull Algorithm for Convex Hulls. Acm Transactions on Mathematical Software 1996, 22 (4), 469-483. 40. Kučerka, N.; Katsaras, J.; Nagle, J., Comparing Membrane Simulations to Scattering Experiments: Introducing the SIMtoEXP Software. J. Membr. Biol. 2010, 235 (1), 43-50. 41. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V., et al., Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 2011, 12, 2825-2830. 42. Martin Ester, H.-P. K., Jörg Sander, Xiaowei Xu, A Density-based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Proceedings of 2nd International Conference on Knowledge Discovery and Data Mining 1996, 96 (34), 226-231. 43. Zhuang, X. H.; Makover, J. R.; Im, W.; Klauda, J. B., A Systematic Molecular Dynamics Simulation Study of Temperature Dependent Bilayer Structural Properties. Biochimica Et Biophysica Acta-Biomembranes 2014, 1838 (10), 2520-2529. 44. Kucerka, N.; Nieh, M. P.; Katsaras, J., Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochimica Et Biophysica ActaBiomembranes 2011, 1808 (11), 2761-2771. 45. Douliez, J. P.; Leonard, A.; Dufourc, E. J., Restatement of Order Parameters in Biomembranes Calculation of C-C Bond Order Parameters from C-D Quadrupolar Splittings. Biophys. J. 1995, 68 (5), 1727-1739. 46. Seelig, J.; Waespesarcevic, N., Molecular Order in Cis and Trans Unsaturated Phospholipid Bilayers. Biochemistry 1978, 17 (16), 3310-3315. 47. Slotte, J. P.; Yasuda, T.; Engberg, O.; Al Sazzad, M. A.; Hautala, V.; Nyholm, T. K. M.; Murata, M., Bilayer Interactions among Unsaturated Phospholipids, Sterols, and Ceramide. Biophysical Journal 2017, 112 (8), 1673-1681. 48. Veiga, M. P.; Arrondo, J. L. R.; Goni, F. M.; Alonso, A., Ceramides in Phospholipid Membranes: Effects on Bilayer Stability and Transition to Nonlamellar Phases. Biophysical Journal 1999, 76 (1), 342350. 49. Shah, J.; Atienza, J. M.; Duclos, R. I.; Rawlings, A. V.; Dong, Z. X.; Shipley, G. G., Structural and Thermotropic Properties of Synthetic C16-0 (Palmitoyl) Ceramide - Effect of Hydration. Journal of Lipid Research 1995, 36 (9), 1936-1944. 50. Pinto, S. N.; Silva, L. C.; Futerman, A. H.; Prieto, M., Effect of Ceramide Structure on Membrane Biophysical Properties: The Role of Acyl Chain Length and Unsaturation. Biochimica Et Biophysica Acta-Biomembranes 2011, 1808 (11), 2753-2760. 51. Carrer, D. C.; Schreier, S.; Patrito, M.; Maggio, B., Effects of a Short-Chain Ceramide on Bilayer Domain Formation, Thickness, and Chain Mobililty: DMPC and Asymmetric Ceramide Mixtures. Biophysical Journal 2006, 90 (7), 2394-2403. 52. Holopainen, J. M.; Lemmich, J.; Richter, F.; Mouritsen, O. G.; Rapp, G.; Kinnunen, P. K. J., Dimyristoylphosphatidylcholine/C16 : 0-Ceramide Binary Liposomes Studied by Differential Scanning Calorimetry and Wide- and Small-Angle X-ray Scattering. Biophysical Journal 2000, 78 (5), 2459-2469. 53. Wang, E.; Klauda, J. B., Examination of Mixtures Containing Sphingomyelin and Cholesterol by Molecular Dynamics Simulations. J. Phys. Chem. B 2017, 121 (18), 4833-4844. 54. Nagle, J. F.; Tristram-Nagle, S., Structure of Lipid Bilayers. Biochimica et Biophysica ActaReviews on Biomembranes 2000, 1469 (3), 159-195. 36 ACS Paragon Plus Environment

Page 37 of 38

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


55. Janosi, L.; Gorfe, A. A., Simulating POPC and POPC/POPG Bilayers: Conserved Packing and Altered Surface Reactivity. Journal of Chemical Theory and Computation 2010, 6 (10), 3267-3273. 56. Tristram-Nagle, S.; Liu, Y. F.; Legleiter, J.; Nagle, J. F., Structure of Gel Phase DMPC Determined by X-Ray Diffraction. Biophys. J. 2002, 83 (6), 3324-3335. 57. Dahlen, B.; Pascher, I., Molecular Arrangements in Sphingolipids - Thermotropic PhaseBehavior of Tetracosanoylphytosphingosine. Chemistry and Physics of Lipids 1979, 24 (2), 119-133. 58. Paloncyova, M.; Vavrova, K.; Sovova, Z.; DeVane, R.; Otyepka, M.; Berka, K., Structural Changes in Ceramide Bilayers Rationalize Increased Permeation through Stratum Corneum Models with Shorter Acyl Tails. Journal of Physical Chemistry B 2015, 119 (30), 9811-9819. 59. Morrow, M. R.; Helle, A.; Perry, J.; Vattulainen, I.; Wiedmer, S. K.; Holopainen, J. M., Ceramide-1-Phosphate, in Contrast to Ceramide, Is Not Segregated into Lateral Lipid Domains in Phosphatidylcholine Bilayers. Biophysical Journal 2009, 96 (6), 2216-2226. 60. Huang, H. W.; Goldberg, E. M.; Zidovetzki, R., Ceramide Induces Structural Defects into Phosphatidylcholine Bilayers and Activates Phospholipase A(2). Biochemical and Biophysical Research Communications 1996, 220 (3), 834-838. 61. Tristramnagle, S.; Zhang, R.; Suter, R. M.; Worthington, C. R.; Sun, W. J.; Nagle, J. F., Measurement of Chain Tilt Angle in Fully Hydrated Bilayers of Gel Phase Lecithins. Biophysical Journal 1993, 64 (4), 1097-1109. 62. Zhuang, X.; Davila-Contreras, E. M.; Beaven, A. H.; Im, W.; Klauda, J. B., An Extensive Simulation Study of Lipid Bilayer Properties with Different Head Groups, Acyl Chain Lengths, and Chain Saturations. Biochimica Et Biophysica Acta-Biomembranes 2016, 1858 (12), 3093-3104. 63. Faure, C.; Dufourc, E. J., The Thickness of Cholesterol Sulfate-Containing Membranes Depends Upon Hydration. Biochimica Et Biophysica Acta-Biomembranes 1997, 1330 (2), 248-253. 64. Rerek, M. E.; Chen, H. C.; Markovic, B.; Van Wyck, D.; Garidel, P.; Mendelsohn, R.; Moore, D. J., Phytosphingosine and Sphingosine Ceramide Headgroup Hydrogen Bonding: Structural Insights through Thermotropic Hydrogen/Deuterium Exchange. Journal of Physical Chemistry B 2001, 105 (38), 9355-9362. 65. Hub, J. S.; Salditt, T.; Rheinstader, M. C.; de Groot, B. L., Short-Range Order and Collective Dynamics of DMPC Bilayers: A Comparison Between Molecular Dynamics Simulations, X-ray, and Neutron Scattering Experiments. Biophysical Journal 2007, 93 (9), 3156-3168. 66. Pandit, S. A.; Jakobsson, E.; Scott, H. L., Simulation of the Early Stages of Nano-Domain Formation in Mixed Bilayers of Sphingomyelin, Cholesterol, and Dioleylphosphatidylcholine. Biophysical Journal 2004, 87 (5), 3312-3322. 67. Chiu, S. W.; Jakobsson, E.; Subramaniam, S.; Scott, H. L., Application of a Combined Monte Carlo and Molecular Dynamics Method to the Simulation of a Dipalmitoyl PhosphatidylcholineCholesterol Mixed Lipid Bilayers. Biophysical Journal 1999, 76 (1), A58-A58. 68. Berger, O.; Edholm, O.; Jahnig, F., Molecular Dynamics Simulations of a Fluid Bilayer of Dipalmitoylphosphatidylcholine at Full Hydration, Constant Pressure, and Constant Temperature. Biophys. J. 1997, 72 (5), 2002-2013. 69. Ekman, P.; Maula, T.; Yamaguchi, S.; Yamamoto, T.; Nyholm, T. K. M.; Katsumura, S.; Slotte, J. P., Formation of an Ordered Phase by Ceramides and Diacylglycerols in a Fluid Phosphatidylcholine Bilayer - Correlation with Structure and Hydrogen Bonding Capacity. Biochimica Et Biophysica ActaBiomembranes 2015, 1848 (10), 2111-2117. 70. Massey, J. B., Interaction of Ceramides with Phosphatidylcholine, Sphingomyelin and Sphingomyelin/Cholesterol Bilayers. Biochimica Et Biophysica Acta-Biomembranes 2001, 1510 (1-2), 167-184. 71. Holopainen, J. M.; Brockman, H. L.; Brown, R. E.; Kinnunen, P. K. J., Interfacial Interactions of Ceramide with Dimyristoylphosphatidylcholine: Impact of the N-Acyl Chain. Biophysical Journal 2001, 80 (2), 765-775.



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 38 of 38

TOC Graphic


Molecular dynamics simulations of histidine-containing cod antimicrobial peptide paralogs in self-assembled bilayers.

Universal Method for Embedding Proteins into Complex Lipid Bilayers for Molecular Dynamics Simulations.

The Integrin Receptor in Biologically Relevant Bilayers: Insights from Molecular Dynamics Simulations.

Molecular dynamics simulations in biology.

Floating orbital molecular dynamics simulations.

Molecular dynamics simulations: advances and applications.

Microsecond molecular dynamics simulations of lipid mixing.

Molecular dynamics simulations of classical stopping power.

Molecular dynamics simulations of large macromolecular complexes.

Accelerated molecular dynamics simulations of protein folding.

Molecular dynamics simulations of helix denaturation.

Molecular dynamics simulations of glycoproteins using CHARMM.

Interactions of the EGFR juxtamembrane domain with PIP2-containing lipid bilayers: Insights from multiscale molecular dynamics simulations.

AceCloud: Molecular Dynamics Simulations in the Cloud.

New ways to boost molecular dynamics simulations.

Mechanical properties of lipid bilayers from molecular dynamics simulation.

Molecular dynamics simulations: insight into molecular phenomena at interfaces.

Quenched molecular dynamics simulations of tuftsin and proposed cyclic analogues.

Molecular dynamics simulations of dinucleoside and dinucleoside-drug crystal hydrates.

Force spectroscopy of polymer desorption: theory and molecular dynamics simulations.

Molecular Dynamics and Metadynamics Simulations of the Cellulase Cel48F.

Nonadiabatic molecular dynamics simulations: synergies between theory and experiments.

Nanofluidic osmotic diodes: theory and molecular dynamics simulations.

Lipid Models for United-Atom Molecular Dynamics Simulations of Proteins.