Molecular Dynamics Simulation of Skin Lipids - ACS Publications

Nov 21, 2016 - Engineering & Physical Sciences, TCS Research, Tata Research Development and Design Centre, Tata Consultancy Services, 54B,. Hadapsar I...
2 downloads 15 Views 2MB Size
Subscriber access provided by UIC Library

Article

Molecular Dynamics Simulation of Skin Lipids : Effect of Ceramide Chain Lengths on Bilayer Properties Rakesh Gupta, Balarama Sridhar Dwadasi, and Beena Rai J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b08059 • Publication Date (Web): 21 Nov 2016 Downloaded from http://pubs.acs.org on November 25, 2016

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 25

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 of Skin Lipids : Effect of Ceramide Chain Lengths on Bilayer Properties Rakesh Gupta, Balarama Sridhar D and Beena Rai* Engineering & Physical Science, TCSTM Research Tata Research Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar Industrial Estate, Pune – 411013, INDIA

*Corresponding author: beena.rai@tcs.com Fax:

91-20-66086399

Tel:

91-20-66086203

Abstract The Stratum Corneum (SC), outer layer of skin serves as a barrier for pathogens and maintains the trans-epidermal water loss. The lipid matrix of SC is the major diffusion-rate limiting pathway as molecules will have to pass through it. Ceramides play a key role in structuring, and, maintaining the barrier function of the skin. In this study, atomistic molecular dynamics (MD) simulations were employed to systematically investigate the effects of chain length of ceramide tails on the properties of bilayer at skin temperature of 310 K. The barrier properties were examined by means of permeation studies of water through the model membrane using steered MD simulations. Our studies revealed that shorter chains of one leaflet of bilayer do not interdigitate with the chains of other leaflet and lead to more free space in the middle of bilayer, thereby, leading to higher permeability. The dissimilar chain length of CER interdigitate with the other leaflet lipids and increases the electron density in the middle of the bilayer. The bilayer thickness increases with increase in CER chain length. The permeability of smaller chain ceramides was found to be an order of magnitude higher than the longer chain ceramides.

Introduction The human skin which is the largest organ of human body, protects it from the loss of water, attack of foreign pathogens and penetration of toxic compounds, irritants, allergens and microbes1,2. The uppermost layer (15-20 µm thick) of epidermis, known as SC, is mainly responsible for this barrier function.1,2 The multilayer SC is primarily made up of flattened dead cells, known as corneocytes, and intercellular lipid matrix.3 The corneocytes and lipid matrix are packed in brick and mortar assembly, respectively.3 The corneocytes are mostly impermeable as compared to lipid matrix and

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

thus the latter is considered to be the most important pathway for diffusion of substances into the body. The lipid matrix is the major diffusion-rate limiting pathway, as most of the drugs and cosmetics applied topically have to pass through it. This lipid matrix is made up of an approximately equimolar mixture of ceramides, cholesterol (CHOL) and free fatty acid (FFA) with some amount of cholesterol sulphate.3 The high content of CERs, is necessary to maintain the barrier function.3,4 CERs usually contain long saturated tails, unlike phosphatidylglycerols of cell membranes.5 The head group of CERs is rather small as compared to phospholipids and forms an extensive network of hydrogen bonds, which is the basic prerequisite of a competent skin barrier.3,6 The skin barrier function can be mainly attributed to CERs, which is supported by the fact that reduced or altered level of CERs were found in skins that were diseased with atomic dermatitis, ichthyosis and psoriasis.7-12 These lipid barrier abnormalities were corrected to a promising extent by topical application of CERs or their synthetic analogues.13,14 CERs are differentiated on the basis of position, number of hydroxyl (-OH) group and chain length of their tails (both sphingo and fatty acid chain) and based on these almost 15 kinds of subclasses of CERs have been identified.15-17 Out of these CERs, CER-NS is the most abundant and experimentally well studied. J. Novotný et al. synthesized a series of CER analogues with the polar head identical to CER-NS, where sphingosine truncated to 12 carbons, and acyl chain lengths were of 2–12 carbons. It was shown that the skin permeability was unchanged for ceramide analogues with 2C and 12C acyl and CER-NS. Whereas, CER with 4–8C acyl chain increased the permeability of two model permeants up to 10.8 times.18 Barbora et al. reported that long hydrophobic chains in the NS-type CERs are essential for maintaining the skin barrier function, but, this ability was not expressed by their short-chain counterparts.19 It was shown that lower proportion of tight orthorhombic packing in short chain ceramides increased the permeability of the short chain ceramide lipid membranes compared to a long-chain ceramide membrane.19 K. Vávrová et al. showed that simple L-serine based CER analogue 14S24 with the same chain length as the natural CER was able to repair damaged skin barrier14,20, while, its analogue with shorter chains behaved as a permeation enhancer.21 Školová et al. prepared lipid matrix of CER-NS (or its shorter analogue), lignoceric acid, and cholesterol with addition of small amount of cholesterol sulphate and showed that the lipid mixtures arranged themselves into multilayer membrane, and, highest permeability was observed for short ceramide containing membrane.22 Mojumdar et al. also reported that the small chain CER bilayers are more permeable than long chain CER bilayers.23 Despite all these studies, purely experimental approach to design strategies for preventing the water loss, and, excessive skin permeability by using ceramides or its analogue is limited because

ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25

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

