Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple

Feb 21, 2018 - The barrier function of the stratum corneum (SC) is intimately related to the structure of the lipid matrix, which is composed of ceram...
0 downloads 6 Views 1MB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

Article

Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple Interior Stratum Corneum Models Eric Wang, and Jeffery B. Klauda J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b00348 • Publication Date (Web): 21 Feb 2018 Downloaded from http://pubs.acs.org on February 23, 2018

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

Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple Interior Stratum Corneum Models 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

The barrier function of the stratum corneum (SC) is intimately related to the structure of the lipid matrix, which is composed of ceramides (Cer), cholesterol (Chol), and free fatty acid (FFA). In this study, the all-atom CHARMM36 (C36) force field is used to simulate bilayers of N-palmitoylsphingosine (Cer16), N-lignoceroylsphingosine (Cer24), Chol, and lignoceric acid (LA) as simple models of the SC. Equimolar mixtures of Cer, Chol, and LA are replicated from experiment for comparison and validation of the C36 force field, and the effect of lipid diversity and temperature is studied. The presence of Chol and LA have effects on nearly all membrane properties including surface area per lipid, area compressibility moduli, chain order, Chol tilt, bilayer thickness, interdigitation, hydrogen bonding, and lipid clustering while temperature has a more moderate effect. In systems containing Cer16, there is a profound difference in interdigitation between pure Cer and mixed systems while systems containing Cer24 are relatively unaffected. Increasing temperature has the potential to shift hydrogen bonding pairs rather than uniformly decrease bonding, which can lead to greater Cer-Cer bonding at higher temperatures. Comparison with deuterium order parameter experiments demonstrates good agreement which supports further use of these class of lipids and fatty acids for development of more complex SC models.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

Introduction The outer layer of the skin, known as the stratum corneum (SC), is mainly composed of

corneocytes, dead cells containing keratin and water, within a lipid matrix. The corneocyte envelope is an effective barrier against extracellular molecules which leads the majority of transdermal penetration to occur through the lipid matrix1. Barrier function is closely related to the concentrations and types of ceramides (Cer), a class of sphingolipids notable for possessing small headgroups lacking any phosphocholine moieties. For example, lowered Cer concentration has been suggested to induce irregularities such as dermatitis2, 3. Other major components of the SC include free fatty acids (FFA) and cholesterol (Chol), but phospholipids, common in physiological fluid membranes, are notably absent. Lipid organization within the SC has been studied through a variety of methods such as electron microscopy4, Fourier transform infrared spectroscopy5, and most significantly X-ray diffraction1,

6, 7

. From these methods, the existence of two distinct lamellar phases have been

determined. The short periodicity phase (SPP) has a repeat distance of approximately 6.4 nm while the long periodicity phase (LPP) has a repeat distance of approximately 13.4 nm. While the molecular structure of these phases is not concretely known, it is believed that the SPP is characterized as a bilayer7. The LPP is more complicated, and experimental models have proposed a multilayer structure with bilayers surrounding a monolayer center5. A variety of Cer types exist within the SC and differ based on their hydroxylation and fatty acid chain length. The most common base of Cer is sphingosine, and fatty acid chains are usually between 20 and 30 carbons long. For accurate structure determination of the SC, experimental methods such as X-ray diffraction are essential. However, these methods are resource intensive, and structural determination on the atomic level is difficult. To probe properties on the atomic level, molecular dynamics (MD) simulations are commonly used. Over the last few decades, the rapid expansion of computational power has increased the accessibility of MD for the simulation of a wide variety of biochemical systems. In the existing literature, there have already been investigations of the SC’s SPP using pure Cer, Chol/FFA, and Cer/Chol/FFA bilayers with the GROMOS force field8-11. These have found significant effects due to lipid molar ratios and chain length on bilayer properties. However, while numerous simulation studies with custom force field parameters12 and the GROMOS force field8-10 already exist, simulations using the CHARMM36 2 ACS Paragon Plus Environment

Page 3 of 36 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

(C36) force field, which is known to better agree with experiment for sphingolipids13, are not as prevalent. C36 force field has also been used to simulate pure Cer and Cer mixed with phosphatidylcholine (PC) and found better agreement with experiment over GROMOS in properties such as deuterium order parameters and hydrogen bonding14. Therefore, using the C36 force field to simulate more complex systems such as Cer/Chol/FFA is highly relevant for studying the SC. Recently, Moore et al.15 published a random walk molecular dynamics (RWMD) study of the SC’s SPP using the C36 force field which allowed for overcoming kinetic trapping due to biased initial builds. Sphingolipid parameters were not obtained from the parameters in Venable et al.13 but from previous work16, and LA was protonated. Nearly all past simulations of the SC utilize a protonated FFA. However, there is certainly some proportion of deprotonated FFA in the physiological SC, even in the outer acid mantle with pH ranging from 4.5-5.317. The pKa of FFA is around 518, suggesting a relatively evenly balanced ratio of protonated to deprotonated FFA in the acid mantle. Deeper in the SC, there is a pH gradient, reflecting increases in deprotonated FFA, with an inflection at the lower third of the SC. At the interior border of the SC, the pH ranges from 6.8-6.917 in which FFAs are completely deprotonated. Additional evidence of deprotonation exists from experimental work on the magainin antimicrobial peptide. This peptide is traditionally used to permeabilize negatively charged bacterial membranes and is able to similarly interact with and permeabilize the SC in the presence of a chemical enhancer19 due to the presence of negatively charged FFAs. Therefore, it is of interest here to investigate the properties of the SC with a deprotonated FFA as models for the interior SC. Previous simulations placing FFA in a PC membrane found that the chain order increased with deprotonated FFA20 demonstrating the effect of the protonation of FFAs on fluid membrane properties. The effect of protonation on solid, highly ordered SC bilayers remains to be seen. The purposes of this study are: (1) to understand the effect of lipid composition, Cer chain length, and temperature affect bilayer properties and possibly affect barrier function; (2) to verify the accuracy of the C36 force field by comparing C36 results against previous experimental results. This work differs from past simulations in force field, FFA protonation state, membrane analyses, and comparison with recent NMR data on the same exact lipid mixtures21. In the past, neutron scattering results22 have been compared against MD15 which is also done in this work. However, the experimental lipid mixtures are more than those used in 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

past and our simulations, and neutron scattering is less detailed than NMR, so comparison with NMR is valuable. The bilayers of this work consist of N-palmitoylsphingosine (Cer16), Nlignoceroylsphingosine (Cer24), Chol, and lignoceric acid (LA). Pure Cer at 32ºC and mixed equimolar Cer/Chol/LA over a range of temperatures are simulated. Mixed bilayers are compared with experimental conditions from Stahlberg et al.21, and a variety of membrane properties are calculated. The results indicate overall good agreement with experiment and reveal how lipid diversity and temperature affect the SC. 2.

Methods

2.1. System setup and simulation protocol Previous experimental studies have examined properties of Cer/Chol/FFA mixtures, usually in equimolar ratios. In our work, we replicated the systems of Stahlberg et al.21 by creating equimolar mixtures of Cer, Chol, and LA at various temperatures, the specifics of which can be found in Table 1. For simplicity in calculation, an excess of one Cer is added to each leaflet, but the concentrations are still roughly equimolar. Pure bilayers of Cer16 and Cer24 were run to investigate the difference between simple and more complex SC models, and these were simulated at 32ºC, the physiological temperature of the SC. All ternary simulations contained 100 lipids per leaflet with 24 water molecules per lipid as estimated from experimental conditions of 50% wt hydration from Stahlberg et al. while pure Cer simulations contained 50 lipids per leaflet and 24 water molecules per lipid. Based on our previous work of ceramide bilayers14, this water:lipid ratio likely produces fully hydrated systems and is useful for saving computational resources. Potassium counter-ions were placed in ternary mixtures to maintain an electroneutral system while pure Cer systems did not contain any ions. Triplicate replicas of each systems were generated in rectangular boxes and tetragonal cells (X = Y ≠ Z) using CHARMMGUI Membrane Builder and equilibrated using the standard six-step process in NAMD23-26. Randomized initial configurations generated by CHARMM-GUI are justified by the agreement between randomized MD and RWMD15. Productions runs were continued for 300-400 ns such that each ternary replicate had at last 300 ns of equilibrated data while pure Cer systems had at least 150 ns of equilibrated data. Equilibration was evaluated through inspection of plots of surface area per lipid (SA/lip) time series and the cumulative average for lack of drift. Table 1 contains the final production time for each system. 4 ACS Paragon Plus Environment

