Molecular Dynamics Simulation Study of Skin Lipids: Effects of the

Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar Industrial Estate, Pune - 411013, India. J. Phys. Chem. B , 2015, 119 (35)...
1 downloads 13 Views 2MB Size
Subscriber access provided by Georgetown University | Lauinger and Blommer Libraries

Article

Molecular Dynamics Simulation Study of Skin Lipids: Effects of Molar Ratio of Individual Components Over a Wide Temperature Range Rakesh Gupta, and Beena Rai J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b02093 • Publication Date (Web): 14 Aug 2015 Downloaded from http://pubs.acs.org on August 18, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

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

The Journal of Physical Chemistry

Molecular Dynamics Simulation Study of Skin Lipids: Effects of 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

ABSTRACT Atomistic molecular dynamics (MD) simulations were employed to systematically investigate the effects of molar ratio of individual components cholesterol (CHOL), free fatty acid (FFA) and ceramides (CER) on the properties of skin lipid bilayer over a wide temperature range (300 - 400K). Several independent simulations were performed for bilayers comprising of only CER, CHOL or FFA molecules as well as those made up of mixture of CER: CHOL: FFA molecules different molar ratios. It was found that CHOL increases the stability of bilayer since the mixed (CER: CHOL: FFA) 1:1:0, 1:1:1 and 2:2:1 bilayer remained stable till 400 K while 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 mixed bilayer system. The CHOL molecule provided more rigidity to the mixed bilayer and lead to a more ordered phase at elevated temperatures. The CHOL molecule provided fluidity to bilayer below the phase transition temperature of CER and kept the bilayer rigid above the phase transition temperature. The FFA interdigitize with CER molecules and increases the thickness of 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 does forms inter and intra molecular hydrogen bonds while CHOL only forms intermolecular hydrogen bonds.

Keywords: Skin lipids, area and volume per lipid, phase transition, cholesterol interaction, hydrogen bond, area compressibility.

*To whom correspondence should be addressed. [email protected] ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 2 of 38

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 structure-based models5. 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 outer most part of skin epidermis, provides main barrier to penetration of pathogens and against water loss. The SC mainly comprises 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 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 presence of more than 300 different CER species comprising of varying head group 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) 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 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 micro needles but 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 behaviour and structure at 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 sphingo motif is always found to be

ACS Paragon Plus Environment

2

Page 3 of 38

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

The Journal of Physical Chemistry

between 16-18 carbons but the length of fatty acid varies from 16 to 34 carbons, thereby leading to asymmetry in the length of ceramides chains.17-18 Atomistic modeling of such a complex system incorporating all possible ceramides and fatty acids found in SC is far away from current computational capabilities. However, to in order to simulate a realistic SC layer19, we have chosen the most abundant ceramide, CER-NS 24:0 and free fatty acid, FFA 24:0. Holtze et al.20 first time reported the simulation of two component mixtures of fatty acid and cholesterol, Pandit and Scott21 studied single component CER-NS (16:0) bilayer, Notman et al.22 investigated effect of DMSO on pure CER-NS (24:0) bilayer. The first study of three component mixtures of CER NS (24:0), free fatty acid (24:0) and cholesterol was performed by Das et al. 23 Authors investigated effect of molar ratio of individual component on structural properties and followed by diffusion and permeability calculation for CER and mixed bilayer.24 Same authors have also reported on the modelling of 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 Berger26 and CHARMM27 force field. They found that OH group in CER-NP aligned itself perpendicular to the bilayer normal as compare to CER-NS, resulting into more hydrogen bonded network, which leads to high gel-to-liquid phase transition temperature in CER-NP. Guo et al. performed simulations for symmetric CER-NS bilayer which leads to a proper ordering of lipids tail in its interior, resulting into a high gel to liquid phase temperature. To the best of author’s knowledge no studies have been reported on the phase transition of mixed skin lipid bilayers. We report on the phase transition of skin lipids in at various compositions along with calculated area per CER molecule in multicomponent bilayers. The effects of CHOL and FFA in ceramide lipid bilayers over wide temperature range is systematically investigated via atomistic molecular dynamics (MD) simulation in constant NVT and NPT condition. In bilayer simulation the area per lipid is a primary parameter which is directly correlated to the molecular organisation of lipid bilayer.28 It can be compared with experiments as well as it gives information about the equilibration of the system in simulation. We have modified the method of Hofsäß et al. 29 to calculate

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

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

Page 4 of 38

the area per lipid in the mixed bilayer.

1. Simulation details 1.1. Force field The potential function and interaction parameters for CER were taken from Berger force field 26, which is properly parameterized for Palmitoylsphingomyelin (PSM) bilayers and phospholipids. The PSM’s have similar structure as CER-NS except for head group geometry. The head group size and partial charges in PSM are very high as compared to CER.30 The polar hydrogen’s 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 were same as used by Holtze et al.20 The simple point charge (SPC) model was used for water molecule.

32

The united atom carbon in all components had

zero partial charge.