of the vast class of ceramides present in the human skin. With the recent molecular level decoding of SC’s structure, it is now possible to develop an appropriate virtual model to accurately mimic its barrier properties.24 Molecular simulations also reveal phenomena occurring at time scale of picoseconds to nanoseconds, which are challenging or unavailable from any of the current experimental techniques. Phospholipids have been studied for the last several decades. These efforts include development of force field parameters25 and study of permeation of different solutes and drug molecules using umbrella sampling and steered MD simulation.26-32 Researchers have started looking at skin lipids from last few decades and significant progress has been made. First simulation of FFA and CHOL mixtures at skin temperature was performed by Holtze et al.33 Pandit and Scott simulated the CER-NS (16:0) bilayer at high temperature of 368 K.34 Effect of permeation enhancer (oleic acid) on model skin membrane was performed by Hoopes at al.35 Notman et al.36 reported that 0.4 mole fraction of Dimethyl sulfoxide (DMSO) transformed the pure CER-NS (24:0) bilayer into liquid phase from gel phase and increased its permeability. First simulation incorporating all three component of SC, CER-NS (24:0), FFA (24:0) and CHOL was performed by Das et al. and the effect of molar ratio of individual components on structural properties were reported.37 They also calculated diffusion and permeability of water through pure CER and mixed bilayer system.38 Recently, Akinshina et al. reported MD simulation study on effect of concentration of monoolein, monostearin, monoelaidin, oleic acid, stearic acid and linoleic acid with CER bilayer.39 It was shown that at low oil concentration (∼5%), the bilayers consisting oils and CER became more rigid than pure CER bilayers due to more efficient lipid packing.39 Gupta et al.40 reported on the phase transition of skin lipids at various compositions and temperature. They also reported the effect of gold nanoparticle size on skin permeability using coarse grained MD simulations.41 Recently, permeation of hydrophilic and hydrophobic though mixed skin bilayer has been reported.42 A theoretical approach, which provides the molecular level insight on the permeability of long and short chain ceramides as well as their effect on structural properties of bilayer, combined with experiments shall provide a rational design framework. In this study, we have presented the molecular level insight of effect of the fatty acid chain length of CER over lateral packing of CER bilayer, and, water permeation through it. The SC is made up of more than 300 different types of CER, which are classified according to the structure of head group and attached fatty acid chain length. For the detailed structures and types of CER present in the SC, the reader may refer to these articles.43,44 Atomistic simulation of molecular model of heterogeneous mixture of lipid matrix of SC layer incorporating all possible CER, FFA and CHOL is challenging and beyond the current computational capabilities. Here we are only using most abundant ceramide, CER-NS for

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

mimicking the SC bilayer, which is a crude simplification of the bilayer. In this study, we have performed several independent, long united atom MD simulations (~500 ns) and steered MD simulations (~25 ns for each configuration) of bilayer of varying ceramide chain length of 8, 12, 16, 18, 20, 22 and 24 Carbon. The permeability of water through each bilayer was calculated using steered MD simulations. The rest of the article is divided into two sections. In the first section, computational details like force field parameters of ceramide, simulation parameters, method of generation of initial configuration, calculation of permeability and diffusion are discussed. The results are discussed in the second section.

1. Computational Details The force field parameters for the CER were taken from the Berger et al.45 and GROMOS87. The CH2 and CH3 groups were treated as single united atom C entity with zero partial charge. The polar hydrogen’s were included explicitly and charges on polar groups were taken from earlier simulations.37,38,40 The Ryckaert-Bellemans dihedral potential was used for the hydrocarbon chains of CER.46 The simple point charge (SPC) model was used for water molecule.47 The simulations were carried out in NVT and NPT ensembles using the latest GROMACS molecular dynamics package.48-51 The temperature was controlled at 310 K by Berendsen thermostat (in equilibration run) with a time constant of 0.5 ps and coupled separately to lipid molecules and water. Later in production run, the thermostat was changed to Nose-Hover with time constant of 0.5 ps. Pressure was controlled at 1 bar by Berendsen (equilibrium run) and Parrinello-Rahman (production run) barostat with a time constant of 5 ps and compressibility of 4.5 x 10-5 with semi isotropic coupling (XY and Z direction coupled separately) to achieve tensionless bilayer. All the bonds of CER were constrained using LINCS algorithm, while SETTLE algorithm was used for water. A time step of 2 fs was used for all the simulations except for pull simulations, where, 1 fs was used. A cut off of 1.2 nm was used for van der Waal and electrostatic interactions. Using a direct cut off for the electrostatic interaction can lead to some artefacts hence in order to avoid these, long range electrostatic interactions were computed using particle mesh Ewald (PME) method. The atomistic simulation for each bilayer was run for 500 ns. The steered simulation for each window (given z position) was run for 25 ns, out of which last 15 ns run was used for the calculation of free energy and permeability.

ACS Paragon Plus Environment

Page 4 of 25

Page 5 of 25

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

The Journal of Physical Chemistry

Figure 1. Molecular structure of ceramide molecules used in simulations. Chains sn1 and sn2 represent lipid tails of CER. Colours red, blue, white and cyan represent oxygen O, nitrogen N, hydrogen H and carbon C (united atom),respectively. Carbon atoms in chain sn2 vary from 8 to 24. Images were created using visual molecular dynamics VMD software.71

Molecule structure of ceramides used in these simulations, are shown in Figure 1. To remove bad contacts between molecules, energy minimization was carried out in vacuum using steepest decent and conjugate gradient method. The minimized CER was then replicated using “genconf” (Gromacs) in X and Y direction to obtain a single layer. This single layer was then replicated in Z direction to obtain the initial bilayer configuration followed by energy minimization. The 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 310 K, by keeping lipid molecules fixed for proper hydration. The system was further run in NPT condition without any constraints. Simulated annealing was finally performed to remove any artefacts related to the structure generation. All structures from NPT run were heated till 370 K and cooled back again to 300 K in 10 ns run as shown in Figure 2. The structures obtained after simulated annealing were subjected to a 500 ns NPT run to arrive at the realistic bilayer structures. Each bilayer comprised of 128 CER and 5120 water molecules.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Temperature (K)

380 360 340 320 300

CER8 CER12 CER16 CER18 CER20 CER22 CER24

0.42 2

apl (nm )

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

0.4 0.38 0.36 0

2000 4000 6000 8000 time (ps)

Figure 2. Simulated annealing run for each bilayer system. Evolution of temperature and area per lipid (aL) with simulation time.

Figure 3. Side view of simulation snapshot of four water molecules constrained at different xy and z position in given window of each simulation system. The water molecules are separated at least 1.2 nm (in Z direction) from each other. Images were created using visual molecular dynamics VMD software.71

For steered MD simulation, the membrane normal (z) was chosen to be the reaction coordinate, where z = 0 nm corresponds to the center of mass (COM) of CER molecules. To sample the whole phase space of the bilayer, multiple molecules need to be constrained at different xy and z plane. We have thus applied the approach in which four water molecules are constrained in a single window, based on the cut off used for electrostatic and van der Waal interactions, as shown in Figure 3. The starting structures for steered simulations were obtained from the final structures of 500 ns NPT run. For each system, four water molecules were inserted manually inside the upper water layer at the distance of ~ 4.8 nm from the COM of the bilayer. Overlapped water molecules

ACS Paragon Plus Environment

Page 6 of 25

Page 7 of 25

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

