Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple

Article Cite This: J. Phys. Chem. B 2018, 122, 2757−2768

pubs.acs.org/JPCB

Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple Interior Stratum Corneum Models Eric Wang† and Jeffery B. Klauda*,†,‡ †

Department of Chemical and Biomolecular Engineering and ‡Biophysics Graduate Program, University of Maryland, College Park, Maryland 20742, United States

Downloaded via DURHAM UNIV on July 9, 2018 at 07:00:24 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

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 effects of lipid diversity and temperature are 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 this class of lipids and fatty acids for development of more complex SC models.

1. 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 matrix.1 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 dermatitis.2,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 microscopy,4 Fourier transform infrared spectroscopy,5 and most significantly X-ray diffraction.1,6,7 From these methods, the existence of two distinct lamellar phases has been determined. The short periodicity phase (SPP) has a repeat distance of ∼6.4 nm, while the long periodicity phase (LPP) has a repeat distance of ∼13.4 nm. While the molecular structure of these phases is not concretely known, it is believed that the SPP is characterized as a bilayer.7 The LPP is more complicated, and experimental models have proposed a multilayer structure with bilayers surrounding a monolayer center.5 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 © 2018 American Chemical Society

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 field.8−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 (C36) force field, which is known to better agree with experiment for sphingolipids,13 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 bonding.14 Therefore, using the C36 force field to simulate more complex systems such as Cer/Chol/FFA is highly relevant for studying Received: January 11, 2018 Revised: February 14, 2018 Published: February 21, 2018 2757

DOI: 10.1021/acs.jpcb.8b00348 J. Phys. Chem. B 2018, 122, 2757−2768

Article

The Journal of Physical Chemistry B

Table 1. System Parameters for Each Bilayer Systema

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 work,16 and lignoceric acid (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 to 5.3.17 The pKa of FFA is ∼5,18 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 to 6.9,17 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 FFA, 20 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 on bilayer properties and possibly 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 mixtures.21 In the past, neutron scattering results22 have been compared against MD,15 which is also done in this work. However, the experimental lipid mixtures are more than those used in 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), N-lignoceroylsphingosine (Cer24), Chol, and 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.

a

system

No. Cer

No. Chol

No. LA

T (°C)

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

Number of lipids and water are reported for the whole system.

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. On the basis of our previous work of ceramide bilayers,14 this water/lipid ratio likely produces fully hydrated systems and is useful for saving computational resources. Potassium counterions were placed in ternary mixtures to maintain an electroneutral system, while pure Cer systems did not contain any ions. Triplicate replicas of each system were generated in rectangular boxes and tetragonal cells (X = Y ≠ Z) using CHARMM-GUI Membrane Builder and equilibrated using the standard six-step process in NAMD software.23−26 Randomized initial configurations generated by CHARMM-GUI are justified by the agreement between randomized MD and RWMD.15 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. The CHARMM36 lipid force field27 with updates for sphingolipids13 was used with the TIP3P water model.28,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 Nosé−Hoover− Langevin piston,30,31 and constant temperature was maintained through Langevin dynamics. The Lennard-Jones potential with a force-switching function32 over the range from 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 1 × 10−6, with electrostatics being calculated every step. Bond lengths involving hydrogen were kept constant by the RATTLE algorithm.34 The simulations time step was 2 fs, and coordinates were saved every 5 ps, except in Cer16/Chol/LA systems at 65 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

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 2758


Article

The Journal of Physical Chemistry B modulus (KA), component SA/lip, deuterium order parameters (SCD), Chol tilt, gel phase tilt angle of pure Cer systems, electron density profiles (EDPs), bilayer thickness, interdigitation, hydrogen bonding (H bonding), lipid clustering, and twodimensional radial distribution functions (2D RDFs). Most analyses were performed with CHARMM,37,38 and averages with standard errors were calculated from triplicates. SCD was calculated by the formula SCD =

3 2 1 cos θ − 2 2

(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 Nlinked 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). With these coordinates, Quickhull39 constructs a Voronoi tessellation in which each atom is represented by a polygon that 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 KA =

kBT ⟨A⟩ Nσ⟨A⟩2

(2)

where kB is Boltzmann’s constant, T is the absolute temperature, ⟨A⟩ is the average SA of entire leaflet, N is the number of lipids per leaflet, and σ⟨A⟩2 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 from −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 halfmaximums of the water EDP, DHH is the distance between the peaks of the total EDP, and 2DC is the distance between halfmaximums of the hydrophobic EDP. For comparison with experimental neutron scattering length density (NSLD) data, SIMtoEXP40 was used to obtain a lipid NSLD profile. From the profile, the peak-to-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 +L