1.2. Simulation setup All simulations were carried out in NVT and NPT ensemble using the latest GROMACS molecular dynamics package.33-35 The temperature was controlled by Berendsen thermostat with a time constant of 0.5 ps for initial equilibration, later it was changed to Nose-Hoover thermostat with same time constant. Pressure was maintained at 1 bar using Berendsen barostat for initial equilibration with a time constant of 0.5 ps and compressibility of 4.5 x 10-5 bar-1, later in production run it was changes to ParrinelloRahman barostat with a time constant of 5 ps and compressibility of 4.5 x 10-5 bar-1. The semi isotropic coupling (coupled separately in xy and z direction) of barostat was used to simulate tension less bilayer. All the bonds in lipid molecules were constrained using LINCS algorithm36 while for water SETTLE37 algorithm was used. A time step of 2 fs was used for simulation below 360 K and 1 fs for higher temperature. It has been shown earlier that cutoff for electrostatic interaction leads to some artifacts in phospholipids bilayer simulation.38 A cut off of 1.2 nm was used for van der Waals and electrostatic interactions. Long range electrostatic interaction was computed using particle mesh Ewald (PME) method. Equilibration time for all the systems was 10 ns (NVT) followed

ACS Paragon Plus Environment

4

Page 5 of 38

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

The Journal of Physical Chemistry

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 production run. All the properties were averaged over 5 blocks of 2 ns.

1.3. Initial bilayer structures Molecule structure of individual component used in this simulations, are shown in Figure1. We have used hairpin configuration of CER molecules.23 To remove bad contacts between molecules, energy minimization was carried out in vacuum using steepest decent and conjugate gradient method. For single component bilayer system, an individual molecule was first minimized and then replicated using “genconf” (Gromacs) in X and Y direction to obtain single layer. This single layer was then replicated in Z direction to obtain 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 10 ns NVT MD run at 300 K, by keeping lipid molecules fix for proper hydration. For mix bilayer, initial random configurations were generated using the tool- PACKMOL.39-40 Initial and final configurations of mix bilayers mostly correlated in gel phase simulations.20 We further performed simulated annealing to remove these artifacts. All the minimized structures were heated till 360 K and cooled back again to 300 K in 1 ns run. Structures obtained from the simulated annealing, were equilibrated for 100 ns at 300K. 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. Same procedure was followed for the simulations of pure CHOL and FFA bilayers. Table 1 shows the molar ratio used in the simulations and corresponding number of individual molecules. The ratios in the article should be considered in the order of CER: CHOL: FFA until unless it is specified specifically.

2. Area and Volume per molecule In an MD simulation of single component bilayer which has normal along the z direction, the volume per lipid (VPL) could be calculated by subtracting volume of water molecule from the total volume,

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

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

V  Lx Ly Lz

Vlipid 

Page 6 of 38

............ (1)

V  NW VW N lipid

........... (2)

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

Vcer 

V  NW VW  N cholVchol  N ffaV ffa N cer

…………..…. (3)

where Nffa and Nchol are number of FFA and cholesterol molecules in the mix bilayer respectively. Vchol, Vffa are volume per lipid for pure cholesterol and FFA bilayer respectively, which were calculated by separate individual simulations of pure cholesterol and FFA system in NPT ensemble at each temperature

using identical simulation

condition such as time steps and cut-offs. The area per lipid (APL) is a primary parameter which describes the molecular packing of lipid bilayer. It can be compared with experiments as well as it gives information about the equilibration of the system in simulation. In a molecular dynamics simulation of pure lipid bilayer which has normal along the z direction, the APL can be calculated using following equation

ACS Paragon Plus Environment

6

Page 7 of 38

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

The Journal of Physical Chemistry

A

2 Lx Ly

............... (4)

N lipid

Where Lx , Ly is the box length in X and Y direction respectively and N

lipid

is total

number of lipids in the bilayer. Some methods have been proposed for calculating the area per phospholipid molecule and area per CHOL in mix bilayer system at different temperature and CHOL mole fractions.41-42 Chiu et al.42 assumed area per lipid molecule to be independent of mole fraction of CHOL and calculated area per lipid by simple linear regression. Hofsäß et al.29 used volume per lipid and common bilayer thickness to compute individual area per molecule at different mole fraction of CHOL. We have employed similar approach to calculate area per molecule in our mix bilayer system. We have assumed that all the available volume was occupied by lipid molecules that mean there is no free volume inside the bilayer. Average thickness of a bilayer is given by:

havg 

V  NW VW A

…........ (5)

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

Acer 

Vcer havg

............ (6)

using equation (3), (5) and (6)

Acer 

 V  NW VW  N cholVchol  N ffaV ffa  A   V  NW VW  N cer 

Acer  Achol  Affa 

N chol Vchol  N ffaV ffa  A  1   N cer  V  NW VW 

AVchol V  NW VW

……... (7)

............ (8)

AV ffa

…....... (9)

V  NW VW

Where Achol and Affa are the area per molecule of pure cholesterol and FFA bilayer

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

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

Page 8 of 38

respectively. Ncer , Nchol and Nffa are the number of ceramide, cholesterol and FFA molecules in the mix bilayer. Vchol , Vffa are the volume per molecule calculated from the pure CHOL and FFA bilayer simulation.

3. Results and discussion 3.1. Heating run Initially equilibrated bilayer structures at 300 K were used for the heating run. All the bilayers were heated from 300K to 450K slowly at the rate of 5K/ns with 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 rigid structure of the CHOL the pure CHOL bilayer was stable for whole heating run and pure FFA disintegrated near 380 K. The bilayers disintegrated in the order of (CER: CHOL: FFA) 0:0:1, 1:0:1, 1:1:0, 1:0:0, 1:1:1, 2:2:1, 0:1:0.