Page 5 of 36 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

Figure 1. a) Chemical structure of Cer16, Cer24, Chol, and LA. Atom labels include terminal tail carbons, tilt vector carbons, RDF selection atoms, and component surface area representative atoms. b) Snapshot of Cer16/Chol/LA at 80ºC at the end of the simulation. Cer16 molecules are

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

shown in cyan with colored headgroups, Chol in blue, LA in red, water as purple dots, and potassium ions as brown beads. Table 1. System parameters for each bilayer system. Number of lipids and water are reported for the whole system. System

# Cer

# Chol

# LA

T (ºC)

# Water

Time (ns)

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

100 68 68 68 68 100 68 68 68 68

0 66 66 66 66 0 66 66 66 66

0 66 66 66 66 0 66 66 66 66

32 32 50 65 80 32 32 50 65 80

2400 4800 4800 4800 4800 2400 4800 4800 4800 4800

300 350 350 400 400 400 350 350 350 350

The CHARMM36 lipid force field27 with updates for sphingolipids13 was used with the TIP3P water model28,

29

. The constant molecule number, pressure, and temperature (NPT)

ensemble was used in NAMD using a constant pressure of 1 bar and with the temperatures as specified in Table 1. Constant pressure was maintained through the Nosé-Hoover Langevinpiston30, 31 and constant temperature was maintained through Langevin dynamics. The LennardJones potential with a force-switching function32 over the range of 10 to 12 Å was used to describe van der Waals interactions. The particle-mesh Ewald (PME) method33 was used to model long-range electrostatics using an interpolation order of 6 and direct space tolerance of 106

with electrostatics being calculated every step. Bond lengths involving hydrogen were kept

constant by the RATTLE algorithm34. The simulations time step was 2 fs and coordinates were saved every 5 ps, except in Cer16/Chol/LA systems at 65ºC and 80ºC which saved coordinates every 10 ps. Visual Molecular Dynamics35 (VMD) was used to create snapshots, and gnuplot36 was used to create plots. 2.2. Analysis Analysis was performed using the last 300 ns for ternary systems and 150 ns for pure Cer systems. All systems had reached equilibrium by the beginning of analysis based on examination of SA/lip plots. The data analyses presented in this work are overall SA/lip, area compressibility modulus (KA), component SA/lip, deuterium order parameters (SCD), Chol tilt, gel phase tilt 6 ACS Paragon Plus Environment

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

angle of pure Cer systems, electron density profiles (EDPs), bilayer thickness, interdigitation, hydrogen bonding (H bonding), lipid clustering, and two-dimensional radial distribution functions (2D RDFs). Most analyses were performed with CHARMM37, 38, 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. Chol tilt angle distribution was calculated by determining the angle between the C3 and C17 (Figure 1) carbons against the bilayer normal, and the average tilt angle was calculated from the first moment of the distribution. For the gel phase tilt angle of Cer, a similar procedure was used except that the lipid tail vector (defined as the C2-C16 vector for the sphingosine chain and the C2-C14 vector for the N-linked chain) was used. The overall SA/lip was calculated by determining the area of the simulation box and dividing by the number of lipids per leaflet. Component SA/lip values were calculated by determining the X and Y coordinates of representative atoms for each lipid (C2S, C1F, and C4S) for Cer lipids, O3 of Chol, and C1 of LA as shown in Figure 1). Using these coordinates, Quickhull39 constructs a Voronoi tessellation in which each atom is represented by a polygon which can be used to calculate individual areas. The sum of all polygonal areas 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,

 =

 〈〉

(2)

 〈〉

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 calculated by recentering the bilayer to place the bilayer center at Z=0 Å, calculating atomic densities with a slab thickness of 0.2 Å over the range of -50 Å to 50 Å, and combining atom densities 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 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

(NSLD) data, SIMtoEXP40 was used to obtain a lipid NSLD profile. From the profile, the peakto-peak distance is calculated to determine the repeat distance. Interdigitation of lipids was calculated by integration of a density overlap profile according to Das et al.10, )+

 = ,+ 4

!" #$%∗!' #