with constrained water molecules were removed and system was energy minimized. Each system was further equilibrated for 20 ns by keeping constrained molecule fixed at the current position. The solute molecule was then pulled slowly towards the center of the bilayer (-z direction) at the rate of 0.02 nm/ps using pull code of GROMACS. The time step used in the run was 1 fs. As the z distance between the center of mass (COM) of the lipid and water molecules changes by 0.2 nm, the configurations were stored for further simulations. In each window, four water molecules were constrained at different z and xy position as shown in Figure 3. Total 25 windows were generated using the above procedure which spanned the whole phase space from the bulk water (upper leaflet) to the middle of the bilayer. These stored equidistance configurations were ran for another 25 ns. The distance between the COM of water molecules and COM of bilayer was constrained in z direction while the remaining water molecules were free to move in lateral direction. The configurations were stored at every 0.5 ps and constrained force was stored at every 10 fs. First 10 ns run was discarded for equilibration and last 15 ns runs of each simulation were used to calculate diffusion, potential of mean force and permeability of water molecule. To improve the statistical accuracy of sampling, each solute was constrained at four different xy plane for a given z coordinate. Averages of these four simulations for each water molecules are presented here. It was assumed that bilayer is symmetric and result from the one side of bilayer will be the mirror image of the result of other side. The homogenous solubility diffusion model is generally used for the calculation of passive permeability of solutes through membranes. According to this model the solute first dissolves into the membrane, diffuses through the membrane interior, and finally dissolves again in the outer surrounding medium. According to this model, the permeability is given by

=

 

1

where, P is the permeability of the molecules across the bilayer, K is the solute partition function from aqueous to organic phase, D is the diffusion coefficient of solute and d is the thickness of the bilayer. This model has been challenged many times. It has been shown previously that solute size and bilayer composition changes the permeability many folds52 as well as local partition function

K(z) governs the solute permeability.53 Molecular dynamics simulation provides a nice way to calculate both local partition K(z) and local D(z). Marrink et al.26,27 calculated permeability of water molecules through dipalmitoyl phosphatidylcholine (DPPC) bilayer and proposed fourregion model. It was shown that permeability of molecules can be related to their potential of mean force (PMF) across the bilayer. We have adopted the similar strategy to compute the diffusivity, partition function, potential of mean force (PMF) and permeability.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 8 of 25

The permeability for a symmetric bilayer system which has normal in z direction is given by  

/

= /

  



2

where, P is the permeability, d is thickness of the bilayer. K(z) and D(z) is partition function and diffusion coefficient at a given z position from the center of bilayer, which are calculated using following equation

  =

     ∆ , ∆ , 

∆$

  = exp





∆% , & = % , & − < % , & >

∆*  = − ,-./0 < %  + , & >  +

3 4 5 6

where, R is gas constant, T is temperature, F(z,t) is constrained force on solute at a given z and △G is potential of mean force and is averaged force over time.

2. Results and Discussion A. United Atom MD Simulations The snapshots of each bilayer system obtained at the end of 500 ns NPT simulation at 310 K are shown in Figure 4. The mid plane of tail and the mid plane of bilayer correspond to z ~ 1.2-1.4 and z~0 respectively. All the bilayers except CER8 were well ordered in which both sphingo and fatty acid chains of each leaflet were aligned parallel to each other. At physiological temperature (310 K), the lateral packing of CER chains of each of the bilayer was in the gel phase, except CER-8, which had almost liquid like packing. These findings are in good agreement with earlier experimental54-56 and simulation studies.34,37,40 It is interesting to note that in the middle of the bilayer (z ~ 0) the lateral packing density increases with the chain length of the fatty acid. The head groups were found to be packed randomly in each bilayer system. The CER8 bilayer had higher density in the middle of the bilayer as compared to CER12, CER16 and CER18, also, the density in the middle of the bilayer increased with the CER chain length from CER12 onwards. These two factors may explain the experimental phenomena of shorter chain ceramides exhibiting more permeability.18-23

ACS Paragon Plus Environment

Page 9 of 25

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 4. Snapshot of each bilayer system at the end of 500 ns unconstrained simulation. The CER chains, head group oxygen, nitrogen and hydrogen and water are shown in red, blue, white and pink colour, respectively. The mid plane of tail and mid plane of bilayer correspond to z ~ 1.2-1.4 and z~0, respectively. Images were created using visual molecular dynamics VMD software.71 The molecular packing of a lipid bilayer is described by the parameter called area per lipid (aL). It can be compared with the experiments, and it also gives information about the equilibration of the bilayer system during the simulation (please see supplementary information, Figure S1, S2). In a molecular dynamics simulation of a pure lipid bilayer which has the normal along the z direction, the aL can be calculated using the following equation. 23 =

45 46  789:9;

7

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. In similar way, the volume per lipid (vL) could be calculated by subtracting volume of water molecule from the total volume, V = Lx Ly Lz

=> 789:9;

8 9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

where, vL 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.

1.2

0.4

a)

1.1

0.39

3

vL (nm )

2

aL (nm )

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 25

0.38

b)

1 0.9 0.8 0.7

0.37 6 8 10 12 14 16 18 20 22 24 26 Cer chain length

0.6 6 8 10 12 14 16 18 20 22 24 26 Cer chain length

Figure 5. a) Area per lipid (aL) and b) volume per lipid (vL) of each bilayer system calculated in last 300 ns of unconstrained simulation. The aL and vL are shown in Figure 5a and 5b, respectively. The area per lipid was found to be minimum (0.38 nm2) for CER8 and maximum for CER-24 (0.39 nm2). The aL values do not change much for CER12- CER18 and CER20-CER24 bilayers. We have obtained aL value of ~0.38 nm2 for CER16 bilayer. In contrast, Pandit et al.34 reported very high aL value of (0.551±0.007 nm2) for pure CER16 bilayer, however their simulations were carried out at the temperature of 368 K (above the gel to liquid-crystalline phase transition temperature). The aL values of 0.378 nm2 and 0.40 nm2 for CER16 monolayers are reported as determined from surface pressure isotherms experiments.57,58 Notman et al.36 reported aL value of 0.374 nm2 for pure CER24 bilayer at 323 K. Das et al.37 and Gupta et al.40 reported a value of ~0.389 and 0.391 nm2 for pure CER24 bilayer system, respectively. The computed aL for CER24 (~0.39 nm2) at 310K matches very well with x-ray experiment value of ~0.4 nm2 obtained by Dahlen et al.59 as well as those reported in literature.37,40 The volume per lipid increased with the CER chain length. The volume of single SPC water was 0.0305 nm3 at 310 K. The computed vL values and profile are in good agreement with earlier simulation study.60

ACS Paragon Plus Environment

Page 11 of 25

6.5

6

5

a) 5.2 4.99 4.75

4.5

4.47 4.24

4 3.77

3.5 3.07

36

8 10 12 14 16 18 20 22 24

bilayer width d (nm)

O17-O17 C19-C19 O21-O21 N23-N23