3.2. Structures Figure 3 shows snapshot of pure CER bilayer and tails arrangement at each temperature. At physiological temperatures bilayer was found to be in gel phase with proper ordering of lipid tails in both of the leaflets. The lipids tails were found to be in hexagonal phase which are 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 centre of bilayer almost liquid like phase was observed due to asymmetric chain length of CER tails (Figure 3). Figure 4 and 5 show the individual component structure of 1:1:1 and 2:2:1 bilayer at each temperature. At lower temperature (< 340 K), all the components are in proper order for both the systems. Since FFA has lower melting point (~ 357 K) and CER has gel to liquids phase transformation temperature of ~ 365 K, near the 380 K more disordering is observed in both the systems. Due to its rigid nature, CHOL provides stability to bilayer, which is apparent from the structure of CHOL in both the systems at each temperature.

ACS Paragon Plus Environment

8

Page 9 of 38

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

The Journal of Physical Chemistry

3.3. Volume per lipid The pure bilayer of FFA and CHOL were simulated to get the volume per FFA and CHOL molecules and the results are presented in Table 2. We have obtained volume and area per FFA (APF) of 0.604 nm3 and 0.206 nm2, respectively for a pure FFA system, which is in close agreement with the earlier experimental value47-48 of 0.20-0.21 nm2 and simulation studies20 of palmitic fatty acid at 303 K. We have obtained area per cholesterol molecule of 0.365 nm2 at 300 K in pure cholesterol system, which is closer to experimental value49-50 of 0.36-0.37 nm2 and simulation value

20

2

of ~0.37 nm . The melting point of

FFA (lignoceric acid, 24:0) and cholesterol are found to be around ~ 357 K51 and 423 K52, respectively. We also noticed that pure FFA bilayer system was not stable at temperature more than 370 K. Volume per CER molecule calculated using equation 3 for each bilayer system has shown in Figure 6a. The Vcer was least in 1:1:0 bilayer and it increases in presence of FFA.

3.4. Area per lipid Figure 6b shows area per CER (Acer) molecule for each bilayer at different temperature calculated using equation 7. The APL increases with temperature because of more thermal energy which loses the packing of lipid tails. Acer value of ~0.39 nm2 for pure CER bilayer at 300K matches very well with x-ray experiment value of ~0.4 nm2 obtained by Dahlen et al.53 as well as literature.23,24 Acer increases gradually with temperature in all systems. The Acer values are found to be higher in case of mixed 2:2:1 bilayer and lower for 1:0:1 bilayer at each temperature studied. Figure 6c shows area spanned by each cholesterol (Achol ) molecule in the bilayer, calculated from equation 8. It is interesting to note that CHOL spanned less area in each of the mixed bilayer as compare to pure cholesterol bilayer. Figure 6d shows area per free fatty acid (Affa) molecule calculated by equation 9. We found that more of the cholesterol present in the bilayer leads to more of Acer and Affa values. Overall effect of CHOL is to provide more area to CER and FFA in mixed bilayer system.

3.5. Tail order parameter

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry

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

Page 10 of 38

This parameter is used to characterize the order in lipid bilayers and can be measured by deuterium NMR. An order parameter for every CH2 group in the chains is defined as:

SZ 

3 1 cos 2 d    2 2

............ (10)

where dθ is the angle between the bilayer normal (z axis) and vector connecting cn-1 to cn+1 CH2 atom. The value of Sz = 1 represents tails are oriented parallel to bilayer axis, while a value of -0.5 represents orientation perpendicular to bilayer normal. Figure 7 shows the comparison of order parameter in each bilayer system below and above the phase transition temperature. The profiles for all the system are almost similar but order parameter decreases with temperature. Chain sn1 (C1-C16) has little bit more ordering as compare to chain sn2 (C27-C49). The order parameter of chains are very low near the interface (C24 and C16) and increases as we move further down to high lipid density region and decreases towards the centre of the bilayer. The asymmetry in ceramide tails leads to more disordering of sn2 chain length in low lipid density region. The similar chain length of FFA leads to a proper inter-digitization thereby increasing the ordering of ceramide tails. CHOL induces fluidity in phospholipid lipid membranes by decreasing the chain ordering below the gel to liquid crystal transition temperature while it increases the rigidity of bilayers above the transition temperature. This phenomenon has been observed both in experiments54-55 and simulations studies.56-57 Mizushima et al. 58

performed differential scanning calorimetry (DSC) experiments on ternary mixture of

CER, FFA and CHOL and found that CHOL increases the fluidity of alkane chains by reducing ordering. This phenomenon has also been observed in recent simulation study.22 We have noticed both these effects in our simulations as shown in Figure 7. Below the phase transition temperature of CER (~365 K), the order parameter is found to be little bit higher in 1:1:1 system as compare to 2:2:1 system. While at higher temperature (>~365 K) reverse phenomena is observed. The order parameter is little bit more in 2:2:1 system as compare to 1:1:1 system. The order parameter for 1:0:1 is highest for both of the chains at each temperature. The same effect of CHOL is also seen in fatty acid chain as shown in Figure 8. It follows the similar trends as CER chains, below the melting point of FFA (~357 K), the order parameter is higher in 1:1:1 system while above this temperature

