Molecular Dynamics Simulation Study of Skin Lipids - ACS Publications

Aug 14, 2015 - Tata Research Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar Industrial Estate, Pune - 411013, India. •S Sup...
0 downloads 4 Views 8MB Size
Article pubs.acs.org/JPCB

Molecular Dynamics Simulation Study of Skin Lipids: Effects of the Molar Ratio of Individual Components over a Wide Temperature Range Rakesh Gupta and Beena Rai* Tata Research Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar Industrial Estate, Pune - 411013, India S Supporting Information *

ABSTRACT: Atomistic molecular dynamics (MD) simulations were employed to systematically investigate the effects of the molar ratio of the individual components cholesterol (CHOL), free fatty acid (FFA), and ceramides (CER) on the properties of the skin lipid bilayer over a wide temperature range (300−400 K). Several independent simulations were performed for bilayers comprised of only CER, CHOL, or FFA molecules as well as those made up of a mixture of CER:CHOL:FFA molecules in different molar ratios. It was found that CHOL increases the stability of the bilayer, since the mixed (CER:CHOL:FFA) 1:1:0, 1:1:1, and 2:2:1 bilayers remained stable until 400 K while the pure ceramide bilayer disintegrated around ∼390 K. It was also observed that CHOL reduces the volume spanned by ceramide molecules, thereby leading to a higher area per CER and FFA molecule in the mixed bilayer system. The CHOL molecule provided more rigidity to the mixed bilayer and led to a more ordered phase at elevated temperatures. The CHOL molecule provided fluidity to the bilayer below the phase transition temperature of CER and kept the bilayer rigid above the phase transition temperature. The FFA interdigitizes with CER molecules and increases the thickness of the bilayer, while rigid CHOL decreases the bilayer thickness. The presence of CHOL increases the compressibility of the bilayer which is responsible for the high barrier function of skin. The CER molecule forms inter- and intramolecular hydrogen bonds, while CHOL only forms intermolecular hydrogen bonds.

1. INTRODUCTION The accurate prediction of dermal uptake of chemicals is relevant to both transdermal drug delivery as well as topical application of cosmetics. Extensive research has been performed over the last several decades to predict skin permeability of various molecules.1−4 These efforts include the development of empirical approaches such as quantitative structure−permeability relationships5 and porous pathway theories6 as well as the establishment of rigorous structurebased models.5 However, a molecular-level understanding of the skin’s surface layerthe Stratum Corneum (SC)which shall ultimately lead to the development of rapid in-silico screens to predict permeability of a molecule, from knowledge of its molecular structure alone, is still not on the horizon. SC, the outermost part of the skin epidermis, provides the main barrier to penetration of pathogens and against water loss. The SC mainly is comprised of corenocytes (brick) and lipid matrix (mortar).7,8 The extracellular lipid matrix is the key determinant for its barrier functions and is composed of a heterogeneous mixture of long chain ceramides (CER), cholesterol (CHOL), and free fatty acid (FFA) in certain ratios.9,10 The heterogeneous nature of SC arises due to the presence of more than 300 different CER species comprised of varying head groups and fatty chains.11 Each CER is made up of a sphingoside motif which is bound to a fatty acid chain. CERs can be classified into sphingosines (S), phytosphingosines (P), © 2015 American Chemical Society

and hydroxysphingosines (H) according to functional group present in the sphingosine motif. CERs can be further classified as nonhydroxy (N), α-hydroxy (A), and esterified ω-hydroxy fatty acid CERs according to the hydroxyl group attached to the fatty acid chain.12 Detailed structures of different types of ceramides present in the human skin can be found elsewhere.6,13 There are some approaches proposed for breaching the skin barrier such as chemical penetration enhancers, electroporation, sonophoresis, and microneedles, but a complete understanding of these molecular processes is still in its infancy due to the limitation associated with experimental techniques.14 Molecular simulation offers a way to yield important physical insight at molecular level resolution. However, the investigation of a realistic skin behavior and structure at the molecular level is challenging because of the heterogeneity of ceramides and free fatty acids present in SC.15,16 In human SC ceramides, the length of the sphingo motif is always found to be between 16 and 18 carbons but the length of fatty acid varies from 16 to 34 carbons, thereby leading to asymmetry in the length of ceramide chains.17,18 Atomistic modeling of such a complex system incorporating all Received: March 3, 2015 Revised: August 2, 2015 Published: August 14, 2015 11643

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

Hoover thermostat with the same time constant. Pressure was maintained at 1 bar using a Berendsen barostat for initial equilibration with a time constant of 0.5 ps and a compressibility of 4.5 × 10−5 bar−1; later, in a production run, it was changed to a Parrinello−Rahman barostat with a time constant of 5 ps and a compressibility of 4.5 × 10−5 bar−1. The semi-isotropic coupling (coupled separately in the xy and z direction) of the barostat was used to simulate a tensionless bilayer. All of the bonds in lipid molecules were constrained using the LINCS algorithm,36 while for water the SETTLE37 algorithm was used. A time step of 2 fs was used for simulations below 360 K and 1 fs for higher temperature. It has been shown earlier that the cutoff for electrostatic interaction leads to some artifacts in phospholipid bilayer simulation.38 A cutoff of 1.2 nm was used for van der Waals and electrostatic interactions. Long range electrostatic interaction was computed using the particle mesh Ewald (PME) method. The equilibration time for all of the systems was 10 ns (NVT) followed by 100 ns (NPT) equilibration. The final 10 ns run of NPT simulation was used for calculation of equilibrium properties (production run). The configuration was sampled at every 0.5 ps in the production run. All of the properties were averaged over five blocks of 2 ns. 2.3. Initial Bilayer Structures. The molecule structures of the individual components used in these simulations are shown in Figure 1. We have used the hairpin configuration of CER