%$(!" #$%)!' #$%*



-.

(3)

where λov is the interdigitation parameter, ρt is the density of the top leaflet, ρb is the density of the bottom leaflet, z is the axis along the bilayer normal, and L is a length that extends past regions of non-zero density (50 Å in this work). Two-dimensional (2D) RDFs were calculated with a ∆r value of 0.1 Å for Cer-Cer, CerChol, and Chol-LA interactions. The C1 atom of Cer, the O3 of Chol, and C1 of LA were considered for all RDFs. The number of intermolecular hydrogen (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. Intramolecular H bonds were not considered as they were previously found to be insignificant in Cer due to the lack of a flexible phosphocholine group14. Lipid clustering was analyzed using a python scikit-learn software41 with a densitybased spatial clustering of applications with noise (DBSCAN) algorithm42. Results from two cutoff distances are reported. A cut-off distance of 4.5 Å between head groups was used for overall clustering while a distance of 5.0 Å was used for Cer-Cer clustering. As overall and Cer-Cer clustering is identical for pure Cer systems, these have two sets of similar clustering data but with different cut-offs. Only cluster sizes of three lipids and more were considered. The C2S atom of Cer, O3 atom of Chol, and C1 atom of LA are considered. All reported p-values are based on the t-test for equal means. 3.

Results

3.1. Overall Surface Area per Lipid (SA/lip) and Area Compressibility Modulus (KA) The overall SA/lip is a useful yet simple measure of lateral packing and lipid density in the bilayer. In particular, a plot of SA/lip over time, shown in Figure S1, is one of the best measures of system equilibration. All systems show relatively quick equilibration in less than 100 ns and plateau until the end of simulation. The probability distribution (Figure S2) follows a normal distribution which justifies the use of the t-test. The average SA/lip (Table 2) drastically 8 ACS Paragon Plus Environment

Page 9 of 36 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

decreases from the pure Cer system to the ternary system, which is a natural result of the small SA/lip of Chol and LA. At 32ºC, the SA/lip between Cer16 and Cer24 are not significantly different in both mixtures and pure systems (p=0.12). However, ternary mixtures of Cer16 and Cer24 exhibit differences at higher temperatures as Cer16 exhibits slightly higher expansions. For both systems, the SA/lip follows a linear thermal expansion even up to 80ºC (Table 2), increasing around 12% from 32ºC to 80ºC. Similar linear trends have been observed for pure PC systems43, and this appears to hold for more diverse lipid mixtures as well. The maintenance of this linear trend is suggestive of a lack of phase transition, which is also apparent from visual inspection (Figure 1) which shows consistent lamellar structure. Table 2. Overall surface area per lipid (SA/lip) and area compressibility modulus (KA) for all systems. System

T (ºC)

SA/lip (Å2)

KA (N/m)

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

43.6±0.3 32.98±0.01 34.06±0.05 35.5±0.1 37.2±0.1 42.8±0.2 32.81±0.01 33.62±0.02 34.9±0.1 36.76±0.03

1.4±0.2 2.9±0.2 1.69±0.07 0.95±0.07 0.64±0.01 2.0±0.2 3.4±0.1 2.36±0.06 1.18±0.03 0.67±0.02

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

The area compressibility modulus (KA), a measure of a membrane’s rigidity, is known to be increased in the presence of Chol44 and Cer14, and Cer has also been referred to as a “Chol surrogate” because of this ability45. In Table 2, the KA is clearly increased in the lipid mixtures over pure Cer systems, which is due to the presence of Chol as Chol is much more rigid than Cer. Both systems have large KA in the N/m range, which is an order of magnitude greater than values for fluid bilayers. Previous simulations using the Berger GROMOS force field by Gupta et al. have shown that bilayers of pure Chol have much higher KA than pure Cer8 although pure Chol bilayers are not physiological. However, between the Berger GROMOS and C36 force fields, there is significant difference in the KA values obtained as the GROMOS bilayer possesses higher rigidity over all temperatures. GROMOS maintains high KA on the order of 3 N/m even at 87ºC while the C36 bilayer is much less rigid with a KA around 0.6-0.7 N/m at 80ºC. A potential factor affecting the rigidity could be the protonation state of LA which is uncharged in the work of Gupta et al. but charged in this work. The proximity of charged LA could inhibit tight packing and increase permeability. On the other hand, the presence of counter ions provides a screening effect that can allow for charge proximity, and deprotonation of FFA in phospholipid bilayers has been found to increase chain order20. Therefore, the difference in KA could be a force field effect. 3.2. Component Surface Areas The component SA/lip is able to account for the natural lower surface area of Chol and LA and is reported in Table 3. In pure Cer, the SA/lip is around 43 Å2, but the SA/lip significantly increases in the mixed system to near 57 Å2 which is related to the shift in H bonding as discussed in section 3.6. Gupta et al.8 noticed a similar trend although the increase was less extreme, changing from 39 Å2 to 46 Å2. In the mixed systems, temperature has the most significant effect on Cer SA/lip as it increases from 57 Å2 at 32ºC to 65 Å2 at 80ºC. Chol and LA also expand, but less extremely than Cer. Chol has SA/lip in the 25-30 Å2 range, which is standard for sterols, while LA has SA/lip in the 15-18 Å2 range. This is significantly lower than the SA/lip of Chol, which is somewhat expected given the bulkiness of the sterol ring compared to the narrow fatty acid chain. In general, the SA/lip is similar between Cer16 and Cer24, suggesting that the effect of chain length is minimal in comparison to the headgroup identity, 10 ACS Paragon Plus Environment

Page 11 of 36 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

which is not surprising given that headgroup atoms are used to represent the lipid for analysis. The same applies to Chol and LA which do not show significant differences between Cer16 and Cer24 mixtures. Table 3. SA/lip (Å2) of specific lipid types in all systems. System

T (ºC)

Cer (Å2)

Chol (Å2)

LA (Å2)

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

43.6±0.3 57.3±0.6 60.0±0.3 61.9±0.2 65.3±0.2 42.8±0.2 57.3±0.2 58.8±0.3 61.4±0.4 64.4±0.3

-25.4±0.3 25.4±0.2 27.1±0.1 28.3±0.2 -24.6±0.2 25.2±0.1 26.0±0.2 27.64±0.03

-15.9±0.2 16.1±0.1 16.46±0.04 17.26±0.04 -16.05±0.03 16.1±0.1 16.7±0.1 17.40±0.04

3.3. Deuterium Order Parameters (SCD) Previous work has compared C36 deuterium order parameters (SCD) to values obtained from NMR for Cer/PC mixtures and found general good agreement14, although discrepancies appear near the melting temperature as C36 inaccurately describes the phase transition which is significantly blurred by the introduction of Cer. As the systems of this work do not transition through a melting temperature, this problem is unlikely to appear. Furthermore, NMR order parameters and phase composition estimates from Stahlberg et al.21 are available for comparison and discussion. Due to the lack of site-specific deuteration, order parameters are sorted to be monotonically decreasing for comparison purposes.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

Figure 2. Comparison of C36 SCD against NMR for ternary mixtures for the a) fatty acid chain of Cer16, b) fatty acid chain of Cer24, and c) fatty acid chain of LA. Values are sorted to be monotonically decreasing. The C36 SCD is in overall good agreement with NMR, particularly for Cer16 and LA at 65º and 80ºC (Figure 2). However, there is a slight discrepancy in the SCD for Cer24 at 65ºC, which shows a region of greater disorder than suggested by NMR. There are several possible factors for this discrepancy. One possibility is that the NMR signal reported by Stahlberg et al. does not display prominent doublets for Cer24 at 65ºC, drawing the validity of the fitted SCD into question. That said, since the SCD results for Cer16 and LA likewise do not display clear doublets and agree well with C36 data, it would be inconsistent to disregard the results for Cer24. From 12 ACS Paragon Plus Environment

Page 13 of 36 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

NMR, the order of Cer24’s fatty acid tail is similar to LA’s order, which is sensible given the similarity in structure. Comparing SCD between Cer24 and LA for C36, the difference lies after the C10 carbon. Whereas Cer24 exhibits a near linear decrease after C10 and reaches a value of 0.30 at C15, LA exhibits a much shallower decrease in the C10-C15 interval and reaches a value of 0.34 at C15 (Figure S3). The order of Cer24 is low in this range and causes the discrepancy with NMR. Unfortunately, SCD comparison is not possible at physiologically relevant temperatures due to a lack of reported experimental data, and the large system sizes necessary for a complete phase composition analysis are not feasible for simulation. Previous experimental work21 estimated the phase composition of their samples at different temperatures for each system by fitting of free induction decays to time domain data of the NMR spectrum. Three phases were considered in their model: gel/crystalline, fluid, and isotropic. By calculating a superposition of two of three line shapes which are each characteristic of one of these phases (Pake doublets for gel/crystalline/fluid and a single line for isotropic), the proportion of each phase was calculated. These phases indicate varying levels of disorder with gel/crystalline phases being highly ordered and exhibiting little lipid mobility, fluid phases exhibiting faster lipid motion while retaining lamellar structure, and isotropic phases indicating loss of lamellar structure. From C36, the lamellar structure of these systems at 32ºC (Figure S4) qualitatively agrees with a nearly fully crystalline system from the experimental model. However, the maintenance of this lamellar structure at 80ºC contradicts the isotropic model from experiment. Both C36 and NMR report LA SCD around 0.35 at 80ºC which is inconsistent with a highly disordered isotropic structure suggested by the experimental model fit. A tail comparison between all systems is given in Figure S5. Both ternary mixtures of Cer16 and Cer24 show the expected trend of decreasing order with temperature. At 32ºC, both fatty acid and sphingosine chains reach peak values of around 0.45, which is indicative of either the gel or solid-ordered state. Given the significant sterol concentration in the bilayer as well as the lack of chain tilt, the mixture systems are likely in the solid-ordered state while the pure Cer systems are in the gel state. At 80ºC, the values peak at around 0.35 which is still relatively high and bordering the liquid-ordered region. Between Cer16 and Cer24, the order is nearly identical over matching carbons. However, past the C13 carbon, the order of Cer16 quickly drops while 13 ACS Paragon Plus Environment

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

Page 14 of 36

Cer24 more steadily drops. There is a relatively sharp edge between the hydrophobic region and the free volume in the bilayer center for Cer16, which leads to the sharp drop in order as only the terminal carbons are free to rotate. However, due to chain asymmetry in Cer24, this sharp edge is replaced by a region of fluid interdigitating chains that leads to the steady decrease in SCD. Between mixed and pure systems, there is a clear increase in order for mixed systems which is due to the ordering effect of Chol that is stronger than in Cer. This correlates with higher compressibility moduli as previously mentioned but seems to conflict with an increase in component SA/lip. As has been previously observed in sphingolipid simulations14,

44

, an