ACS Paragon Plus Environment

10

Page 11 of 38

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

The Journal of Physical Chemistry

reverse phenomena 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. 3.6. Lipid density In Figure 9, the densities of lipid components are plotted along the bilayer normal (zaxis) for ceramide bilayer at different temperature. The profile is similar to the four region model as proposed by Marrink et al.59-60 Density from z ~ -6.0 (not shown in 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 head groups density increases ( z~ 3.2 to z ~ -2.0). Next region from z ~ -2.0 to z ~ - 0.8, has very ordered and tightly packed lipid tails, which is mainly responsible for barrier properties of skin. In the last region total lipid density has 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 group sit just below the water-lipid interface. The interfacial width is found to be smaller than those of phospholipids membrane because of small size and low partial charge of head group. Next to this region, lipid tail density is very high because of proper ordering of tails, the density in this region is found to be around of ~ 1.3 g/cm3. In the middle of bilayer there is dip in lipid density because of asymmetric chain length, this region has more and more free volume and has density of ~ 0.7 g/cm3. The temperature induces disorder in the system which leads to a small tail density in liquid crystalline phase (~360-370 K) as compare to gel phase (300K). The head group density profile gets broader with temperature because of more fluidity in liquid crystalline phase. Figures 10 shows the density profiles of each component of mix 1:1:1 and mix 2:2:1 bilayer. The density profile is almost similar at all the temperature studied. It is worth noting that the peaks of CER and FFA density are near the interface indicating that head groups of both of the component sit just below the lipid -water interface. The CHOL mostly sits in highly dense region of bilayer. Due to smaller size and low partial charge of head groups in ceramide bilayer as compare to phospholipid CHOL molecules sit much below the head group region, almost between the lipid tails. The CER and FFA profile

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

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

Page 12 of 38

shows a small peak at the centre of lipid bilayer at low temperature which was absent in pure ceramide bilayer. This confirms that similar chain length of CER and FFA leads to proper inter-digitization. These peaks get distorted with increase in temperature because higher thermal energy leads to disordering in the tails enabling fluid like environment. The lipid tail density decreases with increase in temperature because bilayer moves from gel to liquid phase at higher temperature. Figures 11a and 11b show density profile of individual components along the bilayer normal in each bilayer at 300 K and 360 K, respectively. Overlapping of FFA and ceramide density peak and inter-digitization between them can easily be seen. FFA induces more and more ordering in lipid tails and increase the thickness of bilayer but CHOL reduces the ordering as well as thickness of bilayer.16 The lipid density profile shows the similar behaviour at all the temperature studied. Peak to peak distance between lipid density was found to be highest in ceramide bilayer and lowest in mix 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 temperature. The shape of electron density in individual bilayer is almost similar for all the temperature but density gets lower with increase in temperature. It is interesting to note the mixed bilayers which comprise CHOL have two peaks while single peak was observed in CHOL free mixed bilayer. This confirms that CHOL mostly sits in between the CER and FFA chains.

3.7. Interfacial and bilayer width The interface (bilayer/water) width was estimated as the distance over which water density drops to 1/e of its bulk value.22,

23

The bilayer thickness could be defined in

several ways such as distance between phosphate atom in each leaflet in phospholipid membrane61, or distance between peaks of electron density profile.62 As shown in Figures 10 and 11, since CHOL sits in the interior space of bilayer, it is not convenient to calculate bilayer thickness by calculating distance between some of the head group atoms. We have used two different methods to calculate bilayer thickness. Schematic of the first approach has shown in supplementary information. Second approach is based on

ACS Paragon Plus Environment

12

Page 13 of 38

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

The Journal of Physical Chemistry

volume per CER and area per CER using equation 6. Results for bilayer thickness of each system estimated by 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 temperature studied. The bilayer thickness in mixed bilayer was found to be lowest in 1:1:0 bilayer and maximum in 1:0:1 bilayer. This shows that the same chain length fatty acid interdigitize with ceramide molecules chain and increases the chain length on other side being small in size CHOL reduces the bilayer thickness. This phenomenon has also been observed experimentally by Schroeter et al.63 The bilayer thickness decreases with increasing temperature because percentage in area expansion is larger than that of volume. The interfacial width (~0.5 nm) is smaller than that of sphingomyelin bilayers (~0.9 nm) because of small head group of 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. Minimum distance was found in 1:1:0 bilayer while maximum distance was found in 1:0:1 system. Another interesting thing we noticed that above the phase transition temperature of CER, in mixed bilayer distances between these atoms in two leaflets were more in high CHOL containing bilayer which supports the fact that above the phase transition temperature CHOL provides more ordering to the bilayer.

4.7 Area Compressibility Area compressibility of a bilayer whose normal oriented along the Z axis is calculated as K A  kbT

 A  A    A 2 2

………………………………………….. (12)

where A = Lx x Ly, is area of the bilayer in XY plane, kB is Boltzmann constant and T is temperature. The angular brackets denote the ensemble averages taken over the course of simulation Figure 14 shows the area compressibility of each bilayer at different temperature. 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. Presence of FFA reduces the area

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

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

Page 14 of 38

compressibility and makes bilayer soft in nature. The compressibility increases with 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 (supplementary information). The KA decreases with the increase in system size in tensionless bilayer system. This has also been observed in earlier studies of phospholipids bilayer.64

3.8 Hydrogen bond The numbers of hydrogen bonds present in bilayer systems were calculated based on the geometrical criteria. A hydrogen bond exists if the distance between the donor and acceptor atom is ≤ 0.35 nm and acceptor-donor -hydrogen angle is ≤ 60◦. Table 4 shows the number of hydrogen bonds between individual components of the bilayer at different temperature. Earlier study shows that the number of hydrogen bonds per water molecule in liquid water is around 3.5-3.8 at 300K.65 In all lipid bilayer system this number decreases to between 1.5-1.7 at 300K since at the interface CER forms the hydrogen with water. Being having small head group and small partial charge, CHOL is unable to form hydrogen bond with other CHOL molecules although it does forms hydrogen bonds with water molecules at interface. CER molecules form more intramolecular hydrogen bonds as compare to CHOL and FFA. The number of hydrogen bonds decreases for each bilayer system with increase in temperature because of high kinetic energy (thermal fluctuation) at elevated temperature.

Conclusion The effect of temperature on skin lipid bilayers has been studied at molecular scale. We have performed molecular dynamics simulations for different composition ratio of main constituent (CER-NS 24, CHOL and FFA 24:0) of skin's stratum corenum lipid matrix. We have found that CER bilayer exhibits hexagonal gel phase at normal skin temperature which gets converted into liquid crystalline phase at ~365 K. The asymmetric chain length of CER tails induces fluid like environment in the centre of bilayers. The presence of symmetric chain length FFA induces more and more inter-digitization and increases the ordering of CER tails. The CHOL molecule sits in the interior of lipid tails and

ACS Paragon Plus Environment

14

Page 15 of 38

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

The Journal of Physical Chemistry

provides more rigidity to bilayer at higher temperature. CHOL molecules shows duel nature, it decreases the chain ordering mixed bilayer below the phase transition temperature of CER and increases the ordering above the phase transition temperature. The fatty acids of similar chain length of CER, increases the thickness of the bilayer while short and rigid CHOL molecule decreases the thickness of the bilayer. Presence of CHOL increases the area compressibility of bilayer which is responsible for the barrier function of skin while FFA reduces the compressibility. CER molecule forms intra and inter molecule hydrogen bonds while CHOL only forms inter molecule hydrogen bonds. In real, human skin lipids comprise of heterogeneous mixture of poly disperse 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.

Acknowledgement Authors would like to thank High Performance Computing at Tata Consultancy Services (TCS) for providing access to 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 Labs, TCS-TRDDC, Pune for their constant encouragement and support during this project.

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

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

Page 16 of 38

Figures

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 represents lipid tails of CER and FFA respectively. Colour red, blue, white and cyan represents oxygen O, nitrogen N, hydrogen H and carbon C (united atom) respectively.

ACS Paragon Plus Environment

16

Page 17 of 38

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

The Journal of Physical Chemistry

Figure 2. Evolution of a) temperature b) area per lipid and c) box volume with simulation time in heating run. Colour black, red, green, blue, yellow, cyan and magenta represents system at composition (CER: CHOL: FFA) 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1 and 2:2:1 respectively. Please refer to online article for the colour code.

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

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