5.5 Distances (nm)

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

6

b)

5.5 5 4.5 4 3.5 8 10 12 14 18 20 22 24 6 16 26

CER Chain length

CER chain length

Figure 6. a) Distance between the head group atoms and b) bilayer width of each bilayer system calculated in last 300 ns of unconstrained simulation. The atom numbers are marked in the right most image. The image was created using visual molecular dynamics VMD software.71 The average thickness of a bilayer is calculated based on volume and area per lipid using the following equation.

=

 ? @4

10

A4

Figure 6 shows the distance between the head group atoms of both of the leaflet and bilayer width of each bilayer system. The atom names are also shown in the Figure 6. The distance and bilayer width increased with CER chain length. Tail order parameter is used to characterize the ordering in lipid bilayers and can be measured by deuterium NMR. An order parameter for every CH2 group in the chains is defined as C







B = DEF  G −

11

where, θ 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. The overall order parameter was calculated using the following relationship: B=

∑J 9KL I /

12

M

where, n is number of CH2 in the ceramide molecules and Sz is order parameter for ith CH2 of ceramide chain as shown in the Figure 1.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 a)

1

sn1 chain

sn2 chain

0.9

b)

0.8

0.8 0.6

CER8 CER12 CER16 CER18 CER20 CER22 CER24

0.4 0.2 0 2 4 6 8 10 12 14 Atom number

4

8

12 16 20

Atom number

S

0.7 Sz

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 25

0.6 0.5 0.4 0.3 0.2 6 8 10 12 14 16 18 20 22 24 26 Cer chain length

Figure 7. a) Order parameter (Sz) of CER chains sn1 and sn2. b) Overall order parameter S of CER chain in each bilayer system. The chain sn1 and sn2 are shown in Figure 1. The individual order parameter and overall order parameter are shown in Figure 7a and 7b, respectively. The order parameter profile was similar in each bilayer system. CER bilayers with larger chains form highly ordered structures and the ordering decreased with the chain length. The CER-8 exhibited had very low order parameter which confirms that bilayer is in liquid crystalline phase and this could be one of the possible reasons for the higher permeability of shorter chain CER bilayer. Experimentally, FTIR spectra of pure CER bilayer have also shown the decreased ordering in CER8 and CER12.19 It is interesting to note that the when the chain length of fatty acid is equal or greater than sphingo chain (16), the ordering increased for sn1 chain significantly. The overall order parameter shows a bell curve kind of behaviour. It increased with the chain length, reached a maximum for CER16 and then decreased. The possible explanation could be that due to the similar lengths of fatty acid and sphingo chain length, they straighten each other, instead of interdigitizing with the other leaflet chains. In Figure 8 the electron density profile of each bilayer along the bilayer normal is plotted. The profiles are almost similar for each bilayer. As the chain length increases the electron density in the mid of the bilayer plane increases, except for CER-8. For CER-8 bilayer the density in the middle of the bilayer was greater than that of other ceramide bilayer. The dissimilar chain length of CER24 interdigitate with the other leaflet lipids and increase the density in the middle of the bilayer. The difference in the peak positions of the electron density also gives information of the bilayer thickness. The peak positions increased with the CER chain length.

ACS Paragon Plus Environment

Page 13 of 25

Partial density

-3

Electron density (e nm )

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

500 450 400 350 300 250 200 150 100

CER8 CER12 CER16 CER18 CER20 CER22 CER24

-4

-2

0

z (nm)

2

4

Figure 8. Electron density profile along the bilayer normal z in each bilayer system calculated in last 300 ns unconstrained simulation.

B. Steered MD Simulation : Water Permeability The steered simulations were performed to obtain the transport and thermodynamics properties of water permeation through each bilayer system. The Figure 9a shows the free energy of permeation of water through each bilayer at 310 K, calculated using equation 6. The errors in the force calculation are standard errors calculated from the difference of the mean force 〈F(z)〉 in each of individual simulations at each depth from their average. The errors in △G are calculated by propagating the errors in the force F(z). The free energy follows the similar trends for each bilayer system. For each bilayer the free energy was almost zero in the bulk water phase, and as constrained water reach near the head group of the bilayer, the free energy increased. The △G increased due to the loss of favourable electrostatic interactions and hydrogen bonds when leaving from the bulk water phase to the hydrophobic interior of the bilayer. After crossing the head group region, free energy increased till middle of the upper leaflet of bilayer and reached a maximum. The reason could be the highly ordered hydrophobic chains of lipid bilayers. Moving towards the center of the bilayer it decreased but it was more than that of free energy in bulk water phase.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Delta G(z) kJ/mol

40 30 20

CER8 CER12 CER16 CER18 CER20 CER22 CER24

b) 5e-05 4e-05

2

CER8 CER12 CER16 CER18 CER20 CER22 CER24

a)

D(z) cm /s

50

10

3e-05 2e-05 1e-05

0

-4

-2

0

z (nm)

2

4

-4

0

10 -1 10 c) -2 10 -3 10 -4 10 -5 10 -6 10 -7 10 -8 10 -9 10 -4

-2

0

z (nm)

2

CER8 CER12 CER16 CER18 CER20 CER22 CER24

-2

0

z (nm)

2

4

10 12 10 11 10 10 10 9 10 8 10 7 10 6 10 5 10 4 10

d)

-4

4

CER8 CER12 CER16 CER18 CER20 CER22 CER24

13

R(z)

K(z)

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 25

-2

0

z (nm)

2

4

Figure 9. a) Local free energy of permeation △G(z) b) local diffusion coefficient D(z) c) local partition coefficient K(z) and d) resistance of permeation R(z) of water molecule along the bilayer normal z calculated from multiple water molecules in single window constrained simulation. The bilayer was assumed to symmetric and results from one leaflet were replicated for another leaflet. The bilayer center correspond to z=0. For the purpose of clarity, the error bars are not shown here, please refer the supplementary information for the same (Figure S7). One possible explanation could be the lesser density of chains in the middle of the bilayer (please see supplementary information, Figure S5). The free energy profile of CER24 has the similar shape as obtained by the Das et al. at 300 K and 350 K for the water permeation from the pure CER and model skin membrane.38 The △G(z) increased smoothly from head group region and reached to maximum (~ 47.6 kJ/mol) in hydrophobic core and then decreased (~34 kJ/mol) at the center of the bilayer. Das et al.38 reported a maximum free energy of 42.7 kJ/mol and ~28 kJ/mol at the center of the pure CER24 bilayer at 350 K. Our values are higher than these values since the temperature was 310 K in our simulations as compared to Das et al. simulations (350K).38 The