expansion of SA/lip while increasing in chain order is possible if the relevant H bonds concurrently decrease. The order of LA is relatively similar to the order of the fatty acid tail in Cer, although there is greater order of C15-C24 carbons for LA mixed with Cer24 due to the additional chain length restricting bond rotation. 3.4. Tilt Angle of Chol and Gel Phase Chains Chol possesses a natural, slight tilt when oriented in the bilayer which is easily characterized due to its rigid ring structure. The average tilt angle (Table 4) shows a clear dependence on temperature which is correlated with decreasing chain order and greater SA/lip. However, unlike in previous simulations of SM/Chol systems44, an increase in tilt at higher temperatures is not correlated with greater Chol H bonding, so the ultimate mechanism of tilt likely lies in SA/lip and chain order. Naturally, Chol with greater area will have greater tilt in order to occupy that area, and disordered Cer chains allow more space for tilting. Additionally, the greater order in Cer24 at the midplane carbons restricts the tilting of Chol so that it remains more upright compared to systems containing Cer16. In all systems, the tilt angle is small with the largest value being Cer16/Chol/LA at 80ºC with an angle of 13.5º. Therefore, all the mixtures systems have a general lack of strong lipid tilt in the bilayer which supports existence of the solid-ordered phase.

14 ACS Paragon Plus Environment

Page 15 of 36 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

Table 4. Average tilt angle for Chol in mixed systems and Cer chains in pure Cer. Cer chain tilt is calculated as an average between fatty acid and sphingosine chains. System

T (ºC)

Cer (º)

Chol (º)

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

17±1 ----18.5±0.6 -----

-9.0±0.1 10.2±0.1 11.6±0.1 13.5±0.1 -8.9±0.1 9.4±0.1 10.85±0.03 12.80±0.02

Upon visual inspection of the pure Cer systems, it is apparent that unlike the mixed systems, these have more bilayer tilt and exist in the tilted gel phase. However, there are differences between replicates which produce differences in chain tilt, and a comparison is found in Figure S6. As has been previously found with pure Cer systems in C3614, there are disordered lipids within the bilayer that intrude further. Occasionally, ordered lipids protrude farther which can produce a buckling effect, best demonstrated in replica 2 of Cer16. A more detailed discussion of similar systems has been provided in previous work14. Ultimately, the difference between Cer16 and Cer24 is not significant, and both exhibit chain tilt in the 16-18º range. This agrees with previous work on pure Cer systems14, but the low tilt in Cer24 is somewhat unexpected as Cer24 is experimentally determined to have a tilt angle of 29º46. The most likely reason for this discrepancy is the non-lamellar structure that Cer24 takes in crystallography. Due to the inverted cone shape of Cer and low hydration, deviation from the traditional hairpin conformation is feasible and crystallography tends to capture V-shaped packing of Cer24. Since the Cer24 is simulated in a bilayer, comparison of the experimental tilt angle to the simulated tilt angle produces the aforementioned discrepancy. The propensity for V-shaped packing is also a possible cause of the bilayer defects and bilayer buckling that was previously found. Since the bilayer phase is not completely stable, it seems that pure Cer systems are inadequate models of lamellar SC phases. At the same time, it has long been known that the SC is incredibly heterogeneous and that pure Cer is simply a starting model. Taken with the possibility that pure Cer exhibits a mix of phases, some of which may be lamellar, systems of pure Cer still have a 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 16 of 36

place in SC modeling as long as there is an understanding the more complex models will eventually need to take their place. 3.5. Electron Density Profiles (EDPs), Bilayer Thickness, and Interdigitation Groen et al.22 have constructed neutron scattering length density (NSLD) profiles for SPP models using Cer, Chol, and FFA mixtures at pH 5. For comparison purposes, the NSLD profile of Cer24/Chol/LA at 32ºC is given in Figure 3. Despite some discrepancies, the general shape of the NSLD profile is in agreement with experiment. In particular, the appearance of secondary peaks at ±15 Å due to the presence of Chol and a slight bump in the bilayer center are consistent between C36 and experiment. However, the height of the peaks is more accentuated in C36, and the major peaks representing the SPP repeat distance are more closely spaced. The C36 peak-topeak distance is 4.9 nm, which is less than the 5.36 nm from experiment. A number of factors could cause a thinner bilayer. Among these, temperature and hydration differences are less likely to be significant since the temperature difference is less than 10ºC and hydration is suggested by Groen et al. to have minimal effect on repeat distance. The most likely cause is the diversity in Cer type as Groen et al. include a variety of Cer which vary in hydroxylation. These include Nlignoceroylsphingosine

(NS-18:1/24:0),

N-lignocerophytosphingosine

(NP-18:0/24:0),

α-

hydroxy-lignoceroylsphingosine (AS-18:1/24:0), N-palmitoylphytosphingosine (NP-18:0/16:0), and α-hydroxy-lignoceroylphytosphingosine (AP-18:0/24:0) in a 60:19:5:11:6 molar ratio. NS bases have two hydroxyls while NP, AS, and NP bases have three hydroxyls and AP bases have four hydroxyls. Cer in this work are exclusively NS bases, so there is less potential for H bonding which decreases chain order and bilayer thickness.

16 ACS Paragon Plus Environment

Page 17 of 36 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

Figure 3. Comparison of NSLD profiles obtained from C36 for Cer24/Chol/LA at 32ºC and experimental neutron scattering of the SPP. Since Cer24/Chol/LA at 32ºC is the most physiologically accurate system studied in this work, it is the only system that is analyzed for an NSLD profile. Further analysis is done through examination of EDP profiles and calculated bilayer thicknesses. The overall EDP (Figure 4) shows the strongest change in the pure Cer to ternary mixture transition, particularly for Cer16. Pure Cer16 displays closer major peaks, lack of secondary peaks at ±10 Å, and lower density in the bilayer center compared to Cer16 mixtures. Besides the lack of secondary peaks, which are due to the presence of Chol, pure Cer24 exhibits these differences to much lesser degrees than pure Cer16. The lack of interdigitating chains in pure Cer16 (Figure S6) is likely responsible for these differences including the closer major peaks that result from near symmetric Cer chains allowing for tighter packing. Cer24, which has significant interdigitation in all systems, is less affected by the introduction of LA. However, there is a larger bump in the bilayer center that is due to higher interdigitation. In the ternary mixtures, temperature causes consistent effects of lower density and bilayer thinning while maintaining the overall shape. This agrees with the small linear increases in SA/lip previously mentioned.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

Figure 4. Total electron density profiles of pure Cer at 32ºC and mixed systems at various temperatures for a) Cer16 systems and b) Cer24 systems. Component EDPs are given in Figures S6-7. The Cer EDP in Figure S7 displays significant differences between pure Cer, 32ºC mixtures, and 80ºC mixtures. In pure Cer16, the methylene density significantly drops at the bilayer center and is near the methyl density due to a lack of interdigitation. This is also responsible for the drop in the overall EDP previously mentioned. The gap between methylene and methyl densities in the bilayer center is larger for Cer24 systems and Cer16 mixtures due to interdigitation. Additionally, the sphingosine density shifts outward in Cer16 mixtures over pure Cer which is indicative of bilayer thickening. A shift in water density toward the bilayer center is apparent in ternary mixtures which represents greater water penetration. Whereas the water density mainly exists exterior to the sphingosine peak in pure Cer, it crosses near the peak in mixtures. In Cer24/Chol/LA at 32ºC, the methylene and methyl densities demonstrate the appearance of secondary peaks which have previously been 18 ACS Paragon Plus Environment

Page 19 of 36 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