Page 18 of 38

Figure 3. Snapshot of pure ceramide bilayer at different temperature (top), arrangement of ceramide tails at plane ~ 1.2 nm above the bilayer centre (below). Ceramide head group, tails and water are shown in blue, cyan and red respectively.

Figure 4. Snapshots of individual components and lateral packing of lipid tails in the upper leaflet (at plane ~ 1.4 nm above the bilayer centre) of mixed (CER: CHOL: FFA) 1:1:1 bilayer at different temperature. From top snapshots are at 300 K, 340K, 360K,

ACS Paragon Plus Environment

18

Page 19 of 38

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

The Journal of Physical Chemistry

370K, 380K and 400K respectively. CER, CHOL and FFA are shown in light green, pink and yellow colour respectively

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

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

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

Page 20 of 38

Figure 6. a) volume per CER Vcer molecule calculated by equation 3 b) area per CER Acer molecule calculated by equation 7 c) area per cholesterol AChol molecule calculated by equation 8 d) Area per free fatty acid Affa molecule calculated by equation 9 in each system at different temperature. The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1 and 2:2:1 represents the molar ratio of CER: CHOL: FFA in the bilayer.

ACS Paragon Plus Environment

20

Page 21 of 38

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

The Journal of Physical Chemistry

Figure 7. Tail order parameter Sz of CER sn1 and sn2 chains at 300K and above the phase transition temperature of 370K. The chain 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 represents the molar ratio of CER: CHOL: FFA in the bilayer.

Figure 8 Order parameter (Sz) of FFA chain in mixed 1:1:1 and 2:2:1 bilayer at different

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry

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

Page 22 of 38

temperature. FFA chain is shown in Figure 1. The 1:1:1 and 2:2:1 represents the molar ratio of CER: CHOL: FFA in the bilayer.

Figure 9. Density of individual component of pure ceramide bilayer along the bilayer normal at 300 K (from top), 340 K, 360 and 370 K. System, water, head group, tail and ceramide densities are shown in brown, black, cyan, violet and magenta colour respectively. Please refer to online article for the colour code.