λov =

∫−L

4

ρt (z) × ρb (z) (ρt (z) + ρb (z))2

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 shown in cyan with colored headgroups, Chol in blue, LA in red, water as purple dots, and potassium ions as brown beads.

Two-dimensional (2D) RDFs were calculated with a Δr value of 0.1 Å for Cer−Cer, Cer−Chol, 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

dz (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 nonzero density (50 Å in this work). 2759


Article

The Journal of Physical Chemistry B group.14 Lipid clustering was analyzed using a python scikitlearn software41 with a density-based spatial clustering of applications with noise (DBSCAN) algorithm.42 Results from two cutoff distances are reported. A cutoff 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 cutoffs. 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.

of pure Chol have much higher KA than pure Cer,8 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 value of ∼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. On the one hand, the proximity of charged LA could inhibit tight packing and increase permeability. On the other hand, the presence of counterions provides a screening effect that can allow for charge proximity, and deprotonation of FFA in phospholipid bilayers has been found to increase chain order.20 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 ∼43

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)

Table 3. SA/lip (Å2) of Specific Lipid Types in All Systems

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)


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

system

T (°C)


32 32 50 65 80 32 32 50 65 80

Cer (Å2) 43.6 57.3 60.0 61.9 65.3 42.8 57.3 58.8 61.4 64.4

± ± ± ± ± ± ± ± ± ±

0.3 0.6 0.3 0.2 0.2 0.2 0.2 0.3 0.4 0.3

Chol (Å2)

LA (Å2)

± ± ± ±

15.9 ± 0.2 16.1 ± 0.1 16.46 ± 0.04 17.26 ± 0.04

25.4 25.4 27.1 28.3

0.3 0.2 0.1 0.2

24.6 ± 0.2 25.2 ± 0.1 26.0 ± 0.2 27.64 ± 0.03

16.05 ± 0.03 16.1 ± 0.1 16.7 ± 0.1 17.40 ± 0.04

Å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 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, 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. 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 agreement,14 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,

drastically 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 ∼12% from 32 to 80 °C. Similar linear trends have been observed for pure PC systems,43 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. The area compressibility modulus (KA), a measure of a membrane’s rigidity, is known to be increased in the presence of Chol44 and Cer,14 and Cer has also been referred to as a “Chol surrogate” because of this ability.45 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 Newton per meter 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 2760


Article


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 value of ∼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 ∼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 solidordered state, while the pure Cer systems are in the gel state. At 80 °C, the values peak at ∼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 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, because of 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 simulations,14,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 systems,44 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 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. Upon visual inspection of the pure Cer systems, it is apparent that, unlike the mixed systems, these have more bilayer tilt and

order parameters are sorted to be monotonically decreasing for comparison purposes. The C36 SCD is in overall good agreement with NMR, particularly for Cer16 and LA at 65 and 80 °C (Figure 2).

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.

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 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 nearly 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/ 2761


Article

The Journal of Physical Chemistry B Table 4. Average Tilt Angle for Chol in Mixed Systems and Cer Chainsa in Pure Cer system

T (°C)


32 32 50 65 80 32 32 50 65 80

Cer (deg)

Chol (deg)

17 ± 1 9.0 ± 0.1 10.2 ± 0.1 11.6 ± 0.1 13.5 ± 0.1 18.5 ± 0.6 8.9 ± 0.1 9.4 ± 0.1 10.85 ± 0.03 12.80 ± 0.02

Figure 3. Comparison of NSLD profiles obtained from C36 for Cer24/Chol/LA at 32 °C and experimental neutron scattering of the SPP.

a