possible ceramides and fatty acids found in SC is far away from current computational capabilities. However, in order to simulate a realistic SC layer,19 we have chosen the most abundant ceramide, CER-NS 24:0, and free fatty acid, FFA 24:0. Holtze et al.20 for the first time reported the simulation of two-component mixtures of fatty acid and cholesterol, Pandit and Scott21 studied the single component CER-NS (16:0) bilayer, and Notman et al.22 investigated the effect of DMSO on the pure CER-NS (24:0) bilayer. The first study of threecomponent mixtures of CER NS (24:0), free fatty acid (24:0), and cholesterol was performed by Das et al.23 The authors investigated the effect of the molar ratio of an individual component on structural properties and followed by diffusion and permeability calculation for CER and the mixed bilayer.24 The same authors have also reported on the modeling of the corneocyte wall and proposed that lamellar layering is induced by the patterned corneocyte wall.13 Recently, Guo et al.25 performed both united atom and all atom molecular dynamics simulation of single component CER-NS (16:0) and CER-NP (16:0) using the Berger26 and CHARMM27 force fields. They found that the OH group in CER-NP aligned itself perpendicular to the bilayer normal as compared with CERNS, resulting in a more hydrogen bonded network, which leads to a high gel-to-liquid phase transition temperature in CER-NP. Guo et al. performed simulations for a symmetric CER-NS bilayer which lead to a proper ordering of lipid tails in its interior, resulting in a high gel to liquid phase temperature. To the best of the authors’ knowledge, no studies have been reported on the phase transition of mixed skin lipid bilayers. We report on the phase transition of skin lipids at various compositions along with the calculated area per CER molecule in multicomponent bilayers. The effects of CHOL and FFA in ceramide lipid bilayers over a wide temperature range is systematically investigated via atomistic molecular dynamics (MD) simulation in constant NVT and NPT conditions. In bilayer simulation, the area per lipid is a primary parameter which is directly correlated to the molecular organization of the lipid bilayer.28 It can be compared with experiments, and it gives information about the equilibration of the system in simulation. We have modified the method of Hofsäß et al.29 to calculate the area per lipid in the mixed bilayer.

2. SIMULATION DETAILS 2.1. Force Field. The potential function and interaction parameters for CER were taken from the Berger force field,26 which is properly parametrized for palmitoylsphingomyelin (PSM) bilayers and phospholipids. The PSMs have a similar structure to CER-NS except for the head group geometry. The head group size and partial charges in PSM are very high as compared to CER.30 The polar hydrogens were included explicitly, and charges on polar groups were taken from earlier simulations reported in the literature.22−24 The Ryckaert− Bellemans dihedral potential was used for lipid hydrocarbon tails.31 The topology for free fatty acid and cholesterol was the same as that used by Holtze et al.20 The simple point charge (SPC) model was used for the water molecule.32 The united atom carbon in all components had zero partial charge. 2.2. Simulation Setup. All simulations were carried out in the NVT and NPT ensemble using the latest GROMACS molecular dynamics package.33−35 The temperature was controlled by a Berendsen thermostat with a time constant of 0.5 ps for initial equilibration; later, it was changed to a Nosé−

Figure 1. Molecular structure of individual skin lipid molecules used in simulations. From left: Ceramide-NS (CER), cholesterol (CHOL), and free fatty acid (FFA). Chain sn1, chain sn2 and Ffa chain sn1 represent lipid tails of CER and FFA, respectively. The colors red, blue, white, and cyan represent oxygen O, nitrogen N, hydrogen H, and carbon C (united atom), respectively.

molecules.23 To remove bad contacts between molecules, energy minimization was carried out in a vacuum using the steepest descent and conjugate gradient method. For a singlecomponent bilayer system, an individual molecule was first minimized and then replicated using “genconf” (Gromacs) in the X and Y direction to obtain a single layer. This single layer was then replicated in the Z direction to obtain the initial bilayer configuration followed by energy minimization. This minimized bilayer was placed in a bigger box, and the box was filled with SPC water molecules. Further, the system was energy minimized and subjected to a 10 ns NVT MD run at 300 K, by keeping lipid molecules fixed for proper hydration. 11644

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

where Nffa and Nchol are the number of FFA and cholesterol molecules in the mixed bilayer, respectively. Vchol and Vffa are the volume per lipid for pure cholesterol and FFA bilayer, respectively, which were calculated by separate individual simulations of pure cholesterol and FFA system in the NPT ensemble at each temperature using identical simulation conditions such as time steps and cutoffs. The area per lipid (APL) is a primary parameter which describes the molecular packing of a lipid bilayer. It can be compared with experiments, and it gives information about the equilibration of the system in simulation. In a molecular dynamics simulation of a pure lipid bilayer which has the normal along the z direction, the APL can be calculated using the following equation

For a mixed bilayer, initial random configurations were generated using the tool PACKMOL.39,40 Initial and final configurations of mixed bilayers mostly correlated in gel phase simulations.20 We further performed simulated annealing to remove these artifacts. All the minimized structures were heated until 360 K and cooled back again to 300 K in a 1 ns run. Structures obtained from the simulated annealing were equilibrated for 100 ns at 300 K. These equilibrated structures were further used as a starting structure for a given temperature and ran for another 80−100 ns at 1 bar. The final 10 ns run of these simulations was used for analysis and computation of properties. The same procedure was followed for the simulations of pure CHOL and FFA bilayers. Table 1 shows Table 1. Number of Individual Molecules Used in Simulation for Each Molar Ratio

A=

systema

CER

CHOL

FFA

water

1:0:0 0:1:0 0:0:1 1:1:0 1:0:1 1:1:1 2:2:1

128 0 0 64 84 52 56

0 128 0 64 0 50 56

0 0 512 0 84 52 32

5120 5120 9020 5120 5120 5120 5120

The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratios of CER:CHOL:FFA in the bilayer.

the molar ratio used in the simulations and the corresponding number of individual molecules. The ratios in the Article should be considered in the order of CER:CHOL:FFA unless it is specified otherwise specifically.

3. AREA AND VOLUME PER MOLECULE In an MD simulation of a single-component bilayer which has the normal along the z direction, the volume per lipid (VPL) could be calculated by subtracting the volume of the water molecule from the total volume V = LxLyLz (1) V − NW VW Nlipid

havg =

V − NW VW − NcholVchol − NffaVffa Ncer

(4)

V − NW VW A

(5)

where V and A are the volume and area given by eqs 1 and 4, respectively. Therefore, the area per ceramide (APC) molecule can be written as

(2)

Acer =