seen in mixtures containing Chol44. These represent a retardation of orthogonal movement due to rigidifying presence of Chol. Cer16/Chol/LA also exhibits these peaks but to a somewhat lesser extent than Cer24. Unsurprisingly, these peaks disappear at 80ºC as disorder is enhanced. The peaks are replaced by regions of plateauing density, especially in Cer24/Chol/LA. Chol/LA densities, shown in Figure S8, confirm that the ±10 Å peaks are due to the presence of Chol. Similarly to Cer, Chol possesses EDPs with secondary peaks at 32ºC that smooth out at 80ºC, although these are more subtle than in Cer. Chol in Cer16 reveals greater density in the bilayer center compared to Chol in Cer24 due to greater interdigitation of Chol in Cer16. The carboxyl of LA is always located exterior to the Chol hydroxyl, suggesting that it is located nearer to Cer headgroups due to the longer fatty acid chain. The bilayer thicknesses (DHH, DB, 2DC) provide quantitative justification for trends seen in EDPs (Table 5). The most notable effect is the significantly smaller DHH in pure Cer16 compared to Cer16 mixtures which agrees with the overall EDPs. The DHH of Cer24 is significantly thinner in mixed systems over pure (p=0.03), and Cer24 is always thicker than Cer16 due to the longer fatty acid chain. Interestingly, although LA has significant effects on thickness in the absence of Cer24, Cer24 is the more dominant factor despite both Cer24 and LA having fatty acid chains that are 24 carbons long. Temperature induces the expected effect of bilayer thinning over all thicknesses, and it also induces higher water penetration which can be measured as DB-DHH. DB-DHH consistently decreases with temperature, suggesting deeper penetration into the bilayer. In the case of Cer16/Chol/LA at 80ºC, DB-DHH is even slightly negative (-0.3 Å) which is unusual for most bilayers. The effect of the lipid mixture also contributes to this effect because the Chol headgroups are localized deeper in the membrane. This provides favorable interactions for water to penetrate deeper compared to pure Cer. However, since Chol also significantly orders the bilayer, the hydrophobic region is thickened (2DC increases by 5 Å) which inhibits water permeation. Along with 2DC, DHH also increases by 5 Å which indicates that this is mainly due to a propagation of the thickening of the hydrophobic region. For systems of Cer24, 2DC and DHH remain constant while DB decreases, so the ordering effect of Chol is outweighed by its deeper location in the bilayer.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

Table 5. 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. System

T (ºC)

DHH

DB

2DC

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

39.4±0.3 44.7±0.1 44.5±0.2 43.9±0.2 43.0±0.0 50.5±0.4 49.1±0.1 48.7±0.1 48.1±0.1 46.8±0.0

43.9±0.4 45.64±0.02 45.0±0.1 44.0±0.1 42.7±0.1 54.2±0.3 50.26±0.03 49.87±0.01 48.9±0.1 47.32±0.03

35.4±0.4 40.0±0.1 39.43±0.04 38.7±0.1 37.68±0.04 45.1±0.3 44.4±0.1 44.21±0.03 43.54±0.03 42.31±0.03

The interdigitation of each lipid is reported in Table 6. The values are obtained using an integration of the density profile according to Eq. 3. Regions of no overlapping density between leaflets will have values of 0 due to the product term while regions of complete overlap will have values of 1 due to the numerical prefactor. Integration of this distribution can be thought of as the average length of lipid that is interdigitating, so larger values in Table 6 indicate greater interdigitation. Most notably, the interdigitation of Cer increases 3-4 fold from Cer16 to Cer24 in mixture systems and 2 fold in pure Cer. This is natural given the chain asymmetry in Cer24, but the strong difference between pure and mixed Cer16 is somewhat unexpected. The strong decrease in interdigitation is due to the increase of chain ordering and bilayer thickness in mixed systems. Since the sphingosine and fatty acid chains are similar in length for Cer16, the chain order increases throughout the whole fatty acid chain which decreases interdigitation, but the asymmetry in Cer24 produces little difference in order for the terminal carbons where interdigitation is most prominent. In all systems, Chol has relatively low interdigitation and LA has higher interdigitation as expected given the long chain structure. LA interdigitates less in Cer24 due to competition with the long Cer24 chains which preferentially interdigitate over LA. The natural preference of Cer24 over LA explains why Cer24 is more influential on bilayer thickness as previously mentioned. Similar to LA, Chol generally interdigitates less in Cer24. Overall, the effect of temperature on interdigitation is not straightforward as it can increase or decrease in different systems. For example, LA and Cer in Cer16 interdigitate more above 32ºC 20 ACS Paragon Plus Environment

Page 21 of 36 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

but interdigitate less above 32ºC in Cer24. High temperatures induce greater interdigitation in Cer16 and less interdigitation in Cer24 as the natural tendencies to interdigitate are opposed. LA, due to its fluidity, is carried along by the Cer and follows the same trends for increasing or decreasing interdigitation. Table 6. Interdigitation (Å) of individual lipids for each system. System

T (ºC)

Cer (Å)

Chol (Å)

LA (Å)

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

3.9±0.3 1.75±0.02 1.81±0.05 2.09±0.02 2.37±0.03 7.64±0.09 7.46±0.04 7.02±0.04 7.08±0.04 7.21±0.04

-2.1±0.1 2.23±0.09 2.43±0.02 3.2±0.2 -1.01±0.08 1.17±0.03 1.4±0.04 1.74±0.01

-8.4±0.1 8.4±0.1 8.51±0.03 8.81±0.02 -7.21±0.04 6.7±0.2 6.65±0.07 6.71±0.02

3.6. Hydrogen Bonding (H bonding) Cer are notable for possessing high H-bonding potential due to the presence of both donor and acceptor groups on the sphingosine backbone. The overall H-bonding frequency is reported in Table 7. In general, Cer has the most H-bonding as the Cer frequency ranges from 0.3-0.4 while Chol and LA fluctuate around 0.2. This is not surprising given the greater number of H-bonding groups on Cer compared to Chol and LA. There is a clear increase in Cer H bonding when transitioning from pure Cer to mixed systems due to the presence of H-bonding groups on the low SA/lip Chol and LA. This rise in H bonding contributes to the rise in order parameters and compressibility moduli in ternary mixtures, although there is also a rise in SA/lip. This increase suggests that the type of bonding atoms, discussed in the next paragraph, is very influential on SA/lip. For Cer H bonding, the frequency is not very dependent on chain length which is somewhat expected as the headgroups are identical, but it does decrease with temperature, around 10% from 32ºC to 80ºC. LA shows an even steeper 30% decrease while Chol shows a weaker 5% decrease which could be related to the greater fluidity of LA making it more susceptible to temperature effects. Overall, the transition from 65ºC to 80ºC decreases H bonding the most, which is likely due to the drop in KA. Although the drop itself is not very

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

large, the KA at 80ºC is around 0.6 N/m which is approaching values for fluid bilayers (0.3 N/m) and makes differences more profound. Table 7. 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. Lipid

System

T (ºC)

NHB/lip

Cer

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80 32 32 50 65 80 32 32 50 65 80 32 32 50 65 80 32 32 50 65 80

0.350±0.003 0.44±0.02 0.413±0.001 0.414±0.003 0.369±0.001 0.380±0.006 0.42±0.03 0.43±0.02 0.40±0.01 0.371±0.004 -0.192±0.006 0.21±0.01 0.194±0.002 0.186±0.003 -0.195±0.008 0.192±0.002 0.199±0.002 0.186±0.001 -0.22±0.02 0.176±0.004 0.171±0.002 0.158±0.001 -0.224±0.01 0.218±0.007 0.185±0.003 0.1756±0.0003

Chol

LA

Individual bond frequencies for the most significant bonds are reported in Tables S1-2. As is common in sphingolipids, the amide-carbonyl bond is the most significant bond with high 22 ACS Paragon Plus Environment

Page 23 of 36 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

frequencies in the 0.3-0.4 range for pure Cer. However, due to competition with bonding groups on Chol and LA, the frequency drops significant in mixed systems to around 0.15. While the overall Cer H bonding increases due to the addition of other bonds, the loss of the amidecarbonyl bond has significant effects on Cer properties. As previously mentioned, the Cer SA/lip expands in ternary mixtures. As Cer SA/lip is especially sensitive to Cer-Cer H bonding, the loss of the amide-carbonyl bond, which is by far the most significant Cer-Cer bond, allows an expansion of SA/lip. Interestingly, temperature can cause an increase in Cer-Cer bonding from 32ºC to 65ºC. Since the overall Cer bonding either stagnates or decreases, this is due to a shifting of bonds toward Cer-Cer bonds from Cer-LA bonds. Bonds involving LA are evenly split between the two oxygen atoms. This is expected given that the two oxygens are identical in a resonance hybrid, but it is still fortunate since MD does not fully capture the resonance phenomenon. Chol H bonding is not sensitive to temperature and is generally slightly weaker than Cer H bonding. Chol-Cer bonds have frequencies around 0.1, and Chol-LA bonds have frequencies around 0.04. The weaker bonding could be due to improper alignment as Chol is located deeper in the membrane than Cer and LA.