ACS Paragon Plus Environment

Page 15 of 25

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

maximum values of △G(z) (~22.9-26 kJ/mol) for the water permeation through DPPC bilayer has been reported earlier.26,29 In our simulation we obtained almost 2 times more free energy as compared to these values. The permeability decreases exponentially with free energy (eq. 2) so the obtained values clarify that why skin provides high barrier for water permeation as compared to biological membrane. Also, Paloncýová et al.61 carried out MD simulations of p-aminobenzoic acid (PABA) and its ethyl (benzocaine) and butyl (butamben) esters with gel phase CER and fluid state dioleoylphosphatidylcholine (DOPC) bilayer and reported that these drug molecules exhibited penetration barriers almost 3 times higher in the gel phase CER24 bilayer. The lower permeability of the drugs were explained in terms of elevated barrier and slower diffusion in the CER24 bilayer, which were caused by the high ordering of CER2 lipid chains.61 The free energy profile of water across the CER8 bilayer were found to be similar to water free energy across the fluid phase DPPC bilayer at 323 K.62 Our profiles are well aligned with the earlier MD simulation study of DPPC bilayer of different chain length of 12,14 and 16 and showed that the free energy barrier increased with the chain length.63 The local density of the CER24 bilayer interior was found to be around ~750 kg/m3 which is comparable to density of hexadecane of 765.31 kg/m3 at 308 K.64 The free energy of solvation of water in liquid hexadecane at 310 K is ~20 kJ/mol65 which is lesser than the bilayer mid-plane values of ~ 34 kJ/mol. The observed higher free energy in the bilayer than in the hexadecane may be due to restricted freedom of water to move in the interior of ordered bilayer as compared to no such hindrance in the bulk hexadecane solution. Recent MD simulation study of permeation of hydrophobic molecules through skin lipids reported similar phenomena.42 Experimentally, it has been shown that small chain ceramides are more permeable than longer chain ceramides.18-23 This phenomenon can be explained on the basis of change in free energy of permeation, which is the difference between the peak value and the value at the centre of bilayer. Lower values indicate that the bilayer is relatively more permeable. This shows that, CER8 is more permeable than any of the higher chain ceramides. Diffusion coefficient D(z) profiles for water in each bilayer, calculated using equation 3, are plotted in Figure 9b. The errors in the diffusivity calculation are standard errors calculated from the difference of the diffusion coefficient in each of individual simulations at each z position from their average. The shape of D(z) profile for each bilayer was found to be similar in nature. The D(z) values in bulk water was found to be greater than that of bilayer interior because of higher density of the lipid environment inside the bilayer. As the water molecule moves down towards the bilayer from the water phase, the D(z) value decreases first and reaches a minimum value and then it increases till the middle of the bilayer. This may be due to less density at the middle of the

ACS Paragon Plus Environment

The Journal of Physical Chemistry

bilayer as shown in Figure S5 (please see supplementary information). It is interesting to note that the sudden drop in the D(z) for each bilayer happened at different z position. This is because of the smaller thickness of the short chain CER bilayers. The D(z) values for water in bulk is ~4.5 x 10-5 cm2/s, which is comparable to the value (~5.0 x 10-5 cm2/s) obtained by Das et al. at 300K.38 The experimental value of water self-diffusion coefficient is ~3.0x10-5 cm2/s.66 SPC water model predicts diffusivity almost ~1.3 times more than that of experimental values.26 Thus our diffusion values are in good agreement with previous experimental and computational studies. -2

10

10

-5

-3

10 P (cm/s)

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 25

-4

10

10

-6

15

20

25

-5

10

-6

10

6 8 10 12 14 16 18 20 22 24 26 Cer Chain Length

Figure 10. Permeability of water molecule in each ceramide bilayer calculated using multiple water molecules constrained in single window constrained simulation. Figure 9c shows the local partition coefficient of water along the bilayer normal calculated using the equation 4. The shape of K(z) profile is complementary to that of free energy profile. The K(z) was found to be maximum in the bulk water and minimum in bilayer interior for each bilayer system. Resistance to permeation R(z) profiles for each water molecule, calculated using R(z) = 1/K(z)D(z), are plotted in Figure 9d. The error bars are not shown for the purposes of clarity. The error bars are higher due to exponential dependence of K(z) or △G (z), even a small change in K(z) amplified the values of R(z). The R(z) profiles follows the similar pattern as △G(z), because of exponential relationship between them. The main resistance of permeation of water was observed in the interior of the bilayer because of less favourable solute partition and diffusion coefficient as compared to water phase. The permeation resistance mostly depends upon the free energy values and the effect of diffusion coefficients looks almost negligible. Bempoard et al.29 obtained values of R(z) of water in DPPC bilayer in the order of ~108, whereas we have obtained in the order of

ACS Paragon Plus Environment

Page 17 of 25

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

~1010-1013 for water, which shows the less permeability of skin membrane as compared to biological membrane. The permeability values (P) calculated using integration of R(z) along the bilayer normal, are plotted in Figure 10. The water permeability for CER24 bilayer was found to be (1.12 ± 0.78) x 10-6 cm/s at 310 K. The experiments performed on porcine and human stratum-corneum were also found to give P values in order of ~ 10-7 cm/s.67, 68 Bemporad et al.29 also obtained an order higher water permeability for DPPC bilayer at 350 K, as compared to experiment. It is worth noting that the human skin has almost ~300 different types of lipids and chain length as well as head group of these lipids change the lamellar and lateral organization of lipid layer which further effects the barrier function of skin significantly.18-23 The permeability for the water was found to be mostly dependent on the free energy. Because of exponential dependence, the observed permeability has order of difference with respect to chain length. The smaller chain length CER8 bilayer had very high permeability as compared to that of CER24 bilayer. The CER8 bilayer was in liquid crystalline phase and the water permeability order was in agreement with water permeability for DPPC bilayer obtained at 350 K above the gel to liquid phase transition temperature.29

Conclusion In this article, long molecular dynamics simulations of pure ceramide- NS bilayers of different fatty acid chain length (from 8 to 24) were performed to calculate the equilibrium and the water transport properties through these bilayer. The inhomogeneous solubility-diffusion model used in this study incorporates both the equilibrium (free energy) and the dynamic (diffusion coefficient) behavior of water in the bilayer. The main findings are that the lipid packing is significantly affected by the CER chain length. The inter leaflet interdigitization between the CER chain were absent in the smaller chain CER, because of that CER8 bilayer was in disordered and almost like liquid crystalline phase. Moreover, the area per lipid has not changed much with the CER chain length while volume per lipid and bilayer thickness increased with the CER chain length. The main resistance of permeation for water was obtained in the interior of the bilayer for each system. The permeability of water was found to be mostly dependent upon the free energy. The obtained permeability of water through skin lipid membrane is in qualitative agreement with experiments. There could be many possible reasons for the same. First, the skin model used here only comprises most abundant single kind of CER (CER-NS) but in reality skin lipid matrix is very heterogeneous. The lipid matrix of the skin is the key determinant for the barrier functions and is mainly composed of heterogeneous mixture of long chain ceramides (CER), cholesterol (CHOL) and free

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 18 of 25