where Vlipid is the volume per lipid, V is the total volume of the box calculated by simulation, Nw is the number of water molecules, and Vw is the volume span by a single water molecule under similar simulation conditions such as temperature, van der Waals, and Coulombic cutoffs. Hofsäß et al.29 calculated the volume per dipalmitoylphosphatidylcholine (DPPC) molecule in a mixed bilayer of DPPC and CHOL by assuming a constant volume of CHOL at all concentrations. We have used a similar procedure but modified the above equations for a mixed component system by assuming that the volumes of cholesterol and FFA molecules do not change with their number (molar concentration) although they change with temperature. We also assumed that all the volume has been occupied by molecules, i.e., there is no free volume in the bilayer. Considering the above assumption, the volume per CER molecule in a mixed bilayer system can be written as Vcer =

Nlipid

where Lx and Ly are the box lengths in the X and Y directions, respectively, and Nlipid is the total number of lipids in the bilayer. Some methods have been proposed for calculating the area per phospholipid molecule and area per CHOL in a mixed bilayer system at different temperatures and CHOL mole fractions.41,42 Chiu et al.42 assumed the area per lipid molecule to be independent of the mole fraction of CHOL and calculated the area per lipid by simple linear regression. Hofsäß et al.29 used the volume per lipid and common bilayer thickness to compute the individual area per molecule at different mole fractions of CHOL. We have employed a similar approach to calculate the area per molecule in our mixed bilayer system. We have assumed that all of the available volume was occupied by lipid molecules, which means there is no free volume inside the bilayer. The average thickness of a bilayer is given by

a

Vlipid =

2LxLy

Vcer havg

(6)

using eqs 3, 5, and 6 Acer =

⎛ V − NW VW − NcholVchol − NffaVffa ⎞ A ⎟ ⎜ V − NW VW ⎝ Ncer ⎠

Acer =

N V + NffaVffa ⎞ A ⎛ ⎟ ⎜1 − chol chol Ncer ⎝ V − NW VW ⎠

Achol = A ffa =

AVchol V − NW VW AVffa V − NW VW

(7)

(8)

(9)

where Achol and Affa are the area per molecule of pure cholesterol and FFA bilayer, respectively. Ncer, Nchol, and Nffa are the number of ceramide, cholesterol, and FFA molecules in

(3) 11645

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B the mixed bilayer. Vchol and Vffa are the volume per molecule calculated from the pure CHOL and FFA bilayer simulation.

4. RESULTS AND DISCUSSION 4.1. Heating Run. Initially equilibrated bilayer structures at 300 K were used for the heating run. All the bilayers were

Figure 4. Snapshots of individual components and lateral packing of lipid tails in the upper leaflet (at a plane ∼1.4 nm above the bilayer center) of a mixed (CER:CHOL:FFA) 1:1:1 bilayer at different temperatures. From the top, the snapshots are at 300, 340, 360, 370, 380, and 400 K, respectively. CER, CHOL, and FFA are shown in light green, pink, and yellow color, respectively.

Figure 2. Evolution of temperature, area per lipid, and box volume with simulation time in heating run. The colors black, red, green, blue, yellow, cyan, and magenta represent the systems at compositions (CER:CHOL:FFA) of 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1, and 2:2:1, respectively.

Figure 5. Snapshots of individual components and lateral packing of lipid tails in the upper leaflet (at a plane ∼1.4 nm above the bilayer center) of a mixed (CER:CHOL:FFA) 2:2:1 bilayer at different temperatures. From the top, the snapshots are at 300, 340, 360, 370, 400, 410, and 420 K, respectively. CER, CHOL, and FFA are shown in light green, pink, and yellow color, respectively.

heated from 300 to 450 K slowly at a rate of 5 K/ns with a time step of 1 fs. Figure 2 shows the evolution of temperature, area per lipid, and box volume with simulation time. The constant slope of the temperature curve confirmed the uniform heating. Due to the rigid structure of the CHOL, the pure CHOL bilayer was stable for the whole heating run and pure FFA

Figure 3. Snapshots of the pure ceramide bilayer at different temperatures (top); arrangement of ceramide tails at a plane ∼1.2 nm above the bilayer center (below). The ceramide head group, tails, and water are shown in blue, cyan, and red, respectively.

11646

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B Table 2. Volumea and Areab per FFA/CHOL Molecule at Different Temperatures Calculated from Pure FFA (0:0:1) and CHOL (0:1:0) Bilayer Simulations, Respectively system

Nlipids

Nwater

temp (K)

Vwc (nm3)

Vchol, Vffa (nm3)

Affa, Achol (nm2)

FFA FFA FFA FFA CHOL CHOL CHOL CHOL CHOL CHOL

512 512 512 512 128 128 128 128 128 128

9020 9020 9020 9020 5120 5120 5120 5120 5120 5120

300 340 360 370 300 340 360 370 380 400

0.0302 0.0313 0.0319 0.0323 0.0302 0.0313 0.0319 0.0323 0.0327 0.0336

0.604 0.634 0.643 0.663 0.626 0.631 0.632 0.636 0.638 0.640

0.206 0.217 0.221 0.253 0.366 0.370 0.372 0.375 0.377 0.380

a

Calculated using eq 2. bCalculated using eq 4. cThe water volume Vw was calculated from the separate simulations of 2048 water molecules at the same simulation conditions used for the respective bilayer system.

disintegrated near 380 K. The bilayers disintegrated in the order of (CER:CHOL:FFA) 0:0:1, 1:0:1, 1:0:0, 1:1:1, 2:2:1,1:1:0, and 0:1:0. 4.2. Structures. Figure 3 shows snapshots of the pure CER bilayer and tail arrangement at each temperature. At physiological temperatures, the bilayer was found to be in the gel phase with proper ordering of lipid tails in both of the leaflets. The lipid tails were found to be in the hexagonal phase which is in agreement with earlier experimental studies.43−45 At higher temperature, this hexagonal phase disappears and we observe a phase change around 360−370 K. Earlier simulations20,21 and experimental studies43,46 have shown that CER-NS exhibits a gel to liquid crystalline phase transformation over a temperature range of 353−363 K. It is important to note that in the center of the bilayer an almost liquidlike phase was observed due to the asymmetric chain length of the CER tails (Figure 3).

Figure 7. Tail order parameter Sz of CER sn1 and sn2 chains at 300 K and above the phase transition temperature of 370 K. The chains sn1 and sn2 are shown in Figure 1. The 1:0:0, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratios of CER:CHOL:FFA in the bilayer.