ACS Paragon Plus Environment

22

Page 23 of 38

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

The Journal of Physical Chemistry

Figure 10. Density of individual component of mixed (CER:CHOL:FFA) 1:1:1 and 2:2:1 bilayer along the bilayer normal at 300 K (from top), 340 K, 360K , 380 K and 400 K. CER, CHOL, FFA, lipid, water and system densities are shown in black, violet, cyan, magenta , yellow and red colour respectively. Please refer to online article for the colour code.

ACS Paragon Plus Environment

23

The Journal of Physical Chemistry

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

Page 24 of 38

Figure 11. Density of a) CER b) CHOL c) FFA and d) lipid along the bilayer normal in each bilayer system at 300 K (left) and 360 K (right). Please refer to online article for the colour code.

ACS Paragon Plus Environment

24

Page 25 of 38

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

The Journal of Physical Chemistry

Figure 12. Electron density of each bilayer system along the bilayer normal z at different temperature. The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1 and 2:2:1 represents the molar ratio of CER: CHOL: FFA in the bilayer. Please refer to online article for the colour code.

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry

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

Page 26 of 38

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

Figure 14. Dependence of Area compressibility kA over the temperature 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 represents the molar ratio of CER: CHOL: FFA in the bilayer.

ACS Paragon Plus Environment

26

Page 27 of 38

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

The Journal of Physical Chemistry

Tables Table1. Number of individual molecules used in simulation for each molar ratio. System*

CER

CHOL

FFA

Water

1:0:0

128

0

0

5120

0:1:0

0

128

0

5120

0:0:1

0

0

512

9020

1:1:0

64

64

0

5120

1:0:1

84

84

168

5120

1:1:1

52

50

52

5120

2:2:1

56

56

32

5120

*The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1 and 2:2:1 represents the molar ratio of CER: CHOL: FFA in the bilayer.

Table 2. Volume# and area$ per FFA/CHOL molecule at different temperature calculated from pure FFA (0:0:1) and CHOL (0:1:0) bilayer simulations respectively. System

Nlipids

Nwater

Temp

Vw*

Vchol, Vffa

Affa, Achol

K

nm3

nm3

nm2

FFA

512

9020

300

0.0302

0.604

0.206

FFA

512

9020

340

0.0313

0.634

0.217

FFA

512

9020

360

0.0319

0.643

0.221

FFA

512

9020

370

0.0323

0.663

0.253

CHOL

128

5120

300

0.0302

0.626

0.366

CHOL

128

5120

340

0.0313

0.631

0.370

CHOL

128

5120

360

0.0319

0.632

0.372

CHOL

128

5120

370

0.0323

0.636

0.375

CHOL

128

5120

380

0.0327

0.638

0.377

CHOL

128

5120

400

0.0336

0.640

0.380

*The water volume Vw was calculated from the separate simulations of 2048 water molecules at the same simulation condition used for respective bilayer system. # calculate using equation 2. $ calculated using equation 4.

ACS Paragon Plus Environment

27

The Journal of Physical Chemistry

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

Page 28 of 38

Table3. Distances between head group atoms# in between the two leaflets of the bilayer and average bilayer thickness of each bilayer at different temperature. System

havg1 & (method 1) Nm

havg2 & ( method 2) Nm

Temp

CER O21*

CHOL O6*

FFA O1*

K

Nm

Nm

Nm

1:0:0

300

5.182

-

-

5.65

5.55

1:0:0

340

5.136

-

-

5.61

5.53

1:0:0

360

5.135

-

-

5.59

5.52

1:0:0

370

5.123

-

-

5.57

5.44

0:1:0

300

-

3.378

-

3.423

3.36

0:1:0

340

-

3.367

-

3.411

3.342

0:1:0

360

-

3.355

-

3.401

3.324

0:1:0

370

-

3.354

-

3.396

3.321

0:1:0

380

-

3.346

-

3.389

3.319

0:1:0

400

-

3.326

-

3.372

3.318

0:0:1

300

-

-

5.635

5.866

5.796

0:0:1

340

-

-

5.54

5.847

5.787

0:0:1

360

-

-

5.528

5.826

5.726

0:0:1

370

-

-

5.497

5.240

5.612

1:1:0

300

4.413

4.1

-

4.693

4.583

1:1:0

340

4.372

4.089

-

4.683

4.573

1:1:0

360

4.305

4.083

-

4.644

4.554

1:1:0

370

4.263

4.001

-

4.645

4.545

1:1:0

380

4.240

3.911

-

4.631

4.511

1:1:0

400

4.101

3.815

-

4.662

4.492

ACS Paragon Plus Environment

28

Page 29 of 38

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

The Journal of Physical Chemistry

1:0:1

300

5.275

-

5.827

5.926

5.896

1:0:1

340

5.176

-

5.711

5.829

5.879

1:0:1

360

5.171

-

5.747

5.888

5.858

1:0:1

370

5.156

-

5.711

5.856

5.836

1:1:1

300

4.681

4.334

5.005

5.119

5.095

1:1:1

340

4.592

4.27

4.962

5.071

4.98

1:1:1

360

4.541

4.257

4.913

5.033

4.94