fatty acid (FFA) in certain ratios.43,44 The CHOL and FFA changes arrangement and packing of CER molecules significantly,37,40 which could change the permeability significantly.38 Second, slight changes in the head group of CER (either position or number of OH group) can change the morphology of bilayer significantly.69 Recently, Uchiyama et al. reported that even fatty acid chain distribution also plays important role in order to determine the skin permeability.70 There are 11 different kinds of CERs present in human skin, but, atomistic simulation by incorporating all kind of skin lipids is beyond the current computational capabilities. Coarse-grained simulation could be attempted to simulate a full skin lipid bilayer. Recently, a coarse-grained MD simulation studies on skin lipids have been reported.41 Currently, we are developing coarse grained model of other CERs and FFAs, results of which will be communicated in separate article.

Supporting Information. S1. Evolution of Potential energy; S2. Evolution of area and volume per lipid; S3. Structure of each bilayer in the end of 500 ns unconstrained simulation; S4. Evolution of distance between head group atoms; S5. Density profiles of each components of bilayer along the bilayer normal z.; S6. Lateral diffusion coefficient of ceramide molecule in each bilayer system; S7. Free energy of permeation and diffusion along the bilayer normal z; S8. Bilayer width based on electron density profile; S9. Effect of system size on structural and transport properties of bilayer.

Acknowledgement Authors would like to thank High Performance Computing Business Unit at Tata Consultancy Services (TCS) for providing access to EKA Super computer. They would also like to thank Mr. K Ananth Krishanan, CTO, TCS, for his constant encouragement and support during this project.

References 1. Menon, G. K.; Cleary, G. W.; Lane, M. E. The Structure and Function of the Stratum Corneum. Int. J. Pharm. 2012, 435, 3–9. 2. Elias, P. M. Epidermal Lipids, Barrier Function, and Desquamation. J. Invest. Dermatol.

1983, 80, 44s-49s. 3. Wertz P.W. Lipids and Barrier Function of the Skin. Acta Dermato VenereologicaSupplement. 2000, 208, 7-11. 4. Wertz P.W. Current Understanding of Skin Biology Pertinent to Skin Penetration: Skin Biochemistry. Skin Pharmacol Physiol. 2013, 26, 217-226.

ACS Paragon Plus Environment

Page 19 of 25

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

5. Van M. G.; Voelker D.R.; Feigenson G.W. Membrane Lipids: Where They are and How They Behave. Nat. Rev. Mol. Cell Biol. 2008, 9, 112-124. 6. Pascher I. Molecular Arrangements in Sphingolipids Conformation and Hydrogen Bonding of Ceramide and Their Implication on Membrane Stability and Permeability. BBA Biomem.

1976, 455, 433-451. 7. Bouwstra J.A.; Ponec M. The Skin Barrier in Healthy and Diseased State. BBA Biomem.

2006, 1758, 2080-2095. 8. Choi M.J.; Maibach H.I. Role of Ceramides in Barrier Function of Healthy and Diseased Skin. Am. J. Clin. Dermatol. 2005, 6, 215-223. 9. Janssens M.; van Smeden J.; Gooris G.S.; Bras W.; Portale G.; Caspers P.J.; Vreeken R.J.; Kezic S.; Lavrijsen A.P.; Bouwstra J.A. Lamellar Lipid Organization and Ceramide Composition in the Stratum Corneum of Patients with Atopic Eczema. J. Invest. Dermatol.

2011, 131, 2136-2138. 10. Janssens M.; van Smeden J.; Gooris G.S.; Bras W.; Portale G.; Caspers P.J.; Vreeken R.J.; Kezic S.; Lavrijsen A.P. Increase in Short-Chain Ceramides Correlates with an Altered Lipid organization and Decreased Barrier Function in Atopic Eczema Patients. J. Lipid Res. 2012, 53, 2755-2766. 11. Ishikawa, J.; Narita, H.; Kondo, N.; Hotta, M.; Takagi, Y.; Masukawa, Y.; Kitahara, T.; Takema, Y.; Koyano, S.; Yamazaki, S. et al. Changes in the Ceramide Profile of Atopic Dermatitis Patients. J. Invest. Dermatol. 2010, 130, 2511-2514. 12. Eckl, K.; Tidhar, R.; Thiele, H.; Oji, V.; Hausser, I.; Brodesser, S.; Preil, M.; Önal-Akan, A.; Stock, F.; Müller, D. et al. Impaired Epidermal Ceramide Synthesis Causes Autosomal Recessive Congenital Ichthyosis and Reveals the Importance of Ceramide Acyl Chain Length. J. Invest. Dermatol. 2013, 133, 2202-2211. 13. Chamlin, S.; Kao, J.; Frieden, I.; Sheu, M.; Fowler, A.; Fluhr, J.; Williams, M.; Elias, P. Ceramide-Dominant Barrier Repair Lipids Alleviate Childhood Atopic Dermatitis: Changes in Barrier Function Provide a Sensitive Indicator of Disease Activity. J. Am. Acad. Dermatol. 2002, 47, 198-208. 14. Vávrová, K.; Hrabálek, A.; Mac-Mary, S.; Humbert, P.; Muret, P. Ceramide Analogue 14S24 Selectively Recovers Perturbed Human Skin Barrier. Br. J. Dermatol. 2007, 157, 704-712. 15. Masukawa, Y.; Narita, H.; Shimizu, E.; Kondo, N.; Sugai, Y.; Oba, T.; Homma, R.; Ishikawa, J.; Takagi, Y.; Kitahara, T. et al. Characterization of Overall Ceramide Species in Human Stratum Corneum. J. Lipid Res. 2008, 49, 1466-1476.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 20 of 25

16. van Smeden, J.; Hoppel, L.; van der Heijden, R.; Hankemeier, T.; Vreeken, R.; Bouwstra, J. LC/MS Analysis of Stratum Corneum Lipids: Ceramide Profiling and Discovery. J. Lipid Res. 2011, 52, 1211-1221. 17. T’Kindt, R.; Jorge, L.; Dumont, E.; Couturon, P.; David, F.; Sandra, P.; Sandra, K.. Profiling