Figures 4 and 5 show the individual component structure of the 1:1:1 and 2:2:1 bilayers at each temperature. At lower temperature ( ∼365 K), a

reverse phenomenon is observed. The order parameter is a little bit more in the 2:2:1 system as compared to the 1:1:1 system. The order parameter for 1:0:1 is the highest for both of the chains at each temperature. The same effect of CHOL is also seen in fatty acid chains, as shown in Figure 8. It follows similar trends as CER chains: below the melting point of FFA (∼357 11649

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

Figure 12. Electron density of each bilayer system along the bilayer normal z at different temperatures. The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratio of CER:CHOL:FFA in the bilayer.

K) as compared to the gel phase (300 K). The head group density profile gets broader with temperature because of more fluidity in the liquid crystalline phase. Figure 10 shows the density profiles of each component of the mixed 1:1:1 and 2:2:1 bilayers. The density profile is almost similar at all of the temperatures studied. It is worth noting that the peaks of CER and FFA density are near the interface, indicating that head groups of both components sit just below the lipid−water interface. The CHOL mostly sits in the highly dense region of the bilayer. Due to the smaller size and low partial charge of head groups in the ceramide bilayer as compared to the phospholipid, CHOL molecules sit much below the head group region, almost between the lipid tails. The CER and FFA profile shows a small peak at the center of the lipid bilayer at low temperature which was absent in the pure ceramide bilayer. This confirms that similar chain lengths of CER and FFA lead to proper interdigitization. These peaks get distorted with an increase in temperature because higher thermal energy leads to disordering in the tails, enabling a fluidlike environment. The lipid tail density decreases with an increase in temperature because the bilayer moves from the gel to liquid phase at higher temperature. Parts a and b of Figure 11 show the density profiles of individual components along the bilayer normal in each bilayer at 300 and 360 K, respectively. Overlapping of FFA and ceramide density peaks and interdigitization between them can easily be seen. FFA induces more and more ordering in lipid tails and increases the thickness of the bilayer, but CHOL reduces the ordering as well as the thickness of bilayer.16 The lipid density profile shows similar behavior at all temperatures studied. The peak to peak distance between lipid density was