3.7. Lipid Clustering and 2D Radial Distribution Functions (2D RDFs) Another aspect of lateral packing is the propensity for lipid clustering, which is closely related to H bonding. The number of clustered lipids over time is given in Figure S9 and demonstrates convergence of clustering. The overall number of clustered lipids (NLT) and number of clusters (NCT) (Table 8) shows a dependence on temperature as both decrease at higher temperatures. This agrees with an increase in SA/lip and decrease in H bonding. Since NLT and NCT decrease approximately proportionately, the size of the clusters does not significantly change with temperature. There is an increase in both NLT and NCT from pure to mixed systems which is due to the naturally small headgroups of Chol and LA (Figure 5). Although a cutoff distance of 5 Å is more appropriate for cluster analysis with larger headgroups, this produces errors in which the clusters become too large and lipids become duplicated in single clusters. Therefore, a cutoff distance of 4.5 Å is used for overall analysis, and a distance of 5 Å is used for Cer-only clustering. The results of Cer-only clustering demonstrate more clustered lipids in pure Cer due to the larger cutoff and fewer lipids in ternary mixtures due to the 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

lower Cer concentration. However, Cer-only clustering recapitulates the phenomenon of H bond shifting previously mentioned. The number of clustered Cer (NLCER) and number of Cer clusters (NCCER) increase in ternary mixtures from 32ºC to 65ºC which is generally unexpected for an increase in temperature. NCCER is most affected by Cer-Cer H bonding, and there is an increase in amide-carbonyl bonding between these systems which drives the increase in clustering. The lipid ratios in clusters in reported in Table S3. The number of Cer per cluster decreases from 3.54.0 in pure systems to around 1.5 in mixed systems, which is again due to competition with Chol and LA. There is little effect on lipid ratios due to temperature differences, but Cer is generally overrepresented in clusters relative to the bulk concentration while LA is underrepresented. Chol nearly perfectly matches the bulk concentration. Table 8. Average number of clustered lipids (NLT) and number of clusters (NCT) for overall clustering with a cutoff distance of 4.5 Å. Number of clustered lipids (NLCER) and number of clusters (NCCER) for Cer-only clustering with a cutoff distance of 5.0 Å is also reported. Cluster sizes of 3 lipids and greater were considered. System

T (ºC)

NLT

NLCER

NCT

NCCER

Cer16 Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer16/Chol/LA Cer24 Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA Cer24/Chol/LA

32 32 50 65 80 32 32 50 65 80

35.0±0.2 48.3±0.4 47.0±0.5 42.9±0.3 39.9±0.3 40±2 50.1±0.7 46.9±0.7 44.3±0.3 40.2±0.1

60±1 7.2±0.4 7.47±0.03 8.3±0.6 7.4±0.1 66±2 6.9±0.5 7.8±0.6 7.4±0.3 7.2±0.1

9.52±0.06 12.0±0.2 11.9±0.1 11.04±0.05 10.49±0.05 10.44±0.06 12.6±0.2 11.8±0.1 11.3±0.1 10.51±0.02

13.4±0.1 2.05±0.06 2.18±0.01 2.4±0.2 2.16±0.03 13.5±0.2 2.1±0.1 2.3±0.1 2.18±0.08 2.12±0.03

24 ACS Paragon Plus Environment

Page 25 of 36 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

Figure 5. Visual representation of lipid clustering for top leaflet in a) pure Cer16 and b) Cer16/Chol/LA at 32ºC. Yellow circles represent Cer lipids, blue circles represent Chol, and red circles represent LA. A cutoff distance of 4.5 Å is used. Clusters are taken from the final snapshot of the simulation rather than from a time average. RDFs offer another perspective for clustering analysis. The evolution of RDFs from the initial configuration is demonstrated in Figure S10. The strongest clustering is Cer-Cer clustering which reaches maximums near 2.0 (Figure 6). Additionally, there is a slight increase for systems at 65ºC which restates the increase in Cer clustering. This peak is located at a radial distance of 5.0-5.2 Å which justifies the previously used cutoff distance of 5.0 Å for clustering analysis. Past 10 Å, clustering disappears as the RDF converges to the bulk. Chol-Chol RDFs reveal the formation of repeating, evenly spaced peaks that are indicate the formation of up to three Chol shells. The individual peaks are not as strong as in Cer-Cer RDFs as these reach maximums of 1.4. The formation of these Chol shells is common in systems with significant concentrations of 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

Chol and usually have repeat distances in the 5.0-5.5 Å range47, 48. In this work, the peaks are located at distances of 6.0, 10.8, and 15.9 Å. The first peak is slightly further spaced compared to previous systems, perhaps due to the lack of crowding phosphates, but is not drastically different. LA-LA RDFs (Figure S11) are weaker than both Cer-Cer and Chol-Chol RDFs, reaching peaks values of near 1.2, which agrees with the relatively weak clustering. After the primary peak, the RDF quickly drops and converges to the bulk without dipping, suggesting the lack of an empty region after the first peak. RDFs between different lipids are also reported in Figure S12-14. The strongest peak is in Chol-LA RDFs which reach maximums at 1.6, suggesting that these interactions could be the most significant for clustering. After the primary peak, secondary peaks also appear similarly to the Chol-Chol RDFs which form due to the Chol shells. In contrast, CerLA and Cer-Chol interactions never far exceed the bulk and are likely less significant than CholLA interactions.

Figure 6. Comparison of 2D RDFs for pure Cer at 32ºC and mixed systems at varying temperatures. (a,b) represent Cer-Cer (C2S-C2S) RDFs and (c,d) represent Chol-Chol (O3-O3) RDFs. (a,c) represent systems involving Cer16 and (b,d) represent systems involving Cer24. Lighter shades of each color represent ±standard error.

26 ACS Paragon Plus Environment

Page 27 of 36 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

4.

Discussion and Conclusion Due to the difficulty and cost of simulating the extremely heterogeneous lipid matrix of

the SC, especially of the multilayered LPP, most atomistic simulations rely on a simple bilayer model of the SC to understand interactions between Cer, Chol, and LA. Bilayers are expected to best model the SPP of the SC. Gupta et al. provided a systematic study of Cer24/Chol/LA with a protonated LA over a variety of lipid compositions and temperatures8 and also studied chain length effects in pure Cer9 using the Berger GROMOS force field. In mixed systems, they found greater stability at high temperatures at high temperatures and greater Cer SA/lip in systems with Chol over pure Cer. These are in qualitative agreement with the results of this work, but the SA/lip are greater in C36 even for pure Cer systems as Cer2 has a SA/lip of 39.3 Å2 in GROMOS but 42.8±0.2 Å2 in C36. Additionally, the effect of transitioning to mixed systems has different effects on Cer chain order. In C36, this increases the SCD from 0.35 to 0.45 while it decreases from 0.45 to 0.42 in GROMOS (SZ was converted to SCD by multiplying by a ½ factor). The significant difference thus lies in the order of the gel phase pure Cer systems. Although C36 systems are evidently in a gel phase, the calculated order parameters are lower than GROMOS which could be due to the presence of residual disordered Cer and ultimately related to instability of the inverted cone shape. Das et al.10 simulated pure and mixed systems, similar to the systems of this work, using GROMOS and found the existence of a fluid like environment in the bilayer center which agrees with both the lowered SCD and snapshots of this work. A discrepancy worth noting arises in the determination of SC lipid composition. Das et al. assert that the physiological composition is a Cer/FFA/Chol ratio of 2:2:1, but X-ray results from Gooris and Bouwstra49 found that the SC is approximately equimolar. Since the 2:2:1 ratio reported by Das et al. is not empirically determined and is uncited, the equimolar ratio determined by Gooris and Bouwstra is likely the more accurate estimate. Of course, the SC is very heterogeneous so local concentrations are likely to vary. For example, Chol is known to be asymmetrically distributed50. Notman et al.51 have used the GROMOS force field to study the effect of DMSO on SC bilayers. We have attempted to replicate these effects in C36 but have found an inability for either DMSO or ethanol to penetrate highly ordered bilayers at physiological temperatures (unpublished work). This could be a result of C36 DMSO52 being parameterized with respect to its interaction with

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