and

Characterizing

Skin

Ceramides

Using

Reversed-Phase

Liquid

Chromatography–Quadrupole Time-Of-Flight Mass Spectrometry. Anal. Chem. 2012, 84, 403-411. 18. Novotný, J.; Janůšová, B.; Novotný, M.; Hrabálek, A.; Vávrová, K. Short-Chain Ceramides Decrease Skin Barrier Properties. Skin Pharmacol. Physiol. 2009, 22, 22–30. 19. Janůšová, B.; Zbytovská, J.; Lorenc, P.; Vavrysová, H.; Palát, K.; Hrabálek, A.; Vávrová, K. Effect of Ceramide Acyl Chain Length on Skin Permeability and Thermotropic Phase Behavior of Model Stratum Corneum Lipid Membranes. BBA Mol. Cell. Biol. L. 2011, 1811, 129–137. 20. Vávrová, K.; Zbytovská, J.; Palát, K.; Holas, T.; Klimentová, J.; Hrabálek, A.; Doležal, P. Ceramide

Analogue

14S24

((S)-2-Tetracosanoylamino-3-Hydroxypropionic

Acid

Tetradecyl Ester) is Effective in Skin Barrier Repair In Vitro. Eur. J. Pharm. Sci. 2004, 21, 581-587. 21. Vávrová, K.; Hrabálek, A.; Doležal, P.; Holas, T.; Zbytovská, J. L-Serine and Glycine Based Ceramide Analogues as Transdermal Permeation Enhancers: Polar Head Size and Hydrogen Bonding. Bioorg. Med. Chem. Lett. 2003, 13, 2351-2353. 22. Školová, B.; Janůšová, B.; Zbytovská, J.; Gooris, G.; Bouwstra, J.; Slepička, P.; Berka, P.; Roh, J.; Palát, K.; Hrabálek, A. et al. Ceramides in the Skin Lipid Membranes: Length Matters. Langmuir. 2013, 29, 15624-15633. 23. Mojumdar, E.; Kariman, Z.; van Kerckhove, L.; Gooris, G.; Bouwstra, J. The Role of Ceramide Chain Length Distribution on the Barrier Properties of the Skin Lipid Membranes. BBA Biomem. 2014, 1838, 2473-2483. 24. Iwai, I.; Han, H.; Hollander, L.; Svensson, S.; Öfverstedt, L.; Anwar, J.; Brewer, J.; Bloksgaard, M.; Laloeuf, A.; Nosek, D. et al. The Human Skin Barrier is Organized as Stacked Bilayers of Fully Extended Ceramides with Cholesterol Molecules Associated with the Ceramide Sphingoid Moiety. J. Invest. Dermatol. 2012, 132, 2215-2225. 25. Egberts, E.; Marrink, S.; Berendsen, H. Molecular Dynamics Simulation of a Phospholipid Membrane. Eur. Biophys. J. 1994, 22, 423-436. 26. Marrink, S.Berendsen, H. Simulation of Water Transport through a Lipid Membrane. Phys. Chem. 1994, 98, 4155-4168.

ACS Paragon Plus Environment

J.

Page 21 of 25

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

27. Marrink, S.Berendsen, H. Permeation Process of Small Molecules across Lipid Membranes Studied by Molecular Dynamics Simulations. J. Phys. Chem. 1996, 100, 16729-16738. 28. Tieleman, D.; Marrink, S.; Berendsen, H. A Computer Perspective of Membranes: Molecular Dynamics Studies of Lipid Bilayer Systems. BBA Rev Biomem. 1997, 1331, 235-270. 29. Bemporad, D.; Essex, J.; Luttmann, C. Permeation of Small Molecules through a Lipid Bilayer: A Computer Simulation Study. J. Phys. Chem. B 2004, 108, 4875-4884. 30. Bemporad, D.; Luttmann, C.; Essex, J. Behaviour of Small Solutes and Large Drugs in a Lipid Bilayer from Computer Simulations. BBA Biomem. 2005, 1718, 1-21. 31. Nademi, Y., Amjad Iranagh, S., Yousefpour, A., Mousavi, S. and Modarress, H. Molecular Dynamics Simulations and Free Energy Profile of Paracetamol in DPPC and DMPC Lipid Bilayers. J Chem Sci. 2014, 126, 637-647. 32. Orsi, M.Essex, J. Permeability of Drugs and Hormones through a Lipid Bilayer: Insights from Dual-Resolution Molecular Dynamics. Soft Matter. 2010, 6, 3797-3808. 33. Höltje, M.; Förster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Höltje, H. Molecular Dynamics Simulations of Stratum Corneum Lipid Models: Fatty Acids and Cholesterol. BBA Biomem. 2001, 1511, 156-167. 34. Pandit, S.Scott, H. Molecular Dynamics Simulation of a Ceramide Bilayer. J. Chem. Phys. 2006, 124, 014708. 35. Hoopes, M.; Noro, M.; Longo, M.; Faller, R. Bilayer Structure and Lipid Dynamics in a Model Stratum Corneum with Oleic Acid. J. Phys. Chem. B 2011, 115, 3164-3171. 36. Notman, R.; den Otter, W.; Noro, M.; Briels, W.; Anwar, J. The Permeability Enhancing Mechanism of DMSO in Ceramide Bilayers Simulated by Molecular Dynamics. Biophys. J. 2007, 93, 2056-2068. 37. Das, C.; Noro, M.; Olmsted, P. Simulation Studies of Stratum Corneum Lipid Mixtures. Biophys. J. 2009, 97, 1941-1951. 38. Das, C.; Olmsted, P.; Noro, M. Water Permeation through Stratum Corneum Lipid Bilayers from Atomistic Simulations. Soft Matter. 2009, 5, 4549-4555. 39. Akinshina, A.; Das, C.; Noro, M. Effect of Monoglycerides and Fatty Acids on a Ceramide Bilayer. Phys. Chem. Chem. Phys. 2016, 18, 17446-17460. 40. Gupta, R.; Rai, B. Molecular Dynamics Simulation Study of Skin Lipids: Effects of the Molar Ratio of Individual Components over a Wide Temperature Range. J. Phys. Chem. B

2015, 119, 11643-11655.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 22 of 25