K), the order parameter is higher in the 1:1:1 system, while, above this temperature, a reverse phenomenon is observed. Another interesting point to note is that the phase transition temperature is found to be lowest in the pure 0:0:1 bilayer and highest in 0:1:0 bilayer. This can be explained by CHOL’s high melting temperature and its rigid structure. 4.6. Lipid Density. In Figure 9, the densities of lipid components are plotted along the bilayer normal (z-axis) for the ceramide bilayer at different temperatures. The profile is similar to the four-region model as proposed by Marrink et al.59,60 The density from z ∼ −6.0 (not shown in the figure) to z ∼ −3.2 is constant, and this flat region is called bulk water; as we move further near the interface, the density of water decreases and the head group density increases (z ∼ −3.2 to z ∼ −2.0). The next region from z ∼ −2.0 to z ∼ −0.8 has very ordered and tightly packed lipid tails, which is mainly responsible for the barrier properties of skin. In the last region, the total lipid density has the minimum value because of loose and random packing of lipid tails. There is a sharp peak in ceramide density near the water− lipid interface, and all the head groups sit just below the water− lipid interface. The interfacial width is found to be smaller than those of the phospholipid membrane because of the small size and the low partial charge of the head group. Next to this region, the lipid tail density is very high because of proper ordering of tails, and the density in this region is found to be ∼1.3 g/cm3. In the middle of the bilayer, there is a dip in lipid density because of the asymmetric chain length; this region has more and more free volume and has a density of ∼0.7 g/cm3. The temperature induces disorder in the system which leads to a small tail density in the liquid crystalline phase (∼360−370 11650

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

Table 3. Distances between Head Group Atomsa in between the Two Leaflets of the Bilayer and Average Bilayer Thickness of Each Bilayer at Different Temperatures system

temp (K)

CER O21b (nm)

1:0:0 1:0:0 1:0:0 1:0:0 0:1:0 0:1:0 0:1:0 0:1:0 0:1:0 0:1:0 0:0:1 0:0:1 0:0:1 0:0:1 1:1:0 1:1:0 1:1:0 1:1:0 1:1:0 1:1:0 1:0:1 1:0:1 1:0:1 1:0:1 1:1:1 1:1:1 1:1:1 1:1:1 ̀ 1:1:1 1:1:1 2:2:1 2:2:1 2:2:1 2:2:1 2:2:1 2:2:1

300 340 360 370 300 340 360 370 380 400 300 340 360 370 300 340 360 370 380 400 300 340 360 370 300 340 360 370 380 400 300 340 360 370 380 400

5.182 5.136 5.135 5.123

CHOL O6b (nm)

FFA O1b (nm)

3.378 3.367 3.355 3.354 3.346 3.326 5.635 5.54 5.528 5.497 4.413 4.372 4.305 4.263 4.240 4.101 5.275 5.176 5.171 5.156 4.681 4.592 4.541 4.491 4.355 4.179 4.542 4.470 4.425 4.459 4.425 4.325

4.1 4.089 4.083 4.001 3.911 3.815 5.827 5.711 5.747 5.711 5.005 4.962 4.913 4.925 4.662 4.578 4.831 4.838 4.823 4.729 4.704 4.645

4.334 4.27 4.257 4.089 3.96 3.86 4.186 4.166 4.125 4.068 4.021 3.922

havg1c (method 1) (nm)

havg2c (method 2) (nm)

5.65 5.61 5.59 5.57 3.423 3.411 3.401 3.396 3.389 3.372 5.866 5.847 5.826 5.240 4.693 4.683 4.644 4.645 4.631 4.662 5.926 5.829 5.888 5.856 5.119 5.071 5.033 4.962 5.023 4.889 4.957 4.911 4.929 4.859 4.829 4.758

5.55 5.53 5.52 5.44 3.36 3.342 3.324 3.321 3.319 3.318 5.796 5.787 5.726 5.612 4.583 4.573 4.554 4.545 4.511 4.492 5.896 5.879 5.858 5.836 5.095 4.98 4.94 4.90 4.85 4.85 4.83 4.80 4.78 4.77 4.76 4.66

The error in the data is around ±0.05% (T < 350 K) to ±0.15% (T > 360 K). bCERO21, CHOLO6, and FFAO1 represent the atom name O21, O6, and O1 in ceramide, cholesterol, and fatty acid molecules, respectively. Atoms are shown in Figure 1. chavg1 and havg2 are calculated from the volume to area per lipid ratio and the drop in water density profile, respectively. Details of the methods are given in the Supporting Information (Figure S5). a

membrane,61 or the distance between peaks of the electron density profile.62 As shown in Figures 10 and 11, since CHOL sits in the interior space of the bilayer, it is not convenient to calculate the bilayer thickness by calculating the distance between some of the head group atoms. We have used two different methods to calculate the bilayer thickness. The schematic of the first approach is shown in the Supporting Information. The second approach is based on the volume per CER and the area per CER using eq 6. The results for the bilayer thickness of each system estimated by the above methods are summarized in Table 3. The bilayer thickness was found to be highest in 0:0:1 and lowest in 0:1:0 layers at all temperatures studied. The bilayer thickness in mixed bilayers was found to be lowest in the 1:1:0 bilayer and maximum in the 1:0:1 bilayer. This shows that the same chain length fatty acid interdigitizes with the ceramide molecule chain and increases the chain length on other side being small in size CHOL reduces the bilayer thickness. This phenomenon has

found to be highest in the ceramide bilayer and lowest in the mixed 2:2:1 bilayer. Also, overall density reduction can be seen in FFA and CHOL with increasing temperature, as illustrated earlier. Figure 12 shows the electric density of lipid bilayers along the bilayer normal in each bilayer system at different temperatures. The shape of the electron density in an individual bilayer is almost similar for all of the temperatures, but the density gets lower with an increase in temperature. It is interesting to note the mixed bilayers which comprise CHOL have two peaks, while a single peak was observed in the CHOL free mixed bilayer. This confirms that CHOL mostly sits in between the CER and FFA chains. 4.7. Interfacial and Bilayer Width. The interface (bilayer/ water) width was estimated as the distance over which the water density drops to 1/e of its bulk value.22,23 The bilayer thickness could be defined in several ways such as the distance between phosphate atoms in each leaflet in the phospholipid 11651

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

Figure 13. Distance between the head group atoms of CER (N23), CHOL (O6), and FFA (O1) in two leaflets. N23, O6, and O1 are shown in Figure 1. The 1:0:0, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratios of CER:CHOL:FFA in the bilayer.

4.8. Area Compressibility. The area compressibility of a bilayer whose normal is oriented along the Z axis is calculated as KA = k bT

⟨A⟩ ⟨A ⟩ − ⟨A⟩2 2

(12)

where A = Lx × Ly is the area of the bilayer in the XY plane, kB is the Boltzmann constant, and T is the temperature. The angular brackets denote the ensemble averages taken over the course of the simulation. Figure 14 shows the area compressibility of each bilayer at different temperatures. The compressibility was found to be highest in the pure CHOL bilayer, which supports the fact of rigidity of cholesterol and less permeability of the skin. The pure FFA bilayers are softer than pure CER and CHOL bilayers. The presence of FFA reduces the area compressibility and makes the bilayer soft in nature. The compressibility increases with an increase in the CHOL concentration for a fixed temperature. Additional bigger bilayer patches were simulated to check the effect of system size over the area compressibility. The details are given in Table S1 (Supporting Information). KA decreases with the increase in system size in a tensionless bilayer system. This has also been observed in earlier studies of a phospholipid bilayer.64 4.9. Hydrogen Bond. The numbers of hydrogen bonds present in bilayer systems were calculated on the basis of the geometrical criteria. A hydrogen bond exists if the distance between the donor and acceptor atoms is ≤0.35 nm and the acceptor−donor−hydrogen angle is ≤60°. Table 4 shows the number of hydrogen bonds between individual components of the bilayer at different temperatures. An earlier study showed that the number of hydrogen bonds per water molecule in liquid water is around 3.5−3.8 at 300 K.65 In all lipid bilayer systems, this number decreases to between 1.5 and 1.7 at 300 K, since at the interface CER forms the hydrogen with water. Having a small head group and a small partial charge, CHOL is unable to form hydrogen bonds with other CHOL molecules, although it does form hydrogen bonds with water molecules at the interface. CER molecules form more intramolecular hydrogen bonds as compared to CHOL and FFA. The number of hydrogen bonds decreases for each bilayer system with an increase in temperature because of

Figure 14. Dependence of area compressibility kA over the temperature range for each bilayer system. The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratios of CER:CHOL:FFA in the bilayer.

also been observed experimentally by Schroeter et al.63 The bilayer thickness decreases with increasing temperature because the percentage in area expansion is larger than that of the volume. The interfacial width (∼0.5 nm) is smaller than that of sphingomyelin bilayers (∼0.9 nm) because of the small head group of the ceramide layer.30 Figure 13 shows the distance between the head group atoms of CER (N23), CHOL (O6), and FFA (O1) in two leaflets. Inclusion of fatty acid increases the distance of N23 between the two leaflets, while CHOL decreases this distance. The minimum distance was found in the 1:1:0 bilayer, while the maximum distance was found in the 1:0:1 system. Another interesting thing we noticed was that, above the phase transition temperature of CER, in mixed bilayers, distances between these atoms in two leaflets were more in high CHOL containing bilayers, which supports the fact that above the phase transition temperature CHOL provides more ordering to the bilayer. 11652

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

Table 4. Number of Hydrogen Bondsa Formed in Different Components of Each Bilayer System at Different Temperatures systemb

temp (K)

CER-CER

CER-CHOL

CER-FFA

CER-Wc

CHOL-CHOL

CHOL-FFA

CHOL-W

FFA-FFA

FFA-W

W-W

1:0:0

300 340 360 370 380 300 340 360 370 380 400 300 340 360 370 300 340 360 370 380 400 300 340 360 370 300 340 360 370 380 400 300 340 360 370 380 400

1.52 1.42 1.36 1.34 1.27 0 0 0 0 0 0 0 0 0 0 0.78 0.75 0.72 0.68 0.66 0.64 1.24 1.19 1.14 1.08 0.76 0.72 0.71 0.73 0.66 0.65 0.80 0.75 0.71 0.69 0.65 0.63

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.293 0.280 0.272 0.265 0.256 0.247 0 0 0 0 0.31 0.28 0.25 0.23 0.23 0.22 0.34 0.33 0.33 0.29 0.26 0.25

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.59 0.57 0.54 0.50 0.19 0.15 0.14 0.12 0.14 0.13 0.13 0.11 0.09 0.08 0.08 0.07

2.13 2.01 1.95 1.90 1.88 0 0 0 0 0 0 0 0 0 0 3.35 3.00 2.94 2.82 2.72 2.51 2.0 1.84 1.75 1.64 2.93 2.59 2.56 2.49 2.43 2.38 2.72 2.68 2.63 2.61 2.59 2.47

∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0 ∼0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.12 0.10 0.10 0.10 0.08 0.07 0 0 0 0 0 0

0 0 0 0 0 2.03 1.92 1.86 1.82 1.78 1.71 0 0 0 0 0.74 0.64 0.56 0.50 0.48 0.43 0 0 0 0 0.88 0.79 0.77 0.76 0.72 0.69 0.86 0.83 0.76 0.74 0.72 0.63

0 0 0 0 0 0 0 0 0 0 0 0.28 0.27 0.24 0.20 0 0 0 0 0 0 0 0 0 0 1.06 1.44 1.61 1.54 1.68 1.62 0.85 0.71 0.61 0.62 0.65 0.59

0 0 0 0 0 0 0 0 0 0 0 1.66 1.56 1.49 1.45 0 0 0 0 0 0 0.26 0.24 0.23 0.20 1.78 1.76 1.71 1.69 1.63 1.61 1.81 1.75 1.68 1.67 1.66 1.57

1.70 1.60 1.55 1.51 1.50 1.70 1.60 1.55 1.52 1.49 1.44 1.67 1.57 1.52 1.49 1.69 1.59 1.54 1.51 1.49 1.43 1.71 1.61 1.57 1.54 1.69 1.59 1.54 1.51 1.48 1.46 1.67 1.59 1.54 1.52 1.48 1.44

0:1:0

0:0:1

1:1:0

1:0:1

1:1:1

2:2:1

The error in the data is around ±1% (T < 350 K) to ±3% (T > 360 K). bThe 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1, and 2:2:1 represent the molar ratio of CER:CHOL:FFA in the bilayer. cW stands for the water molecule.

a

high kinetic energy (thermal fluctuation) at elevated temperature.



ature. The fatty acids of a similar chain length of CER increase the thickness of the bilayer, while the short and rigid CHOL molecule decreases the thickness of the bilayer. The presence of CHOL increases the area compressibility of the bilayer which is responsible for the barrier function of skin, while FFA reduces the compressibility. The CER molecule forms intra- and intermolecule hydrogen bonds, while CHOL only forms intermolecule hydrogen bonds. In reality, human skin lipids are comprised of a heterogeneous mixture of polydisperse FFA (in terms of chain length) and ceramide acid motif (in terms of chain length and functional head group). This work is only related to the most abundant ceramide (CER-NS) and fatty acid found in the human skin. We are currently working on incorporating the effect of polydispersity in our simulations, and results will be presented in a later communication.

CONCLUSION

The effect of temperature on skin lipid bilayers has been studied at the molecular scale. We have performed molecular dynamics simulations for different composition ratios of the main constituents (CER-NS 24, CHOL, and FFA 24:0) of skin’s stratum corenum lipid matrix. We have found that the CER bilayer exhibits a hexagonal gel phase at normal skin temperature which gets converted into a liquid crystalline phase at ∼365 K. The asymmetric chain length of CER tails induces a fluidlike environment in the center of bilayers. The presence of symmetric chain length FFA induces more and more interdigitization and increases the ordering of CER tails. The CHOL molecule sits in the interior of lipid tails and provides more rigidity to the bilayer at higher temperature. CHOL molecules show a dual nature; they decrease the chain ordering mixed bilayer below the phase transition temperature of CER and increase the ordering above the phase transition temper11653

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B



(11) Notman, R.; Anwar, J. Breaching the Skin Barrier–Insights from Molecular Simulation of Model Membranes. Adv. Drug Delivery Rev. 2013, 65 (2), 237−50. (12) Ishikawa, J.; et al. Changes in the Ceramide Profile of Atopic Dermatitis Patients. J. Invest. Dermatol. 2010, 130, 2511−2514. (13) Das, C.; Noro, M. G.; Olmsted, P. D. Lamellar and Inverse Micellar Structures of Skin Lipids: Effect of Templating. Phys. Rev. Lett. 2013, 111, 148101. (14) Paudel, K. S.; Milewski, M.; Swadely, L. C.; Brogden, K. N.; Ghosh, P.; Stinchcomb, A. L. Challenges and Opportunities in Dermal/Transdermal Delivery. Ther. Delivery 2010, 1 (1), 109−131. (15) Farwanah, H.; Wohlrab, J.; Neubert, R. H.; Raith, K. Profiling of Human Stratum Corneum Ceramides by means of normal phase LC/ APCI-MS. Anal. Bioanal. Chem. 2005, 383, 632−637. (16) Norlén, L.; Nicander, I.; Lundsjö, A.; Cronholm, T.; Forslind, B. A New HPLC-Based Method for the Quantitative Analysis of Inner Stratum Corneum Lipids with Special Reference to the Free Fatty Acid Fraction. Arch. Dermatol. Res. 1998, 290, 508−516. (17) Stewart, M. E.; Downing, D. T. A New 6-Hydroxy-4Sphingenine-Containing Ceramide in Human Skin. J. Lipid Res. 1999, 40, 1434−1439. (18) Robson, K. J.; Stewart, M. E.; Michelsen, S.; Lazo, N. D.; Downing, D. T. 6-Hydroxy-4-Sphingenine in Human Epidermal Ceramides. J. Lipid Res. 1994, 35, 2060−2068. (19) Mojumdar, E. H.; Kariman, Z.; Kerckhove, L.; Gooris, G. S.; Bouwstra, J. A. The Role of Ceramide Chain Length Distribution on the Barrier Properties of the Skin Lipid Membranes. Biochim. Biophys. Acta, Biomembr. 2014, 1838 (10), 2473−83. (20) Holtje, M.; Forster, T.; Brandt, B.; Engels, T.; Rybinski, W. Molecular Dynamics Simulations of Stratum Corneum Lipid Models: Fatty Acids and Cholesterol. Biochim. Biophys. Acta, Biomembr. 2001, 1511, 156−167. (21) Pandit, S. A.; Scott, H. L. Molecular-Dynamics Simulation of a Ceramide Bilayer. J. Chem. Phys. 2006, 124 (1), 014708. (22) Notman, R.; Otter, W. K.; Noro, M. G.; Briels, W. J.; Jamshed, A. The Permeability Enhancing Mechanism of DMSO in Ceramide Bilayers Simulated by Molecular Dynamics. Biophys. J. 2007, 93, 2056−2068. (23) Das, C.; Olmsted, P. D.; Noro, M. G. Water Permeation through Stratum Corneum Lipid Bilayers from Atomistic Simulations. Soft Matter 2009, 5, 4549−4555. (24) Das, C.; Olmsted, P. D.; Noro, M. G. Simulation Studies of Stratum Corneum Lipid Mixtures. Biophys. J. 2009, 97, 1941−1951. (25) Guo, S.; Moore, T. C.; Iacovella, C. R.; Strickland, A. L.; 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. (26) Berger, O.; Edholm, O.; Jahnig, F. Molecular Dynamics Simulations of a Fluid Bilayer of Dipalmitoylphosphatidylcholine at Full Hydration, Constant Pressure, and Constant Temperature. Biophys. J. 1997, 72, 2002−2013. (27) Pastor, R. W.; MacKerell, A. D. Development of the CHARMM Force Field for Lipids. J. Phys. Chem. Lett. 2011, 2, 1526−1532. (28) Alwarawrah, M.; Dai, J.; Haung, J. A Molecular View of the Cholesterol Condensing Effect in DOPC Lipid Bilayers. J. Phys. Chem. B 2010, 114 (22), 7516−7523. (29) Hofsäß, C.; Lindahl, E.; Edholm, O. Molecular Dynamics Simulations of Phospholipid Bilayers with Cholesterol. Biophys. J. 2003, 84 (4), 2192−2206. (30) Mombelli, E.; Morris, R.; Taylor, W.; Fraternali, F. HydrogenBonding Propensities of Sphingomyelin in Solution and in a Bilayer Assembly: A Molecular Dynamics Study. Biophys. J. 2003, 84, 1507− 1517. (31) Ryckaert, J. P.; Bellemans, A. Molecular Dynamics of Liquid Butane near its Boiling Point. Chem. Phys. Lett. 1975, 30, 123−125. (32) Berendsen, H.; Postma, J.; Gunsteren, W.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces 1981, 331−342.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b02093. Figure S1 (equilibration of area per lipid), Figure S2 (equilibration of volume per lipid), Figure S3 (equilibration of potential energy), Figure S4 (structures of equilibrated bilayer at 300 K), Figure S5 (method of calculation of bilayer width), and Table S1 (effect of system size on the bilayer properties) (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +91-02066086203. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors would like to thank High Performance Computing at Tata Consultancy Services (TCS) for providing access to the EKA Super computer. We would also like to thank Mr. K Ananth Krishanan, CTO, TCS, and Dr. Pradip, Vice President and Head Process Engineering Innovation Laboratories, TCSTRDDC, Pune, for their constant encouragement and support during this project.



ABBREVIATIONS CER, ceramide; CHOL, cholesterol; FFA, free fatty acid; VPL, volume per lipid; APL, area per lipid; APC, area per ceramide



REFERENCES

(1) Michaels, A. S.; Chandrasekaran, S. K.; Shaw, J. E. Drug Permeation through Human Skin: Theory and In-Vitro Experimental Measurement. AIChE J. 1975, 21, 985−996. (2) Christian, U.; Hansen, C. M.; Dyk, V.; John, W.; Jensen, P. O. Permeability of Commercial Solvents through Living Human Skin. Am. Ind. Hyg. Assoc. J. 1995, 56 (7), 651−60. (3) Bartosova, L.; Bajgar, J. Transdermal Drug Delivery In-Vitro Using Diffusion Cells. Curr. Med. Chem. 2012, 19, 4671−4677. (4) Torin Huzil, J.; Sivaloganathan, S.; Kohandel, M.; Foldvari, M. Drug Delivery through the Skin: Molecular Simulations of Barrier Lipids to Design More Effective Noninvasive Dermal and Transdermal Delivery Systems for Small Molecules, Biologics, and Cosmetics. WIREs NanomeNanobiotechnol. 2011, 3, 449−462. (5) Jepps, O. G.; Dancik, Y.; Anissimov, Y. G.; Roberts, M. S. Modeling the Human Skin Barrier- Towards a Better Understanding of Dermal Absorption. Adv. Drug Delivery Rev. 2013, 65, 152−168. (6) Mitragotri, S.; et al. Mathematical Models of Skin Permeability: An Overview. Int. J. Pharm. 2011, 418 (1), 115−129. (7) Williams, M. L.; Elias, P. M. The Extracellular Matrix of Stratum Corneum: Role of Lipids in Normal and Pathological Function. CRC Crit. Rev. Ther. Drug Carrier Syst. 1987, 3, 95−122. (8) Wertz, P. W.; Downing, D. T. The Nature of the Epidermal Barrier: Biochemical Aspects. Adv. Drug Delivery Rev. 1996, 18, 283− 294. (9) Norlen, L.; Nicander, I.; Rozell, B. L.; Ollmar, S.; Forslind, B. Inter- and Intra-Individual Differences in Human Stratum Corneum Lipid Content Related to Physical Parameters of Skin Barrier Function In Vivo. J. Invest. Dermatol. 1999, 112, 72−77. (10) Weerheim, A.; Ponec, M. Determination of Stratum Corneum Lipid Profile by Tape Stripping in Combination with HighPerformance Thin-Layer Chromatography. Arch. Dermatol. Res. 2001, 293, 191−199. 11654

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655

Article

The Journal of Physical Chemistry B

(55) Los, D. A.; Murata, N. Membrane Fluidity and its Roles in the Perception of Environmental Signals. Biochim. Biophys. Acta, Biomembr. 2004, 1666, 142−157. (56) Mills, T. T.; Toombes, G. M.; Nagle, S. T.; Smilgies, D. M.; Feigenson, W. G.; Nagley, J. F. Order Parameters and Areas in FluidPhase Oriented Lipid Membranes Using Wide Angle X-Ray Scattering. Biophys. J. 2008, 95 (2), 669−681. (57) Smondyrev, A. M.; Berkowitz, M. L. Structure of Dipalmitoylphosphatidyl-choline/Cholesterol Bilayer at Low and High Cholesterol Concentrations: Molecular Dynamics Simulation. Biophys. J. 1999, 77 (4), 2075−2089. (58) Mizushima, H.; Fukusawa, 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, 361−367. (59) Egberts, E.; Marrink, S. J.; Berendsen, H. J. C. Molecular Dynamics Simulation of a Phospholipid Membrane. Eur. Biophys. J. 1994, 22, 423−436. (60) Marrink, S. J.; Berendsen, H. J. C. Simulation of Water Transport through a Lipid Membrane. J. Phys. Chem. 1994, 98 (15), 4155−4168. (61) Reddy, A. S.; Warshaviak, D. T.; Chachisvilis, M. Effect of Membrane Tension on the Physical Properties of DOPC Lipid Bilayer Membrane. Biochim. Biophys. Acta, Biomembr. 2012, 1818 (9), 2271− 2281. (62) Garidel, P. Structural Organisation and Phase Behaviour of a Stratum Corneum Lipid Analogue: Ceramide 3A. Phys. Chem. Chem. Phys. 2006, 8 (19), 2265−2275. (63) Schroeter, A.; Kiselevb, M. A.; Haußc, T.; Dantec, S.; Neubert, R. H. H. Evidence of Free Fatty Acid Interdigitation in Stratum Corneum Model Membranes Based on Ceramide [AP] by Deuterium Labelling. Biochim. Biophys. Acta, Biomembr. 2009, 1788, 2194−2203. (64) Marrink, S. J.; Mark, A. E. Effect of Undulations on Surface Tension in Simulated Bilayers. J. Phys. Chem. B 2001, 105, 6122−6127. (65) Rastogi, A.; Ghosh, A. K.; Suresh, S. J. Hydrogen Bond Interactions between Water Molecules in Bulk Liquid, Near Electrode Surfaces and Around Ions. Thermodynamics-Physical chemistry of aqueous systemm 2011, 291−306.

(33) Berendsen, H.; Spoel, D. V.; Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (34) Hess, B.; Kutzner, C.; Spoel, D. V.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (35) Pronk, S.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−54. (36) Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (37) Apol, E.; et al. GROMACS User Manual; 2006; Ver. 4.5.6. (38) Patra, M.; Karttunen, M.; Hyvonen, M. T.; Falck, E.; Lindqvist, P.; Vattulainenv, I. Molecular Dynamics Simulations of Lipid Bilayers: Major Artifacts Due to Truncating Electrostatic Interactions. Biophys. J. 2003, 84, 3636−3645. (39) Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30 (13), 2157−2164. (40) Martinez, J. M.; Martinez, L. Packing Optimization for Automated Generation of Complex System’s Initial Configurations for Molecular Dynamics and Docking. J. Comput. Chem. 2003, 24 (7), 819−825. (41) Falck, E.; Patra, M.; Karttunen, M.; Hyvonen, M.; Vattulainen, I. Lessons of Slicing Membranes: Interplay of Packing, Free Area, and Lateral Diffusion in Phospholipid/Cholesterol Bilayers. Biophys. J. 2004, 87, 1076−1091. (42) Chiu, S. W.; Jakobsson, E.; Mashl, R. J.; Scott, H. J. CholesterolInduced Modifications in Lipid Bilayers: A Simulation Study. Biophys. J. 2002, 83, 1842−1853. (43) Moore, D. J.; Rerek, M. E.; Mendelsohn, R. FTIR Spectroscopic Studies of the Conformational Order and Phase Behavior of Ceramides. J. Phys. Chem. B 1997, 101, 8933−8940. (44) Shah, J.; Atienza, J. M.; Rawlings, A. V.; Shipely, G. G. Structural and Thermotropic Properties of Synthetic C16:0 (Palmitoyl) Ceramide: Effect of Hydration. J. Lipid Res. 1995, 36, 1936−1944. (45) Shah, J.; Atienza, J. M.; Rawlings, A. V.; Shipely, G. G. Physical Properties of Ceramides: Effect of Fatty Acid Hydroxylation. J. Lipid Res. 1995, 36, 1945−1955. (46) Chen, H. C.; Mendelsohn, R.; Rerek, M. E.; Moore, D. J. Fourier Transform Infrared Spectroscopy and Differential Scanning Calorimetry Studies of Fatty Acid Homogeneous Ceramide 2. Biochim. Biophys. Acta, Biomembr. 2000, 1468, 293−303. (47) Kanicky, J. R.; Shah, D. Effect of Degree, Type, and Position of Unsaturation on the pKa of Long-Chain Fatty Acids. J. Colloid Interface Sci. 2002, 256, 201−207. (48) Hifeda, Y. M.; Rayfield, G. W. Phase Transitions in Fatty Acid Monolayers Containing a Single Double Bond in the Fatty Acid Tail. J. Colloid Interface Sci. 1985, 104, 209. (49) Craven, B. M. Crystal Structure of Cholesterol Monohydrate. Nature 1976, 260, 727−729. (50) Samby, J. M.; Momsen, M.; Kulkarni, V. S.; Brown, R. E. The Interfacial Elastic Packing Interactions of Galactosylceramides, Sphingomyelins and Phosphatidylcholines. Biophys. J. 1996, 70, 868−877. (51) Beare-Rogers, J. L.; Dieffenbacher, A.; Holm, J. V. Lexicon of Lipid Nutrition (IUPAC Technical Report). Pure Appl. Chem. 2001, 73 (4), 684−744. (52) Lewis, R. J. Hawly’s condensed chemical dictionary; John Wiley & Sons, Inc.: New York, 1997. (53) Dahlen, B.; Pasher, I. Molecular Arrangements in Sphingolipids. Thermotropic Phase Behavior of Tetracosanoylphytosphingoisne. Chem. Phys. Lipids 1979, 24, 119−133. (54) Demel, R. A.; Kruijff, B. D. The Function of Sterols in Membranes. Biochim. Biophys. Acta, Rev. Biomembr. 1976, 457, 109− 132. 11655

DOI: 10.1021/acs.jpcb.5b02093 J. Phys. Chem. B 2015, 119, 11643−11655