`1:1:1

370

4.491

4.089

4.925

4.962

4.90

1:1:1

380

4.355

3.96

4.662

5.023

4.85

1:1:1

400

4.179

3.86

4.578

4.889

4.85

2:2:1

300

4.542

4.186

4.831

4.957

4.83

2:2:1

340

4.470

4.166

4.838

4.911

4.80

2:2:1

360

4.425

4.125

4.823

4.929

4.78

2:2:1

370

4.459

4.068

4.729

4.859

4.77

2:2:1

380

4.425

4.021

4.704

4.829

4.76

2:2:1

400

4.325

3.922

4.645

4.758

4.66

*CERO21, CHOLO6 and FFAO1 represents the atom name O21, O6 and O1 in ceramide, cholesterol and fatty acid molecules respectively. Atoms are shown in the Figure 1. # The error in the data is around of ± 0.05% (T < 350 K) to ±0.15 % (T > 360 K). & The havg1 and havg2 are calculated from volume to area per lipid ratio and drop in water density profile respectively. Details of the methods are given in the supplementary information (figure S5).

ACS Paragon Plus Environment

29

The Journal of Physical Chemistry

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

Page 30 of 38

Table4. Number of hydrogen bonds* formed in different components of each bilayer system at different temperature. System$

Temp

CERCER

CERCHOL

CERFFA

CERW#

CHOLCHOL

CHOLFFA

CHOLW

FFAFFA

FFAW

W-W

300

1.52

0

0

2.13

~0

0

0

0

0

1.70

340

1.42

0

0

2.01

~0

0

0

0

0

1.60

360

1.36

0

0

1.95

~0

0

0

0

0

1.55

370

1.34

0

0

1.90

~0

0

0

0

0

1.51

380

1.27

0

0

1.88

~0

0

0

0

0

1.50

300

0

0

0

0

~0

0

2.03

0

0

1.70

340

0

0

0

0

~0

0

1.92

0

0

1.60

360

0

0

0

0

~0

0

1.86

0

0

1.55

370

0

0

0

0

~0

0

1.82

0

0

1.52

380

0

0

0

0

~0

0

1.78

0

0

1.49

400

0

0

0

0

~0

0

1.71

0

0

1.44

300

0

0

0

0

~0

0

0

0.28

1.66

1.67

340

0

0

0

0

~0

0

0

0.27

1.56

1.57

360

0

0

0

0

~0

0

0

0.24

1.49

1.52

370

0

0

0

0

~0

0

0

0.20

1.45

1.49

300

0.78

0.293

0

3.35

~0

0

0.74

0

0

1.69

340

0.75

0.280

0

3.00

~0

0

0.64

0

0

1.59

360

0.72

0.272

0

2.94

~0

0

0.56

0

0

1.54

370

0.68

0.265

0

2.82

~0

0

0.50

0

0

1.51

380

0.66

0.256

0

2.72

~0

0

0.48

0

0

1.49

400

0.64

0.247

0

2.51

~0

0

0.43

0

0

1.43

300

1.24

0

0.59

2.0

~0

0

0

0

0.26

1.71

K 1:0:0

0:1:0

0:0:1

1:1:0

1:0:1

ACS Paragon Plus Environment

30

Page 31 of 38

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

The Journal of Physical Chemistry

1:1:1

2:2:1

340

1.19

0

0.57

1.84

~0

0

0

0

0.24

1.61

360

1.14

0

0.54

1.75

~0

0

0

0

0.23

1.57

370

1.08

0

0.50

1.64

~0

0

0

0

0.20

1.54

300

0.76

0.31

0.19

2.93

~0

0.12

0.88

1.06

1.78

1.69

340

0.72

0.28

0.15

2.59

~0

0.10

0.79

1.44

1.76

1.59

360

0.71

0.25

0.14

2.56

~0

0.10

0.77

1.61

1.71

1.54

370

0.73

0.23

0.12

2.49

~0

0.10

0.76

1.54

1.69

1.51

380

0.66

0.23

0.14

2.43

~0

0.08

0.72

1.68

1.63

1.48

400

0.65

0.22

0.13

2.38

~0

0.07

0.69

1.62

1.61

1.46

300

0.80

0.34

0.13

2.72

~0

0

0.86

0.85

1.81

1.67

340

0.75

0.33

0.11

2.68

~0

0

0.83

0.71

1.75

1.59

360

0.71

0.33

0.09

2.63

~0

0

0.76

0.61

1.68

1.54

370

0.69

0.29

0.08

2.61

~0

0

0.74

0.62

1.67

1.52

380

0.65

0.26

0.08

2.59

~0

0

0.72

0.65

1.66

1.48

400

0.63

0.25

0.07

2.47

~0

0

0.63

0.59

1.57

1.44

* The error in the data is around of ±1 % (T < 350 K) - ±3 % (T > 360 K) # W stands for the water molecule. $ The 1:0:0, 0:1:0, 0:0:1, 1:1:0, 1:0:1, 1:1:1 and 2:2:1 represents the molar ratio of CER: CHOL: FFA in the bilayer.

Supplementary Information 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). This material is available free of charge via the Internet at http://pubs.acs.org.

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

ACS Paragon Plus Environment

31

The Journal of Physical Chemistry

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

Page 32 of 38

AUTHOR INFORMATION Corresponding Author Email: [email protected]

Phone: +91-02066086203

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. American Industrial Hygiene Association Journal. 1995, 56(7), 651-60 3) Baratosova, L.; Bajgar, J. Transdermal Drug Delivery In-Vitro Using Diffusion Cells. Current Medical Chemistry. 2012, 19, 4671-4677 4) Huzil, J. T.; 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. DOI: 10.1002/wnan.147 5) Jepps, O. G.; Dancik, Y.; Anissimov, Y. G.; Roberts, M. S. Modeling the Human Skin Barrier- Towards a Better Understanding of Dermal Absorption. Advanced Drug Delivery reviews. 2013, 65, 152-168 6) Mitragotri, S. et al. Mathematical Models of Skin Permeability: An Overview. International Journal of Pharmaceutics. 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 Critical Reviews in Therapeutic Drug Carrier Systems.1987, 3, 95-122 8) Wertz, P. W.; Downing, D. T. The Nature of the Epidermal Barrier: Biochemical Aspects. Adv. Drug Del. Rev. 1996, 18, 283-294

ACS Paragon Plus Environment

32

Page 33 of 38

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

The Journal of Physical Chemistry

9) Norlen, L.; Nicander, I.; Rozell, B. L.; Ollmar, S.; Forslind, B. Inter- and IntraIndividual 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

High-Performance

Thin-Layer

Chromatography. Arch. Dermatol. Res. 2001, 293, 191-199 11) Notman, R.; Anwar, J. Breaching the Skin Barrier--Insights from Molecular Simulation of Model Membranes. Adv. Drug. Deliv. Rev. 2012, 65(2), 237-50 12) Ishikawa, J. et al. Changes in the Ceramide Profile of Atopic Dermatitis Patients. J. Invest. Derm. 2010, 130, 2511-2514 13) Das, C.; Massimo, G. N.; Olmsted, P. D. Lamellar and Inverse Micellar Structures of Skin Lipids: Effect of Templating. Physical Review Letters. 2013, 111, 148101 14) Paude, K. S.; Milewski, M.; Swadely, L. C.; Brogden, K. N.; Ghosh, P.; Stinchcomb, A. L. Challenges and Opportunities in Dermal/Transdermal Delivery Ther Deliv. 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 HPLCBased 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-4-Sphingenine-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. 6Hydroxy-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. Biochimica et Biophysica Acta. 2014, 1838(10), 2473-83.

ACS Paragon Plus Environment

33

The Journal of Physical Chemistry

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

Page 34 of 38

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. 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, 45494555. 24) Das, C.; Olmsted, P. D.; Noro, M. G. Simulation Studies of Stratum Corneum Lipid Mixtures. Biophysical Journal. 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), 75167523 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. Hydrogen-Bonding 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.

ACS Paragon Plus Environment

34

Page 35 of 38

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

The Journal of Physical Chemistry

32) Berendsen, H.; Postma, J.; Gunsteren, W.; Hermans J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces. 1981, 331-342 33) Berendsen, H.; Spoel D. V.; Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comp. Phys. Comm. 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, 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. Comp. 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.

Journal of

Computational Chemistry. 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. Journal of Computational Chemistry. 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. Cholesterol-Induced Modifications in Lipid Bilayers: A Simulation Study. Biophys. J. 2002, 83, 18421853 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.

ACS Paragon Plus Environment

35

The Journal of Physical Chemistry

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

Page 36 of 38

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, 19451955. 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. 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. Journal of colloid and interface science. 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. Biochemistry. 1996, 35, 5696. 51) Rogers J. B.; Dieffenbacher, A.; Holm, J. V. Lexicon of Lipid Nutrition (IUPAC Technical Report). Pure and applied chemistry. 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.1976, 457, 109-132.

ACS Paragon Plus Environment

36

Page 37 of 38

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

The Journal of Physical Chemistry

55) Dmitry, A. L.; Murata, N. Membrane Fluidity and its Roles in the Perception of Environmental Signals. Biochimica et Biophysica Acta (BBA) – Biomembranes. 2004, 1666, 142-157. 56) Thalia, T. M.; Toombes, G. M.; Nagle, S. T.; Smilgies, D. M.; Feigenson, W. G.; Nagley, J. F. Order Parameters and Areas in Fluid-Phase Oriented Lipid Membranes Using Wide Angle X-Ray Scattering. Biophysical Journal. 2008, 95(2), 669- 681. 57) Smondyrev, A. M.; Berkowitz, M. L. Structure of Dipalmitoylphosphatidylcholine/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 Pseudo-Ceramide: A Study of the Function of Cholesterol. Journal of Lipid Research. 1996, 37, 361-367. 59) Egberts, E.; Marrink, S. J.; Berendsen, H. J. C. Molecular Dynamics Simulation of a Phospholipid Membrane. Eur Biophys J. 1993, 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. Biochimica et Biophysica Acta (BBA) Biomembranes. 2012, 1818(9), 2271-2281. 62) Garidel, P. Structural Organisation and Phase Behaviour of a Stratum Corneum Lipid Analogue: Ceramide 3A. Physical Chemistry Chemical Physics. 2006, 8(19), 2265-2275. 63) Schroetera, 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. Biochemica et. Biophysica Acta. 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.

ACS Paragon Plus Environment

37

The Journal of Physical Chemistry

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

Page 38 of 38

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.

Table of Contents

ACS Paragon Plus Environment

38