Cer chain tilt is calculated as an average between fatty acid and sphingosine chains.

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 that varies in hydroxylation. These include N-lignoceroylsphingosine (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 has exclusively NS bases, so there is less potential for H bonding, which decreases chain order and bilayer thickness. 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

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 C36,14 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 work.14 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 systems,14 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 nonlamellar structure that Cer24 takes in crystallography. Because of 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 place in SC modeling as long as there is an understanding that 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 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-to-peak 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

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


Article

The Journal of Physical Chemistry B 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. Component EDPs are given in Figures S6 and S7. 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 seen in mixtures containing Chol.44 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

Table 5. Thicknesses (Å) of Each System for Different Regionsa system

T (°C)


32 32 50 65 80 32 32 50 65 80

DHH 39.4 44.7 44.5 43.9 43.0 50.5 49.1 48.7 48.1 46.8

± ± ± ± ± ± ± ± ± ±

0.3 0.1 0.2 0.2 0.0 0.4 0.1 0.1 0.1 0.0

DB

2DC

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

a

DHH is the head-to-head distance, DB is the overall bilayer thickness, and 2DC is the hydrophobic thickness.

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. The interdigitation of each lipid is reported in Table 6. The values are obtained using an integration of the density profile Table 6. Interdigitation (Å) of Individual Lipids for Each System system

T (°C)

Cer (Å)


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

Chol (Å)

LA (Å)

2.1 ± 0.1 2.23 ± 0.09 2.43 ± 0.02 3.2 ± 0.2

8.4 ± 0.1 8.4 ± 0.1 8.51 ± 0.03 8.81 ± 0.02

1.01 ± 0.08 1.17 ± 0.03 1.4 ± 0.04 1.74 ± 0.01

7.21 ± 0.04 6.7 ± 0.2 6.65 ± 0.07 6.71 ± 0.02

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 by three- to fourfold from Cer16 to Cer24 in mixture systems and by twofold 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 2763


Article


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, ∼10% from 32 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 to 80 °C decreases H bonding the most, which is likely due to the drop in KA. Although the drop itself is not very large, the KA at 80 °C is ∼0.6 N/m, which is approaching values for fluid bilayers (0.3 N/m) and makes differences more profound. Individual bond frequencies for the most significant bonds are reported in Tables S1 and S2. As is common in sphingolipids, the amide-carbonyl bond is the most significant bond with high frequencies in the 0.3−0.4 range for pure Cer. However, because of competition with bonding groups on Chol and LA, the frequency drops significantly in mixed systems to ∼0.15. While the overall Cer H bonding increases due to the addition of other bonds, the loss of the amide-carbonyl 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 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

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 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. 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 Table 7. Hydrogen Bonda Count for Overall Lipids in Each Systemb lipid

system

T (°C)

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

Chol

LA

NHB/lip 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.192 0.199 0.186

± ± ± ±

0.008 0.002 0.002 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

a

Defined as when the distance between donor and acceptor is less than 2.4 Å and the angle is greater than 150°. bValues are reported as the fraction of lipids with the considered hydrogen bond.

general, Cer has the most H-bonding as the Cer frequency ranges from 0.3 to 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 2764


Article


Table 8. Average Number of Clustered Lipids (NLT) and Number of Clusters (NCT) for Overall Clusteringa with a Cutoff Distance of 4.5 Å system

T (°C)

NLT

NLCER

NCT

NCCER


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

a

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 three lipids and greater were considered.

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

Figure 5. Visual representation of lipid clustering for top leaflet in (a) pure Cer16 and (b) Cer16/Chol/LA at 32 °C. (yellow ○) Cer lipids, (blue ○) Chol, and (red ○) 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.

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 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 Chol and usually have repeat distances in the 5.0−5.5 Å range.47,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 it is not drastically different. LA−LA RDFs (Figure S11) are weaker than both Cer−Cer and Chol−Chol RDFs, reaching peaks values of nearly 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 Figures S12−S14. 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, Cer−LA and Cer−Chol interactions never far exceed

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 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.5 to 4.0 in pure systems to ∼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. 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 2765


Article


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 ∼(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 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 time scales 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 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 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

the bulk and are likely less significant than Chol−LA interactions.

4. DISCUSSION AND CONCLUSION Because of 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 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 0.5 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 Bouwstra 49 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 distributed.50 Notman et al.51 have used the GROMOS force field to study the effect of dimethyl sulfoxide (DMSO) on SC bilayers. We attempted to replicate these effects in C36 but 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 parametrized with respect to its interaction with 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 that 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 2766


Article


Maryland Advanced Research Computing Center. The authors would also like to thank Dr. D. Huster in providing more detail on pH conditions of his experimental measures.

with nonlamellar 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 nonlamellar phases. Furthermore, complete exploration of the hightemperature 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 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. On the basis of 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 bilayers,14 supports the use of C36 for further simulations of mixed bilayers and even multilayers as models of the SC.

■

■

(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. Biochim. Biophys. Acta, Biomembr. 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. Arch. Dermatol. Res. 1990, 282 (8), 549−551. (3) Bouwstra, J. A.; Ponec, M. The Skin Barrier in Healthy and Diseased State. Biochim. Biophys. Acta, Biomembr. 2006, 1758 (12), 2080−2095. (4) Breathnach, A. S.; Goodman, T.; Stolinski, C.; Gross, M. FreezeFracture 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. J. Invest. Dermatol. 1989, 92 (2), 251−257. (6) White, S. H.; Mirejovsky, D.; King, G. I. Structure of Lamellar Lipid Domains and Corneocyte Envelopes of Murine StratumCorneum - An X-ray Diffraction Study. Biochemistry 1988, 27 (10), 3725−3732. (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. J. Lipid Res. 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. J. Phys. Chem. 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. J. Phys. Chem. B 2016, 120 (49), 12536−12546. (10) Das, C.; Noro, M. G.; Olmsted, P. D. Simulation Studies of Stratum Corneum Lipid Mixtures. Biophys. J. 2009, 97 (7), 1941− 1951. (11) Höltje, M.; Förster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Höltje, H.-D. Molecular Dynamics Simulations of Stratum Corneum Lipid Models: Fatty Acids and Cholesterol. Biochim. Biophys. Acta, Biomembr. 2001, 1511 (1), 156−167. (12) Pandit, S. A.; Scott, H. L. Molecular-Dynamics Simulation of a Ceramide Bilayer. J. Chem. Phys. 2006, 124 (1), 7. (13) 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. Biophys. J. 2014, 107 (1), 134− 145. (14) Wang, E.; Klauda, J. B. Molecular Dynamics Simulations of Ceramide and Ceramide-Phosphatidylcholine Bilayers. J. Phys. Chem. 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. J. Chem. Theory Comput. 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.; Wiley: New York, 1986.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b00348. 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 (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jeffery B. Klauda: 0000-0001-8725-1870 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research is supported by the NSF Grant MCB-1149187 and by a grant to the Univ. 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 Univ. of Maryland. Production runs used the Extreme Science and Engineering Discovery Environment allocation by Grant No. MCB-100139 on Stampede2 and the computational resources at the 2767


Article


molecular 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. (40) Kuč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) Ester, M.; Kriegel, H.-P.; Sander, J.; Xu, X. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. In Proceedings of 2nd International Conference on Knowledge Discovery and Data Mining; AAAI Press, 1996; Vol. 96, Issue 34, pp 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. Biochim. Biophys. Acta, Biomembr. 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. Biophys. J. 2007, 92 (3), 920−927. (46) Dahlen, B.; Pascher, I. Molecular Arrangements in Sphingolipids - Thermotropic Phase-Behavior of Tetracosanoylphytosphingosine. Chem. Phys. 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. J. Phys. Chem. B 2016, 120 (45), 11761−11772. (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. Biochim. Biophys. Acta, Biomembr. https://doi.org/10.1016/j.bbamem. 2017.11.010.201710.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. Biophys. J. 2007, 92 (8), 2785−2795. (50) McIntosh, T. J. Organization of Skin Stratum Corneum Extracellular Lamellae: Diffraction Evidence for Asymmetric Distribution of Cholesterol. Biophys. J. 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. Biophys. J. 2007, 93 (6), 2056−2068. (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 PseudoCeramide: a Study of the Function of Cholesterol. J. Lipid Res. 1996, 37 (2), 361−367. (54) Bjö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. Phys. Rev. Lett. 2013, 111 (14), 5.

(19) Kim, Y. C.; Ludovice, P. J.; Prausnitz, M. R. Transdermal Delivery Enhanced by Magainin Pore-Forming Peptide. J. Controlled 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. J. Phys. Chem. 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 H2 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.; Demé, B.; Bouwstra, J. A. Disposition of Ceramide in Model Lipid Membranes Determined by Neutron Diffraction. Biophys. J. 2011, 100 (6), 1481−1489. (23) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: a Webbased 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.; DavilaContreras, 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. J. Comput. Chem. 2014, 35 (27), 1997−2004. (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. J. Comput. Chem. 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. J. Phys. Chem. 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. J. Chem. Phys. 1995, 103 (11), 4613−4621. (31) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101 (5), 4177− 4189. (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) 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. J. Comput. Phys. 1983, 52 (1), 24−34. (35) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 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. J. Comput. Chem. 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 Macro2768


Simulations of Pure Ceramide and Ternary Lipid Mixtures as Simple

Recommend Documents