water and not bilayers. Additionally, the lack of polarizability could be detrimental to the simulation of permeating molecules. Moore et al.15 have used the C36 force field to simulate similar systems of Cer/Chol/LA bilayers, but with differing sphingolipid parameters and protonated LA. Additionally, they used a RWMD procedure which allows for efficient sampling with highly biased initial builds. From a randomized initial build, the results of MD and RWMD are very similar which supports the validity of RWMD. In this work, the bilayer build process was randomized through the use of CHARMM-GUI Membrane Builder, so our MD results are also unbiased. In general, Moore et al. noticed little change in lateral membrane properties which agrees with the results of our work. This is unsurprising given that lateral membrane properties are dominated by headgroup interactions, and chain properties such as interdigitation are more significantly affected by chain length. They also characterized the tilt angle of pure Cer16 and found the average angle to be around (16±2)º, similar to our value of (16.5±1.1)º. NSLD calculation and comparison to that experimentally obtained by Groen et al.22 was also performed and similarly found a thinner bilayer. However, their NSLD differs as the secondary peaks due to Chol are less distinct, appearing only as shoulders, and the bump in the bilayer center is more accentuated. In our work, the NSLD displays highly distinct secondary peaks and a weaker bump in the center. Although the secondary peaks are too large, the center bump more closely agrees with experiment. Since the NSLD profiles from experiment are based on Fourier fits to form factors, there could also be artifacts in the shape features in the profile so direct comparison should be taken with caution. Moreover, several Cer lipids were used by Groen et al.22 in their Cer/Chol/LA experiments and likely cause the peak-to-peak distances to be greater with the Cer lipids having more hydroxylation. Overall, we argue that the NSLD of our work, using C36 force field parameters from Venable et al.13 and deprotonated LA, is more accurate than that of Moore et al. due to the clear presence of secondary peaks and a weak center bump. Comparison between MD and experiment can be difficult for SC mixtures as simulation boxes are too small to study large effects such as domain segregation and phase separation from experiment. More general properties have been suggested by Mizushima et al.53 who used differential scanning calorimetry to study Cer/Chol/FFA mixtures and found that Chol tightens lateral packing and decreases lipid chain mobility. These agree with lowered SA/lip and greater 28 ACS Paragon Plus Environment

Page 29 of 36 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

order parameters from pure Cer to mixed systems in this work. Using NMR, Bjorklund et al.54 confirmed that heating increases lipid mobility which is in agreement with lowered order parameters at higher temperatures in mixed systems. Lipid dynamics was not probed in this work, but previous results in SM/Chol bilayers suggested that relaxation times of cross-chain correlation functions exceed the hundreds of nanoseconds range which is computationally expensive and exceeds the timescales of this work. Stahlberg et al.21 constructed models of phase compositions in Cer/Chol/LA mixtures at various temperatures as well as calculating order parameters. They found that the double bond in the sphingosine chain of Cer is essential to maintaining the crystalline phase which is likely related to the propensity for the double bond to allow for orthorhombic chain packing. Phytosphingosine Cer contain greater hydroxylation and can likely facilitate more ordered structure through H bonding rather than through physical packing. At 50ºC and 65ºC, they generally found mixes of gel and fluid phases while the vast majority of lipids exist in an isotropic phase at 80ºC. While the results of this work also indicate the greatest changes in the transition from 65ºC to 80ºC, a lamellar structure was always maintained. The isotropic structure determined from experiment is inconsistent with both the calculated SCD and all SC bilayer simulations, including those of this work. As the pH was not fixed in experiment, LA was most likely deprotonated with the potential to protonate with pH drift. If not due to LA protonation, the isotropic phase discrepancy could be due to invalid assumptions associated with the phase composition calculation. Alternatively, if the composition calculations are valid, isotropic structure could be attributed to inverse micelles which have been simulated by Das et al.55 and are consistent with non-lamellar pockets of the SC. As the systems of this work did not see any inverse micelle or inverse hexagonal phases, future work might investigate the effect of temperature and FFA protonation state on the formation of these non-lamellar phases. Furthermore, complete exploration of the high temperature systems of Stahlberg et al. could require a random arrangement of lipids to counter kinetic trapping effects of the initial bilayer build which can bias results toward lamellar phases. The systems studied in this work are pure Cer at 32ºC and Cer/Chol/LA bilayers over a range of temperatures using the C36 force field. The pure-mixed transition has strong effects on interdigitation in Cer16 systems but not Cer24 systems. Additionally, a unique shift in H bonding 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

is found at 65ºC. The overall H bonding does not change, but Cer-LA H bonding decreases and leads to a corresponding increase in Cer-Cer H bonding which also affects lipid clustering. Deuterium order parameters are compared against experimental NMR values and show strong agreement save for Cer24 at 65ºC, which is attributed to decreased chain order after the C10 carbon of the fatty acid chain based on comparison with LA. Based on NSLD data, the qualitative shape of the profile is in general agreement, but the calculated thickness is underestimated in C36 which is either due to inaccuracy of the force field, possibly related to the decreased order in Cer24, or more likely the presence of varying Cer types in experiment. Nevertheless, strong order parameter agreement in Cer16 and LA, along with previous work suggesting greater accuracy in C36 over GROMOS for Cer bilayers14, supports the use of C36 for further simulations of mixed bilayers and even multilayers as models of the SC. 5.

Supporting Information Available Figures of SA/lip as a function of time, SA/lip probability distribution, order parameter

plots, snapshots of pure Cer and mixture bilayers, component electron density plots, clustering over time, and radial distribution functions. Tables of hydrogen bond frequency and clustering fractions. 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 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 Stampede2 and the computational resources at the Maryland Advanced Research Computing Center (MARCC). The authors would also like to thank Dr. Daniel Huster in providing more detail on pH conditions of his experimental measures.

30 ACS Paragon Plus Environment

Page 31 of 36 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

7.

References

1.

Mojumdar, E. H.; Gooris, G. S.; Groen, D.; Barlow, D. J.; Lawrence, M. J.; Deme, B.; Bouwstra,

J. A., Stratum Corneum Lipid Matrix: Location of Acyl Ceramide and Cholesterol in the Unit Cell of the Long Periodicity Phase. Biochimica Et Biophysica Acta-Biomembranes 2016, 1858 (8), 1926-1934. 2.

Melnik, B.; Hollmann, J.; Hofmann, U.; Yuh, M. S.; Plewig, G., Lipid-Composition of Outer

Stratum-Corneum and Nails in Atopic and Control Subjects. Archives of Dermatological Research 1990, 282 (8), 549-551. 3.

Bouwstra, J. A.; Ponec, M., The Skin Barrier in Healthy and Diseased State. Biochimica Et

Biophysica Acta-Biomembranes 2006, 1758 (12), 2080-2095. 4.

Breathnach, A. S.; Goodman, T.; Stolinski, C.; Gross, M., Freeze-Fracture Replication of Cells of

Stratum Corneum of Human Epidermis. Journal of Anatomy 1973, 114 (JAN), 65-81. 5.

Swartzendruber, D. C.; Wertz, P. W.; Kitko, D. J.; Madison, K. C.; Downing, D. T., Molecular-

Models of the Intercellular Lipid Lamellae in Mammalian Stratum Corneum. Journal of Investigative Dermatology 1989, 92 (2), 251-257. 6.

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. 7.

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. 8.

Gupta, R.; Rai, B., Molecular Dynamics Simulation Study of Skin Lipids: Effects of the Molar

Ratio of Individual Components over a Wide Temperature Range. Journal of Physical Chemistry B 2015, 119 (35), 11643-11655. 9.

Gupta, R.; Dwadasi, B. S.; Rai, B., Molecular Dynamics Simulation of Skin Lipids: Effect of

Ceramide Chain Lengths on Bilayer Properties. Journal of Physical Chemistry B 2016, 120 (49), 1253612546. 10.