41. Gupta R, Rai B. Penetration of Gold Nanoparticle through Human Skin: Unraveling Mechanisms at Molecular Scale. J. Phys. Chem. B 2016, 120, 7133–7142. 42. Gupta, R.; Sridhar, D.; Rai, B. Molecular Dynamics Simulation Study of Permeation of Molecules through Skin Lipid Bilayer. J. Phys. Chem. B 2016, 120, 8987-8996. 43. Bartosova, L.Bajgar, J. Transdermal Drug Delivery In Vitro Using Diffusion Cells. CMC. 2012, 19, 4671-4677. 44. Das, C.; Noro, M.; Olmsted, P. Lamellar and Inverse Micellar Structures of Skin Lipids: Effect of Templating. Phys. Rev. Lett. 2013, 111,148101. 45. Berger, O.; Edholm, O.; Jähnig, F. Molecular Dynamics Simulations of a Fluid Bilayer of Dipalmitoylphosphatidylcholine at Full Hydration, Constant Pressure, and Constant Temperature. Biophys. J. 1997, 72, 2002-2013. 46. Ryckaert, J.Bellemans. A. Molecular Dynamics of Liquid N-Butane near its Boiling Point. Chem. Phys. Lett. 1975, 30, 123-125. 47. Berendsen, H. J.; Postma, J. P.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces 1981, 331-342. 48. Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Com. Phys. Comm. 1995, 91, 43-56. 49. 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. 50. Kutzner, C.; Van Der Spoel, D.; Fechner, M.; Lindahl, E.; Schmitt, U. W.; De Groot, B. L.; Grubmüller, H. Speeding up parallel GROMACS on High‐Latency Networks. J. Comput. Chem. 2007, 28, 2075-2084. 51. Pronk, S.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics. 2013, 29, 845−54. 52. Yang, Y.; Sunoqrot, S.; Stowell, C.; Ji, J.; Lee, C.; Kim, J.; Khan, S.; Hong, S. Effect of Size, Surface Charge, and Hydrophobicity of Poly(Amidoamine) Dendrimers on their Skin Penetration. Biomacromolecules. 2012, 13, 2154-2162. 53. Zocher, F.; van der Spoel, D.; Pohl, P.; Hub, J. Local Partition Coefficients Govern Solute Permeability of Cholesterol-Containing Membranes. Biophys. J. 2013, 105, 2760-2770. 54. Moore, D.; Rerek, M.; Mendelsohn, R. FTIR Spectroscopy Studies of the Conformational Order and Phase Behavior of Ceramides. J. Phys. Chem. B 1997, 101, 8933-8940.

ACS Paragon Plus Environment

Page 23 of 25

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. Shah, J.; Atienza, J. M.; Duclos, R. I.; Rawlings, A.V.; Dong, Z. H.; Shipley, G. G. Structural and Thermotropic Properties of Synthetic C16: 0 (Palmitoyl) Ceramide: Effect of Hydration. J. Lipid Res. 1995, 36, 1936-1944. 56. 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. 57. Brockman, H.; Momsen, M.; Brown, R.; He, L.; Chun, J.; Byun, H.; Bittman, R. The 4,5Double Bond of Ceramide Regulates its Dipole Potential, Elastic Properties, and Packing Behavior. Biophys. J. 2004, 87, 1722-1731. 58. Scheffer, L.; Solomonov, I.; Jan Weygand, M.; Kjaer, K.; Leiserowitz, L.; Addadi, L. Structure of Cholesterol/Ceramide Monolayer Mixtures: Implications to the Molecular Organization of Lipid Rafts. Biophys. J. 2005, 88, 3381-3391. 59. Dahlén, B.Pascher, I. Molecular Arrangements in Sphingolipids. Thermotropic Phase Behaviour of Tetracosanoylphytosphingosine. Chem. Phys. Lipids. 1979, 24, 119-133. 60. Paloncýová, M.; Vávrová, K.; Sovová, Ž.; DeVane, R.; Otyepka, M.; Berka, K. Structural Changes in Ceramide Bilayers Rationalize Increased Permeation through Stratum Corneum Models with Shorter Acyl Tails. J. Phys. Chem. B 2015, 119, 9811-9819. 61. Paloncýová, M.; DeVane, R.; Murch, B.; Berka, K.; Otyepka, M. Rationalization of Reduced

Penetration

of

Drugs

through

Ceramide

Gel

Phase

Membrane. Langmuir. 2014, 30, 13942-13948. 62. Qiao, B.Olvera de la Cruz, M. Driving Force for Water Permeation across Lipid Membranes. J. Phys. Chem. Lett. 2013, 4, 3233-3237. 63. Sugii, T.; Takagi, S.; Matsumoto, Y. A Molecular-Dynamics Study of Lipid Bilayers: Effects of the Hydrocarbon Chain Length on Permeability. J. Chem. Phys. 2005, 123, 184714. 64. Domańska, U.; Łachwa, J.; Morawski, P.; Malanowski, S. Phase Equilibria and Volumetric Properties in Binary Mixtures Containing Branched Chain Ethers (Methyl 1,1Dimethylethyl Ether Or Ethyl 1,1-Dimethylethyl Ether Or Methyl 1,1-Dimethylpropyl Ether Or Ethyl 1,1-Dimethylpropyl Ether). J. Chem. Eng. Data. 1999, 44, 974-984. 65. Wu, Z.; Cui, Q.; Yethiraj, A. A New Coarse-Grained Force Field for Membrane–Peptide Simulations. J. Chem. Theory Comput. 2011, 7, 3793-3802. 66. Holz, M.; Heil, S.; Sacco, A. Temperature-Dependent Self-Diffusion Coefficients of Water and Six Selected Molecular Liquids for Calibration in Accurate 1H NMR PFG Measurements. Phys. Chem. Chem. Phys. 2000, 2, 4740-4742.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 24 of 25

67. Potts, R.Francoeur, M. The Influence of Stratum Corneum Morphology on Water Permeability. J. Invest. Dermatol. 1991, 96, 495-499. 68. Blank, I.; Moloney, J.; Emslie, A.; Simon, I.; Apt, C. The Diffusion of Water across the Stratum Corneum as a Function of its Water Content. J. Invest. Dermatol. 1984, 82, 188194. 69. Guo, S.; Moore, T.; Iacovella, C.; Strickland, 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, 5116-5126. 70. Uchiyama, M.; Oguri, M.; Mojumdar, E.; Gooris, G.; Bouwstra, J. Free Fatty Acids Chain Length Distribution Affects the Permeability of Skin Lipid Model Membranes. BBA Biomem. 2016, 1858, 2050-2059. 71. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph.

1996, 14, 33-38.

ACS Paragon Plus Environment

Page 25 of 25

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

TOC

ACS Paragon Plus Environment