Das, C.; Noro, M. G.; Olmsted, P. D., Simulation Studies of Stratum Corneum Lipid Mixtures.

Biophysical Journal 2009, 97 (7), 1941-1951. 11.

Höltje, M.; Förster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Höltje, H.-D., Molecular

Dynamics Simulations of Stratum Corneum Lipid Models: Fatty Acids and Cholesterol. Biochimica et Biophysica Acta (BBA) - Biomembranes 2001, 1511 (1), 156-167. 12.

Pandit, S. A.; Scott, H. L., Molecular-Dynamics Simulation of a Ceramide Bilayer. Journal of

Chemical Physics 2006, 124 (1), 7.

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

13.

Page 32 of 36

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. 14.

Wang, E.; Klauda, J. B., Molecular Dynamics Simulations of Ceramide and Ceramide-

Phosphatidylcholine Bilayers. The Journal of Physical Chemistry B 2017, 121 (43), 10091-10104. 15.

Moore, T. C.; Hartkamp, R.; Iacovella, C. R.; Bunge, A. L.; McCabe, C., Effect of Ceramide Tail

Length on the Structure of Model Stratum Corneum Lipid Bilayers. Biophys J 2018, 114 (1), 113-125. 16.

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. 17.

Ohman, H.; Vahlquist, A., In Vivo Studies Concerning a pH Gradient in Human Stratum

Corneum and Upper Epidermis. Acta Derm Venereol 1994, 74 (5), 375-379. 18.

Riddick, J. A.; Bunger, W. B.; Sakano, T.; Weissberger, A., Organic Solvents : Physical

Properties and Methods of Purification. 4th ed. ed.; Wiley: New York :, 1986. 19.

Kim, Y. C.; Ludovice, P. J.; Prausnitz, M. R., Transdermal Delivery Enhanced by Magainin Pore-

Forming Peptide. J Control Release 2007, 122 (3), 375-383. 20.

Ramos, A. P.; Lague, P.; Lamoureux, G.; Lafleur, M., Effect of Saturated Very Long-Chain Fatty

Acids on the Organization of Lipid Membranes: A Study Combining H-2 NMR Spectroscopy and Molecular Dynamics Simulations. Journal of Physical Chemistry B 2016, 120 (28), 6951-6960. 21.

Stahlberg, S.; Skolova, B.; Madhu, P. K.; Vogel, A.; Vavrova, K.; Huster, D., Probing the Role of

the Ceramide Acyl Chain Length and Sphingosine Unsaturation in Model Skin Barrier Lipid Mixtures by H-2 Solid-State NMR Spectroscopy. Langmuir 2015, 31 (17), 4906-4915. 22.

Groen, D.; Gooris, G. S.; Barlow, D. J.; Lawrence, M. J.; van Mechelen, J. B.; Demé, B.;

Bouwstra, J. A., Disposition of Ceramide in Model Lipid Membranes Determined by Neutron Diffraction. Biophysical Journal 2011, 100 (6), 1481-1489. 23.

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-1865. 24.

Jo, S.; Kim, T.; Im, W., Automated Builder and Database of Protein/Membrane Complexes for

Molecular Dynamics Simulations. Plos One 2007, 2 (9), 9. 25.

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.; Klauda, J. B.; Im, W., CHARMM-GUI Membrane Builder Toward Realistic Biological Membrane Simulations. Journal of Computational Chemistry 2014, 35 (27), 19972004.

32 ACS Paragon Plus Environment

Page 33 of 36 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

26.

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. Journal of Computational Chemistry 2005, 26 (16), 1781-1802. 27.

Klauda, J. B.; Venable, R. M.; Freites, J. A.; O'Connor, J. W.; Mondragon-Ramirez, C.;

Vorobyov, I.; Tobias, D. J.; 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), 7830-7843. 28.

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. 29.

Durell, S. R.; Brooks, B. R.; Bennaim, A., Solvent-Induced Forces between Two Hydrophilic

Groups. J. Phys. Chem. 1994, 98 (8), 2198-2202. 30.

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. 31.

Martyna, G. J.; Tobias, D. J.; Klein, M. L., Constant Pressure Molecular Dynamics Algorithms.

Journal of Chemical Physics 1994, 101 (5), 4177-4189. 32.

Steinbach, P. J.; Brooks, B. R., New Spherical-cutoff Methods for Long-range Forces in

Macromolecular Simulation. Journal of Computational Chemistry 1994, 15 (7), 667-683. 33.

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. 34.

Andersen, H. C., Rattle - a Velocity Version of the Shake Algorithm for Molecular-Dynamics

Calculations. Journal of Computational Physics 1983, 52 (1), 24-34. 35.

Humphrey, W.; Dalke, A.; Schulten, K., VMD: Visual Molecular Dynamics. Journal of

Molecular Graphics & Modelling 1996, 14 (1), 33-38. 36.

Racine, J., gnuplot 4.0: A Portable Interactive Plotting Utility. Journal of Applied Econometrics

2006, 21 (1), 133-141. 37.

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. Journal of Computational Chemistry 2009, 30 (10), 1545-1614. 38.

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. 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.

33 ACS Paragon Plus Environment

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

40.

Page 34 of 36

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.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E., 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.

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. 45.

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. 46.

Dahlen, B.; Pascher, I., Molecular Arrangements in Sphingolipids - Thermotropic Phase-

Behavior of Tetracosanoylphytosphingosine. Chemistry and Physics of Lipids 1979, 24 (2), 119-133. 47.

Boughter, C. T.; Monje-Galvan, V.; Im, W.; Klauda, J. B., Influence of Cholesterol on

Phospholipid Bilayer Structure and Dynamics. Journal of Physical Chemistry B 2016, 120 (45), 1176111772. 48.

Adams, M.; Wang, E.; Zhuang, X.; Klauda, J. B., Simulations of Simple Bovine and Homo

Sapiens Outer Cortex Ocular Lens Membrane Models with a Majority Concentration of Cholesterol. Biochimica et Biophysica Acta (BBA) - Biomembranes In press. https://doi.org/10.1016/j.bbamem.2017.11.010 49.

Gooris, G. S.; Bouwstra, J. A., Infrared Spectroscopic Study of Stratum Corneum Model

Membranes Prepared from Human Ceramides, Cholesterol, and Fatty Acids. Biophysical Journal 2007, 92 (8), 2785-2795. 50.

McIntosh, T. J., Organization of Skin Stratum Corneum Extracellular Lamellae: Diffraction

Evidence for Asymmetric Distribution of Cholesterol. Biophysical Journal 2003, 85 (3), 1675-1681. 51.

Notman, R.; den Otter, W. K.; Noro, M. G.; Briels, W. J.; Anwar, J., The Permeability Enhancing

Mechanism of DMSO in Ceramide Bilayers Simulated by Molecular Dynamics. Biophysical Journal 2007, 93 (6), 2056-2068. 34 ACS Paragon Plus Environment

Page 35 of 36 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

52.

Strader, M. L.; Feller, S. E., A Flexible All-Atom Model of Dimethyl Sulfoxide for Molecular

Dynamics Simulations. J. Phys. Chem. A 2002, 106 (6), 1074-1080. 53.

Mizushima, H.; Fukasawa, J.; Suzuki, T., Phase Behavior of Artificial Stratum Corneum Lipids

Containing a Synthetic Pseudo-Ceramide: a Study of the Function of Cholesterol. Journal of Lipid Research 1996, 37 (2), 361-367. 54.

Björklund, S.; Nowacka, A.; Bouwstra, J. A.; Sparr, E.; Topgaard, D., Characterization of

Stratum Corneum Molecular Dynamics by Natural-Abundance 13C Solid-State NMR. PLOS ONE 2013, 8 (4), e61889. 55.

Das, C.; Noro, M. G.; Olmsted, P. D., Lamellar and Inverse Micellar Structures of Skin Lipids:

Effect of Templating. Physical Review Letters 2013, 111 (14), 5.

35 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36

36 ACS Paragon Plus Environment