Molecular Dynamics Simulation Study of Permeation of Molecules

Aug 12, 2016 - Physical Science Research Area, TCS Research, Tata Research Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar ...
0 downloads 0 Views 1MB Size
Subscriber access provided by Northern Illinois University

Article

Molecular Dynamics Simulation Study of Permeation of Molecules Thorough Skin Lipid Bilayers Rakesh Gupta, Balarama Sridhar Dwadasi, and Beena Rai J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b05451 • Publication Date (Web): 12 Aug 2016 Downloaded from http://pubs.acs.org on August 15, 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 26

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

Permeation of Molecules through Skin Lipid Bilayers. Rakesh Gupta, D B Sridhar and Beena Rai* Physical Science Research Area, TCS Research Tata Research Development and Design Centre, Tata Consultancy Services, 54B, Hadapsar Industrial Estate, Pune – 411013, INDIA *Corresponding author: [email protected] Fax:

91-20-66086399

Tel:

91-20-66086203

Abstract. Stratum Corneum (SC), the outermost layer of skin, is mainly responsible for skin’s barrier function. The complex lipid matrix of SC determines these barrier properties. In this study, the lipid matrix is modeled as equimolar mixture of ceramide (CER), cholesterol (CHOL) and free fatty acid (FFA). The permeation of water, oxygen, ethanol, acetic acid, urea, butanol, benzene, dimethyl sulfoxide (DMSO), toluene, phenol, styrene and ethyl benzene across this layer is studied using constrained MD simulations technique. Several long constrained simulations are performed at skin temperature of 310 K in NPT conditions. The free energy profiles and diffusion coefficients along the bilayer normal have been calculated for each molecule. Permeability coefficients are also calculated and compared with experimental data. The main resistance for the permeation of hydrophilic and hydrophobic permeants has been found to be in the interior of the lipid bilayer and near the lipid-water interface, respectively. The obtained permeability is found to be a few orders of magnitude higher than experimental values for hydrophilic molecules while for hydrophobic more discrepancy was observed. Overall, the qualitative ranking is consistent with the experiments.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

1. Introduction Skin, the largest organ of our body, protects us from the attack of foreign pathogens, provides barrier to the permeation of many harmful molecules and maintains the hydration level of tissues. The outer layer of skin, also known as Stratum Corneum (SC) is mainly responsible for these barrier properties.1 The delivery of drugs through the skin provides a convenient route of administration because of high surface area of skin and can typically be self-administered.1,2 The accurate prediction of dermal uptake of chemicals is therefore relevant to both transdermal drug delivery as well as topical application of cosmetics. In recent past, transdermal delivery of small molecules (hydrophobic/hydrophilic drugs) and macromolecules (proteins, enzymes) has been an attractive as well as a challenging area of research.3-6 The transdermal route provides following advantages over conventional oral/intravenous routes. •

Elimination of first pass metabolism allows lesser amount of drug to be administered thereby reducing its toxicity.



Drugs with poor or less bio-availability can be delivered due to avoidance of the GI tract (a major barrier in oral route).



Pain free and ease of administration (can be self-administered).



Control drug release over a long period thereby eliminating frequent dosing.



Enhanced patience compliance.

The current industry standard, however, both in pharma and cosmetics, is to conduct detailed invitro and in-vivo trials. These obviously incur huge expenses thereby leading to a very few successful candidates that are finally approved by regulatory authority (FDA). The 2-D in vitro cell culture studies do not accurately reflect the complex interactions that occur between the multiple cells present in the 3-D in vivo skin environment. In vivo studies in rodents and other small animals do not translate well to the human situation due to differences in anatomical structures. Though some of the commercially available human skin equivalents EpiSkin® (L’Oreal, Paris) and EpiDermTM (MatTek, Massachusetts) are available in the market, these require highly specialized skills and are very expensive. The European Union (EU) regulation (76/768/EEC, Feb. 2003) prohibits the use of animal or animal-derived substances for the development and testing of cosmetic and pharmaceutical ingredients

ACS Paragon Plus Environment

Page 2 of 26

Page 3 of 26

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 realistic in-silico model of human skin does not exist. Considering the time and costs involved in the development and testing of new drug/cosmetics formulations, it is imperative to supplement/replace some of the elaborate in-vivo/in-vitro tests with in-silico tests. Extensive research has been carried out over many years to predict skin permeability of various drugs and cosmetic molecules. These efforts include experimental studies5,7-9, development of theoretical models3 and empirical methods.10 A molecular-level understanding of SC, which can ultimately lead to the development of rapid in-silico screening tools to predict permeability of a molecule from knowledge of its molecular structure alone, is still not available. 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.11 The solute transport across the SC layer occurs primarily by passive diffusion. The SC is selectively permeable and allows only relatively lipophilic compounds to diffuse into the layers. Therefor breaching of skin barrier for delivery of molecules is necessary and it can be done in many ways and is broadly classified into two categories namely, passive methods such as chemical penetration enhancers and active methods such as electroporation, iontophoresis, sonophoresis and thermophorasis.12 Molecular simulation provides a convenient way to understand these processes and yield important physical insights.13 The SC mainly consists of corneocytes (brick) and lipid matrix (mortar).1,2 The corneocytes are mostly impermeable and molecules penetrate through lipid matrix. These lipid matrix are the key determinant for the barrier functions and is mainly composed of heterogeneous mixture of long chain ceramides (CER), cholesterol (CHOL) and free fatty acid (FFA) in certain ratios14,15 .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 is advised to read these articles

5,16

. Slight changes in the

head group of CER (either position or number of OH group) can change the morphology of bilayer significantly17. Atomistic simulation of molecular model of heterogeneous mixture of skin’s SC layer incorporating all possible CER, FFA and CHOL is challenging and beyond current computational capabilities. However, in order to simulate a realistic SC layer18, we have chosen the most abundant ceramide, CER-NS 24:0 and free fatty acid, FFA 24:0 for molecular model. Here 24 and 0, represent the length and number of double bonds in fatty acid chain, both in CER-NS and FFA, respectively. Since there is only single type of ceramide and fatty acid are

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

present in our simulations, here after CER-NS 24:0 and FFA 24:0 will be referred as CER and FFA, respectively.

Phospholipids have been studied from the last several decades. These efforts include development of force field parameters19 and permeation calculation of different solutes and drug molecules using umbrella sampling and constrained simulation20-26. Researchers have started looking at skin lipids from last few decades and significant progress have been made, both computationally and experimentally. NMR, DSC, WAXD, SAXD, 2H-NMR and neutron diffraction techniques have been used extensively to study the lipid organization27,28, phase behavior29,30 of skin SC lipids at different composition31, chain length18 and associated barrier properties32,33. Effect of penetration enhancers such as ethanol34 and DMSO35 has also been studied. Computationally, first simulation of FFA and CHOL mixtures at skin temperature was performed by Holtze et al.36. Pandit and Scott37 simulated the CER-NS (16:0) bilayer at high temperature of 368 K. Effect of permeation enhancer oleic acid on model skin membrane was performed by Hoopes at al.38 Notman et al.39 found out that at mole fraction of 0.4, DMSO transform the pure CER-NS (24:0) bilayer from gel to liquid phases and increases the permeability. First simulation incorporating all three component of SC, CER NS (24:0), FFA (24:0) and CHOL was performed by Das et al.40 to study the effect of molar ratio of individual component on structural properties. Same group also calculated diffusion and permeability of water through pure CER and mixed bilayer system41. Gupta et al.42 reported on the phase transition of skin lipids at various compositions and temperature. Recently, Martins et al.43 have demonstrated via both experiments and coarse grained MD simulation that high molecular weight proteins coated with hydrophobic surfactant molecules can effectively permeate the skin. In this study, we have performed several independent, long (~upto 40 ns) constrained united atom MD simulations to compute permeability of water, ethanol, DMSO, toluene, benzene, acetic acid, urea, phenol and styrene in equi-molar mixed bilayer of CER, CHOL and FFA.

Simulation details 2.1 Force field The force field parameters for the CER were taken from the Berger et al and GROMOS87 parameters44. The Berger force field has been parameterized for palmitoylsphingomyelin (PSM)

ACS Paragon Plus Environment

Page 4 of 26

Page 5 of 26

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

and phospholipid bilayes45. The PSM has similar structure as CER-NS except for head group geometry. The CH2 and CH3 groups are treated as single united atom C entity with zero partial charge and the polar hydrogen’s are included explicitly. The charges on polar groups were taken from earlier simulations16,

39-42

. The Ryckaert-Bellemans dihedral potential was used for the

hydrocarbon chains of FFA and CER46. The parameters for FFA and CHOL were taken from the Holtze et al36. The simple point charge (SPC) model was used for water molecule47. Gromos 54a7 and 54a6 parameter were used for each permeate. The force field parameters for each molecule are given in the supplementary information.

2.2 Simulation setup The simulations were carried out in NPT ensemble using the latest GROMACS molecular dynamics package48-51. The temperature was controlled at 310 K by Nose-Hoover thermostat with a time constant of 0.5 ps and thermostat coupled separately to lipid molecules and water. Pressure was controlled at 1 bar by Parrinello-Rahman barostat with a time constant of 5 ps and compressibility of 4.5 x 10-5 bar with semi isotropic coupling (XY and Z direction coupled separately) to achieve tensionless bilayer. All the bonds in lipid and solute molecules were constrained using LINCS algorithm52 while SETTLE algorithm was used for water53. A time step of 2 fs was used for all simulations. A cut off of 1.2 nm was used for vander Waal and electrostatic interactions. The systems were periodic in all three directions.

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

Figure 1. Molecular structure of individual skin lipid molecules used in simulations. From left: Ceramide-NS, cholesterol and free fatty acid. Images were created using visual molecular dynamics VMD software.74

Earlier studies on phospholipids showed that cutoff for electrostatic interaction leads to some artifact54. Thus in order to avoid these, long range electrostatic interactions were computed using particle mesh Ewald (PME) method. The equilibration time for all the systems was around (~ 1020 ns), according to the longest relaxation correlation time. The final 30 ns run of constrained simulation was used in calculation of constrained force.

2.3 Initial bilayer and solute structures The individual skin lipid constituents used in this study, are shown in Figure 1. The equimolar bilayer structure which was equilibrated for ~200 ns, was taken from the earlier simulation42. The bilayer consist of 154 lipids (52 CER, 50 CHOL and 52 FFA) and 5120 water molecules. The bilayer size was 4.9 nm x 4.9 nm x 11.7 nm. Figure 2 shows the density profile of individual bilayer constituents along the bilayer normal at 310 K, calculated using 200 ns unconstrained simulation42. The membrane interior is very much heterogeneous and sampling of such multi component interior is challenging. The CHOL mostly sits in between the middle of

ACS Paragon Plus Environment

Page 6 of 26

Page 7 of 26

each leaflet (z ~ -1.8, 1.8), while there is small increase in the density of FFA and CER in the middle of bilayer (z ~0). To sample the whole phase space of the bilayer, multiple solutes needs to be constrained at different xy and z plane. We have applied an approach in which multiple molecules are constrained in single window based on the cutoff used for electrostatic and vander Waal interactions (Figure 3). The reaction coordinate of the system was chosen to be membrane normal z, where z = 0 nm corresponds to the center of mass (COM) of the lipid molecules (CER+CHOL+FFA). For each system multiple solute molecule were inserted manually inside the upper water layer at the distance of ~ 4.8 nm from the COM of the bilayer. Overlapped water molecules with constrained molecules were removed and system was then energy minimized. Each system was then equilibrated for another 20 ns by keeping solute molecules fixed at the current position. The solute molecules was then pulled towards the center of the bilayer (-z direction) at the rate of 0.125 nm/ps using pull code of GROMACS. The time step for this run was kept at 1 fs. As the z distance between the center of mass (COM) of the lipid and solutes molecules changes by 0.2 nm, the configuration was stored and used for constrained simulation. In each window, four solutes were constrained at different z and xy position as shown in Figure 3. Total 25 windows were generated using the above procedure. These windows span the whole phase space from the bulk water from upper leaflet to the middle of the bilayer.

2000

3

Density (kg/m )

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

CER CHOL FFA LIPID water system

1500

1000

500

0

-4

-2

0

2

4

Z(nm) Figure 2: Density of individual skin lipid constituents along the bilayer normal of eqimolar mixed bilayer of CER, CHO and FFA, calculated from 200 ns unconstrained simulation. For the color code please refer to web version of article.

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 26

Figure 3. Side view of simulation snapshot of multiple molecules constrained in single window. The molecules are separated at least 1.2 nm in Z direction from each other. Images were created using visual molecular dynamics VMD software.74

The stored equidistance configuration was further ran for 30-35 ns. The distance between the COM of solute and COM of bilayer was constrained in z direction while solute was free to move in lateral direction. The configuration was stored at every 0.5 ps and constrained force was stored at every 10 fs. Last 25 ns runs of each simulation were used to calculate diffusion and potential of mean force. 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 solute 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.

2.4 Potential of mean force and diffusivity 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, then diffuses through the membrane interior, and finally dissolves again in the outer surrounding medium. The permeability is given by

p=

KD d

ACS Paragon Plus Environment

1

Page 9 of 26

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

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 composition changes the permeability many folds55 as well as local partition function K(z) governs the solute permeability56. Molecular dynamics simulation provides solution to calculate K(z) and D(z) along the bilayer normal ‘z’. Marrienk et al.20,21 calculated permeability of water molecules through DPPC bilayer and proposed four region model. We have adopted the similar strategy to compute the diffusivity, partition function, potential of mean force (PMF) and permeability.

The permeability for a symmetric bilayer system, which has normal in z direction, can be given by d /2

1 1 = ∫ dz P − d /2 K ( z ) D ( z )

2

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

D( z ) =

( RT )





0

2

< (∆F ( z , t )∆F ( z , 0)) >dt

3

 −∆G ( z )  K ( z ) = exp    kT 

4

∆F ( z , t ) = F ( z , t ) − < F ( z , t ) >

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. Now onwards, △G is referred as dG throughout the text, until unless mentioned otherwise.

2. Result and Discussion 3.1 Free energy of permeation

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Free energy of permeation profiles (dG) of the solutes studied, calculated using Eq 6, are plotted in Figure 4. 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 dG are calculated by propagating the errors in the force F(z). A general trend can be highlighted from the Figure 4. For all permeants the free energy is almost zero in bulk water, and as permeants reach near the head group of the bilayer, the free energy for hydrophilic compound slightly decreases, whereas, it slightly increases for the hydrophobic permeants. The reason could be the partially charged head group (hydrophilic in nature) of bilayer favorable for hydrophilic permeants as compared to hydrophobic one. After crossing the head group region the reverse trend was observed in free energy for both hydrophobic and hydrophilic permeants. In this region, free energy increases for hydrophilic compound and reaches a maximum and then slightly decreases, whereas for hydrophobic permeants it decreases and then slightly increases or remains same. For the former one, hydrophobic lipid tails of bilayer provides main barrier for permeation, while, for the latter, favorable condition are present in the interior of the bilayer.

dG (kJ/mol)

60

60

dG (kJ/mol)

60

oxygen

40

20

20

20

20

0

0

0

-20

-20

-20

0

water

-20

2

4

-4 -2 0

2

40

60

ethanol

40

-4 -2 0

4

-4 -2 0

2

40

4

60

40

40

20

20

20

20

0

0

0

-20

-20

-20

0

urea

40

butanol

40

2

4

-4 -2 0

2

4

60

toluene

benzene

4

40

40

-4 -2 0

2

4

60

phenol

40

-4 -2 0

styrene

40

20

20

20

0

0

0

0

-20

-20

-20

-20

2

4

-4 -2 0 2 Z (nm)

4

2

4

60

20

-4 -2 0

2

60

DMSO

-20

-4 -2 0

60

acetic acid

-4 -2 0

60

60

dG (kJ/mol)

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 26

-4 -2 0

ACS Paragon Plus Environment

2

4

ethylbenzene

-4 -2 0

2

4

Page 11 of 26

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. Free energy of permeation dG(z) along the bilayer normal for each permeants calculated from multiple molecule in single window constrained simulation.

The free energy profile for water has the same shape as obtained by the Das et al. at 300 K and 350 K for the water permeation from the model skin membrane41. The dG(z) increased smoothly from head group region and reached to maximum (~ 38.2 kJ/mol) in hydrophobic core and then decreased (~21 kJ/mol) at the center of the bilayer. The dG(z) increased due to the loss of favorable electrostatic interactions and hydrogen bonds when leaving from the bulk water phase to the hydrophobic interior of the bilayer. In contrast, in the middle of the bilayer, dG(z) again decreased due to the lower local lipid density as shown in Figure 2. The maximum values of dG(z) (~22.9-26 kJ/mol) for the water permeation through DPPC bilayer has been reported earlier20, 23. 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. The local density of the bilayer interior was found to be around ~730 kg/m3 which is comparable to density of hexadecane of 765.31 kg/m3 at 308 K57. The free energy of solvation of water in liquid hexadecane at 310 K is ~20 kJ/mol58 which agrees with the bilayer mid-plane values as shown in figure 4. Excess solvation Free energy (from water to hexadecane) of benzene, toluene, and butanol are 12.25, -15.33, and -8.914 kJ/mol respectively59. In our simulation we found ~10.1, -11.2 and -8.4 kJ/mol for benzene, toluene, and butanol, respectively. The simulated values are higher than that of experimental values. One possible explanation could be the effect of chain packing of the bilayers. Smaller permeants like water can easily be accommodated in the interior of the bilayer but same is not possible for bigger permeated like benzene, toluene and phenol etc. An additional work will be required to accommodate these bigger molecules, thus costing more free energy.

3.2 Diffusion coefficient Diffusion coefficient D(z) profiles for each permeate, calculated using Eq 3, are plotted in Figure 5. 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.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

The shape of D(z) profile for each permeants 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 permeate moves down towards the bilayer from the water phase, the Dz value decreases first and reaches a minimum value and then it increases till the middle of the bilayer. This is may be due to more free volume or less density at the middle of the bilayer as shown in Figure 2. Water, the smallest permeate, has highest D(z) values for any given z position. However we have not observed significant difference in the D(z) of bigger permeants like benzene, ethyl benzene and toluene etc. The D(z) values for water in bulk is ~4.8 x 10-5 cm2/s, which is comparable to the value (~5.0 x 10-5 cm2/s) obtained by Das et al. at 300K41. The experimental value of water self-diffusion coefficient is ~3.0x10-5 cm2/s60. SPC water model predicts diffusivity almost ~1.3 times more than that of experimental values 20. -4

-4

1×10

-4

1×10

1×10

acetic acid

ethanol

2

D(z) cm /s

-4

1×10

oxygen -5

-5

1×10

-6

-5

1×10

-6

-4 -2 0 2 4

-4

1×10

1×10

-6

-4 -2 0 2 4

-4

1×10

1×10

-6

-4 -2 0 2 4

-4

1×10

urea

1×10

-4 -2 0 2 4

-4

1×10

1×10

benzene

butanol

DMSO

2

D(z) cm /s

-5

1×10

water 1×10

-5

-5

1×10

-5

1×10

-6

1×10

-4

1×10

1×10

-6

-4 -2 0 2 4

-4

1×10

-6

-4 -2 0 2 4

-4

1×10

1×10

-4 -2 0 2 4

-4

1×10

1×10

styrene

phenol

ethylbenzene

2

toluene

-5

1×10

-6

-4 -2 0 2 4

1×10

D(z) 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 12 of 26

-5

-5

1×10

-6

1×10

-5

1×10

-6

-4 -2 0 2 4

1×10

-5

1×10

1×10

-6

-4 -2 0 2 4 Z (nm)

1×10

-6

-4 -2 0 2 4

ACS Paragon Plus Environment

1×10

-4 -2 0 2 4

Page 13 of 26

Figure 5. Diffusion coefficient D(z) along the bilayer normal for each permeate calculated from multiple molecule in single window constrained simulation. The Diffusion values in water for oxygen, benzene, acetic acid and urea are 2 x 10-5 cm2/s, 1.02 x 10-5 cm2/s, 1.21 x 10-5 cm2/s and 1.38 x 10-5 cm2/s respectively. We have obtained 3.8 x 10-5 cm2/s, 1.25 x 10-5 cm2/s, 2.4 x 10-5 cm2/s and 2.5 x 10-5 cm2/s for oxygen, benzene, acetic acid and urea respectively 61. Benzene diffusion value of 1.52 x 10-5 cm2/s in bulk water was reported earlier by Bempoard et al.23. Thus our diffusion values are in good agreement with previous experimental and computational studies. Figure 6 shows the diffusion coefficients of each permeants in bulk water, bilayer interior and the middle of the bilayer. The bulk water diffusion coefficient is much more size dependent as compared to one in the interior and as well as at the center. The diffusion coefficient decreases with increase in the size of permeate.

6.1e-05 3.1e-05 2

D (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

The Journal of Physical Chemistry

water phase densed lipid chain middle of the bilayer

1.5e-05 7.6e-06 3.8e-06 1.9e-06 10 20 30 40 50 60 70 80 90 100110120 MW

Figure 6. Diffusion coefficient of each permeate at three different position along the bilayer normal, calculated from multiple molecule in single window constrained simulation. The water

ACS Paragon Plus Environment

The Journal of Physical Chemistry

phase, dense lipid region and middle of the bilayer are correspond to z = ~ 4.6-4.8, z = ~ 1.2-1.4 and z = ~ 0.0 respectively. Please refer to web version of the article for the color code.

3.3 Resistance to Permeation Resistance to permeation R(z) profiles for each permeate, calculated using R(z) = 1/K(z)D(z), are plotted in Figure 7. The error bars are not shown for the clarity purposes. The error bars higher due to exponential dependence of K(z) or dG(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. For hydrophilic permeants, the main resistance of permeation was observed in the interior of the bilayer because of less favorable solute partition and diffusion coefficient as compared to water phase.

14

10 12 10 10 10 8 10

14

10 12 10 10 10 8 10

water

6

10 4 10 2 10 0 10

-4 -2

0

2

4

10 12 10 10 10 8 10

6

urea -4 -2

0

2

4

-4 -2

0

2

4

10 12 10 10 10 8 10

toluene

10 12 10 10 10 8 10

butanol

-4 -2

0

2

4

-4 -2

0

2

4

10 4 10 2 10 0 10

-4 -2

0

2

4

10 12 10 10 10 8 10

-4 -2 0 2 Z(nm)

4

-4 -2

0

2

4

2

4

10 12 10 10 10 8 10

benzene

6

-4 -2

0

2

4

10 4 10 2 10 0 10

dmso -4 -2

0

14

10 12 10 10 10 8 10

styrene

6

10 4 10 2 10 0 10

10 4 10 2 10 0 10 14

14

phenol

acitic acid

6

6

6

10 4 10 2 10 0 10

10 4 10 2 10 0 10 14

14

6

10 4 10 2 10 0 10

10 4 10 2 10 0 10

14

10 12 10 10 10 8 10

ethanol

6

6

14

10 12 10 10 10 8 10

10 4 10 2 10 0 10 14

10 12 10 10 10 8 10 10 4 10 2 10 0 10

14

10 12 10 10 10 8 10

oxygen

6

14

R(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 26

ethylbenzene

6

-4 -2

0

ACS Paragon Plus Environment

2

4

10 4 10 2 10 0 10

-4 -2

0

2

4

Page 15 of 26

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. Resistance of permeation R(z) along the bilayer normal for each permeants calculated from multiple molecule in single window constrained simulation.

The effect of diffusion coefficients looks almost negligible. For hydrophobic permeants the R(z) profiles are completely opposite to hydrophilic permeants. The main resistance of permeation was observed at lipid-water interface or near the head group region which is partially charged and polar. Both bilayer interior and center of bilayer significantly contribute to R(z), whereas in hydrophilic it is mostly dependent on G(z). Bempoard et al.23 obtained values of R(z) of water and benzene in DPPC bilayer in the order of ~105 and ~108 respectively, whereas we have obtained in the order of ~106 and ~109 for benzene and water, respectively, which shows the less permeability of skin membrane as compared to biological membrane.

3.4 Permeability The permeability values (P) calculated using integration of R(z) along the bilayer normal, are reported in Table 1, corresponding experimental data is also reported (wherever available). The water permeability was found to be (7.08±1.18) x 10-5 cm/s at 310 K which is almost 2 orders of magnitude higher than experimental values. Bemporad et al.23 also obtained an order higher water permeability for DPPC bilayer at 350 K, as compared to experiment. The experiments performed on porcine and human stratum-corneum were found to give P values in order of ~ 10-7 cm/s

62, 63

. 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 significantly17,18,64. Being highly hydrophilic, urea exhibits lowest permeability and its order is also matching with available experimental data. The obtained permeability of DMSO and ethanol are less as compared to water but reverse phenomena has been observed in experiments (Table 1). The reason could be permeation enhancing properties of DMSO and ethanol. It has been shown earlier that ethanol and DMSO fluidize the skin bilayer and increase the permeability34,35. The permeability of hydrophilic permeants mostly depends upon the dG(z) as compared to D(z).

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 16 of 26

For the hydrophobic molecules obtained P values are very much higher than experimental values. But, the rank ordering of permeability is same as in experiments. The P values were found to be highest for oxygen and minimum for styrene. For hydrophobic permeants, the main resistance lies at the head group and they are easily partitioned between water and lipid interior.

Table 1. Experimental and computed permeability of each permeants at 310 K.

Permeants

logPoct/w

MW

log Pexp

Pexp (cm/s)

cm/s Urea

DMSO

-2.11

-1.35

60

78

Ethanol -0.31

46

Acetic acid

Water

Oxygen

4.07 x 10

-8 a

-7.41b

3.89 x 10-8

-7.93c

1.18 x 10-8 c

-

-

-6.56a

2.72 x 10-7 a

-6.53b

2.95 x 10-7 b

-6.08c

8.32 x 10-7 c

-7.01a

0.98 x 10-7 a

59

-6.53b,c

2.95 x 10-7 b,c

-

18

-6.53b

2.95 x 10-7 b

-6.28c

5.24 x 10-7 c

-

-

-6.16a

6.92 x 10-7 a

-6.08b

8.32 x 10-7 b

-5.70c

0.20 x 10-7 c

-5.64a

2.29 x 10-6 a

-5.52b

3.20 x 10-6 b

-5.27c

5.37 x 10-6 c

-4.51a

3.09 x 10-5 a

-4.41b

3.89 x 10-5 b

-4.27c

5.37 x 10-5 c

-3.56a

2.75 x 10-4 a

-

32

0.839

74

Phenol 1.46

94

Benzene 2.13

78

PMD

cm/s

cm/s

-7.71

b

-0.17

Butanol

Toluene

-7.39

a

logPMD

ACS Paragon Plus Environment

(1.95 ± 0.98) x 10-8

-6.96

(1.01±0.23) x 10-7

-4.95

(1.12±0.18) x 10-5

-4.78 (1.65 ±0.27) x 10-5 -4.15 (7.08±1.18) x 10-5 1.40

25.11 ± 4.18

0.67

4.67 ± 0.78

0.41

2.57 ± 0.43

0.28

1.91 ± 0.32

Page 17 of 26

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

2.69

Styrene

2.95

92

104

Ethyl benzene 3.27

a

2.63 x 10-4 b

-3.64c

2.29 x 10-4 c

-3.75a,b

1.78 x 10-4 a,b

-3.48a

3.13 x 10-4 a

-3.36b

4.36 x 10-4 b

-3.00c

5.20 x 10-4 c

-0.25

0.56 ± 0.19

-0.85

(1.41 ± 0.3) x 10-1

-0.97

(1.07 ± 0.17) x 10-1

from the literature of X.C.Fu et al.65

b c

106

-3.58b

from the literature of Alan R. Katritzky et al.66

from the literature of Michael Abrham et al.67

The combined effect of thermodynamics (dG(z)) and kinetics (D(z)) lead to higher P values of oxygen. From these observation it can be deduced that the permeability of hydrophobic compound depends upon both dG(z) and D(z). Earlier coarse grained MD simulation of hydrophobic nano-particle across model lipid membrane also predicted very high permeability 6870

. Higher magnitudes of permeability could have been observed due to many factors. One of the

most important factor could be force field employed here for lipid molecules but this force field have shown realistic physical picture of the system as well as reproduced experimental quantities16,39-42. The parameters for permeants used from the similar force field which is parameterized for proteins and lipids, seems to be consistent. The value of was found to be sensitive to the position of the permeants in the bilayer and therefore improper sampling of phase space could lead to wrong free energy and permeability. But to ensure the full phase space sampling the 4 molecules constrained in each window of 5 independent simulations. Recently it has been shown that the inadequate sampling may lead to the inaccurate free energy (magnitude) which can further lead to an order of magnitude of discrepancy in permeability (because of exponential dependence).71-72 To check this, we have performed additional constrained simulation of water and benzene systems for 30 more ns (Please see supplementary information, Figure S2). But we have not observed much change in the free energy profile as well as in the values (with in the error bar). One should note that we have sampled only one reaction coordinate (translation along the bilayer normal) here, but in reality there are multiple reaction coordinates in the system which could be a reason for the discrepancy observed in the

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 26

permeability values. Advance techniques like metadynamics should be employed to explore all possible reaction coordinates in the system.73

Conclusion In this article, constrained molecular dynamics simulation of water, oxygen, ethanol, acetic acid, urea, butanol, benzene, DMSO, toluene, phenol, styrene and ethyl benzene, were performed to calculate the permeability coefficients of these permeants through mode skin lipid membrane and compared it with experimental data. The skin lipid membrane was made up of equimolar mixture of CER, CHOL and FFA. The inhomogeneous solubility-diffusion model used in this study incorporates both the equilibrium (free energy) and the dynamic (diffusion coefficient) behavior of permeants in the bilayer. The main resistance of permeation for hydrophilic permeants was obtained in the interior of the bilayer, whereas, for hydrophobic permeants, the same was observed near the lipid-water interface. The diffusion profiles for each permeate followed similar pattern. Smaller molecules exhibited higher D(z) values as compared to bigger one. Calculated permeability values are generally 1-2 orders of magnitude higher for hydrophilic permeants but qualitatively rank ordering of molecules based on octanol-water partition coefficient, is matching. The permeability for hydrophilic and hydrophobic molecules showed different dependency on both free energy and diffusion coefficient. Former is mostly depended upon the free energy while the latter is a function of both free energy and diffusivity but with a higher role of diffusivity. In experiments, membrane is treated as uniform slab barrier which is far from the reality. Using the given constrained MD simulations method, properties like free energy, diffusion, and local resistance in different regions of the lipid bilayer, can be studied at a molecular level. This technique explores only one reaction coordinate of the system and advanced techniques like metadynamics73 could be used for exploring all possible reaction coordinates in the system. This can further be leveraged for the drug or cosmetic design/testing/screening. Presented skin model only considered single type of ceramide (CERNS). We are currently working on development of coarse grained model of other CER’s and on development of an exhaustive skin model incorporating different kinds of CER, which will be communicated in separate article.

ACS Paragon Plus Environment

Page 19 of 26

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

Supporting Information. S1. Initial position of solutes in constrained simulation. S2. Effect of sampling time over free energy profile S3. Force field parameters of each solute.

Acknowledgement The current work was funded by TCS-CTO organization. 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 and Dr. Pradip, Vice president and Head Process Engineering Innovation Labs, TCS-TRDDC, Pune for their constant encouragement and support during this project.

References 1. Elias, P. M. Epidermal Lipids, Barrier Dunction, and Desquamation. J. Invest. Dermatol.

1983,80. 2. Pirot, F.; Kalia, Y.; Stinchcomb, A.; Keating, G.; Bunge, A.; Guy, R. Characterization Of The Permeability Barrier Of Human Skin In Vivo. Proceedings of the National Academy of Sciences 1997, 94, 1562-1567. 3. (2) Mitragotri, S.; Anissimov, Y.; Bunge, A.; Frasch, H.; Guy, R.; Hadgraft, J.; Kasting, G.; Lane, M.; Roberts, M. Mathematical Models Of Skin Permeability: An Overview. Int. J. Pharm. 2011, 418, 115-129. 4. Torin Huzil, J.; Sivaloganathan, S.; Kohandel, M.; Foldvari, M. Drug Delivery Through The Skin: Molecular Simulations Of Barrier Lipids To Design More Effective Noninvasive Dermal And Transdermal Delivery Systems For Small Molecules, Biologics, And Cosmetics. Wiley Interdiscip. Rev.: Nanomed. Nanobiotechnol. 2011, 3, 449-462. 5. Bartosova, L.Bajgar, J. Transdermal Drug Delivery In Vitro Using Diffusion Cells. Curr Med. Chem. 2012, 19, 4671-4677. 6. Notman, R.; Anwar, J. Breaching the Skin Barrier--Insights from Molecular Simulation of Model Membranes. Adv. Drug. Deliv. Rev. 2012, 65, 237-250.

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 26

7. 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. 8. Christian, U.; Hansen, C. M.; Dyk, V.; John, W.; Jensen, P. O. Permeability of Commercial Solvents through Living Human Skin. Am. Ind. Hyg. Assoc. J. 1995, 56, 651−660. 9. Ghosh, B; Reddy, L. H.; Kulkarni, R. V.; Khanam, J. Comparison of Skin Permeability of Drugs in Mice and Human Cadaver Skin. Ind. J. Exp. Biol. 2000, 38, 42-45. 10. Chantasart, D.; Li, S. Structure Enhancement Relationship Of Chemical Penetration Enhancers In Drug Transport Across The Stratum Corneum. Pharmaceutics 2012, 4, 71-92. 11. 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. 12. Mathur, V.; Satrawala, Y.; Rajput, M. S. Physical and Chemical Penetration Enhancers in Transdermal Drug Delivery System. Asian J. Pharm. 2010, 4, 173-183. 13. Swift, R. V.; Amaro, R. E.; Back to the Future: Can Physical Models of Passive Membrane Permeability Help Reduce Drug Candidate Attrition and Move Us Beyond QSPR? Chem. Biol. Drug Des. 2013, 81, 61-71. 14. Norlen, L.; Nicander, I.; Rozell, B. L.; Ollmar, S.; Forslind, B. Inter- and Intra-Individual Differences in Human Stratum Corneum Lipid Content Related to Physical Parameters of Skin Barrier Function In Vivo. J. Invest. Dermatol. 1999, 112, 72 77. 15. Weerheim, A; Ponec, M. Arch. Determination of Stratum Corneum Lipid Profile by Tape Stripping in Combination with High-Performance Thin-Layer Chromatography. Dermatol. Res.

2001, 293, 191-199. 16. Das, C.; Massimo, G. N.; Olmsted, P.D. Lamellar and Inverse Micellar Structures of Skin Lipids: Effect of Templating. Phys. Rev. Lett. 2013, 111, 148101 17. Engelbrecht, T.; Schroeter, A.; Hauß, T.; Demé, B.; Scheidt, H.; Huster, D.; Neubert, R. The Impact Of Ceramides NP And AP On The Nanostructure Of Stratum Corneum Lipid Bilayer. Part I: Neutron Diffraction And 2H NMR Studies On Multilamellar Models Based On Ceramides With Symmetric Alkyl Chain Length Distribution. Soft Matter 2012, 8, 6599.

ACS Paragon Plus Environment

Page 21 of 26

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

18. Mojumdar, E. H.; Kariman, Z.; van Kerckhove, L.; Gooris, G. S.; Bouwstra, J. A. The Role of Ceramide Chain Length Distribution on the Barrier Properties of the Skin Lipid Membranes. Biochim. Biophys. Acta. 2014, 1838, 2473–2483. 19. Egberts E.; Marrink, S. J.; Berendsen H. J. C. Molecular Dynamics Simulation of a Phospholipid Membrane. Eur Biophys J. 1994, 22, 423 436. 20. Marrink, S. J.; Berendsen H. J. C. Simulation of Water Transport through a Lipid Membrane. J. Phys. Chem. 1994, 98, 4155-4168. 21. Marrink, S. J.; Berendsen H. J. C. Permeation Process of Small Molecules across Lipid Membranes Studied by Molecular Dynamics Simulations. J. Phys. Chem. 1996, 100, 1672916738. 22. Tieleman, D.P.; Marrink, S. J.; Berendsen H. J. C. A Computer Perspective of Membranes: Molecular Dynamics Studies of Lipid Bilayer Systems. Biochim. Biophys. Acta.1997, 1331, 235270. 23. Bemporad, D.; Essex, J. W. Permeation of Small Molecules through a Lipid Bilayer: A Computer Simulation study. J. Phys. Chem. B. 2004, 108, 4875-4884. 24. Bemporad, D.; Luttmann, C.; Essex, J. W. Behaviour of Small Solutes and Large Drugs in a Lipid Bilayer from Computer Simulations. Biochim. Biophys. Acta. 2005, 1718, 1 – 21. 25. Nademi, Y.; Iranagh, S. A.; Yousefpour, A.; Mousavi, S. Z.; Modarress, H. Molecular Dynamics Simulations and Free Energy Profile of Paracetamol in DPPC and DMPC Lipid Bilayers. J. Chem. Sci. 2013, 126, 637–647. 26. Orsi, M.; Essex, J. W. Permeability of Drugs and Hormones through a Lipid Bilayer: A DualResolution Molecular Dynamics Study. Soft Matter 2010, 6, 3797–3808. 27. Groen, D., Gooris G. S., Bouwstra, J. A. New Insights into the Stratum Corneum Lipid Organization by X-ray Diffraction Analysis. Biophys J. 2009, 97, 2242-2249. 28. Groen, D.; Gooris, G. S.; Barlow, D. J.; Lawrence, M. J.; Van Mechelen J. B; Demé, B.; Bouwstra, J. A. Disposition of Ceramide in Model Lipid Membranes Determined by Neutron Diffraction. Biophys J. 2011, 100, 1481-1489. 29. de Jager, M.W.; Gooris, G.S.; Ponec, M.; Bouwstra, J.A. Lipid Mixtures Prepared with welldefined Synthetic Ceramides Closely Mimic the Unique Stratum Corneum Lipid Phase Behavior. J. Lipid. Res. 2005, 46, 2649-2656.

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 26

30. van Smeden, J.; Janssens, M.; Gooris, G.; Bouwstra, J. The Important Role of Stratum Corneum Lipids for the Cutaneous Barrier Function. Biochim. Biophys. Acta, Mol. Cell Biol. Lipids 2014, 295–313 31. Groen, D.; Poole D. S.; Gooris, G. S.; Bouwstra, J. A. Investigating the Barrier Function of Skin Lipid Models with Varying Compositions. Euro. J. Pharm. Biopharm. 2011, 79, 334-42. 32. de Jager, M.; Groenink, W.; van der Spek, J.; Janmaat, C.; Gooris, G.; Ponec, M.; Bouwstra, J. Preparation and Characterization of a Stratum Corneum Substitute for In Vitro Percutaneous Penetration Studies. Biochim. Biophys. Acta, Biomembranes 2006, 1758, 636-644. 33. de Jager, M.; Groenink, W.; Bielsa i Guivernau, R.; Andersson, E.; Angelova, N.; Ponec, M.; Bouwstra, J. A Novel in Vitro Percutaneous Penetration Model: Evaluation of Barrier Properties with P-Aminobenzoic Acid and Two of Its Derivatives. Pharm. Res. 2006, 23, 951-960. 34. Kwak, S.; Brief, E.; Langlais, D.; Kitson, N.; Lafleur, M.; Thewalt, J. Ethanol Perturbs Lipid Organization in Models of Stratum Corneum Membranes: An Investigation Combining Differential Scanning Calorimetry, Infrared and 2H NMR Spectroscopy. Biochim. Biophys. Acta, Biomembranes 2012, 1818, 1410-1419. 35. Marren, K. Dimethyl Sulfoxide: An Effective Penetration Enhancer for Topical Administration of NSAIDs. The Physician and Sportsmedicine 2011, 39, 75-82. 36. 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. Biochim. Biophys. Acta, Biomembranes 2001, 1511, 156-167. 37. Pandit, S.A.; Scott, H.L. Molecular Dynamics Simulation of a Ceramide bilayer. J. Chem. Phys. 2006, 124, 014708. 38. Matthew, I. H,; Massimo, G. N.; Marjorie, L. L.; Roland, F. Bilayer Structure and Lipid Dynamics in a Model Stratum Corneum with Oleic Acid. J. Phys. Chem. B. 2011, 115, 31643171. 39. 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. 40. Das, C.; Olmsted, P. D.; Noro, M. G. Simulation Studies of Stratum Corneum Lipid Mixtures. Biophys J. 2009, 97, 1941-1951.

ACS Paragon Plus Environment

Page 23 of 26

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

41. Das, C.; Olmsted, P. D.; Noro M. G. Water Permeation through Stratum Corneum Lipid Bilayers from Atomistic Simulations. Soft Matter 2009, 5, 4549-4555. 42. Gupta, R.; Rai, B. Molecular Dynamics Simulation Study of Skin Lipids: Effects of the Molar Ratio of Individual Components over a Wide Temperature Range. J. Phys. Chem. B. 2015, 119 (35), 11643–11655. 43. Martins, M.; Azoia, N.; Ribeiro, A.; Shimanovich, U.; Silva, C.; Cavaco-Paulo, A. In Vitro and Computational Studies of Transdermal Perfusion of Nanoformulations Containing a Large Molecular Weight Protein. Colloids Surf., B. 2013, 108, 271-278. 44. 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 45. 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 46. Ryckaert, J. P.; Bellemans, A. Molecular Dynamics of Liquid Butane near its Boiling Point. Chem. Phys. Lett.1975, 30, 123−125 47. Berendsen, H.; Postma, J.; Gunsteren, W.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces 1981, 331−342. 48.

Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message Passing

Parallel Molecular Dynamics Implementation. Comp. Phys. Comm. 1995, 91, 43–56. 49. van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible and Free. J. Comp. Chem. 2005, 26, 1701–1718. 50. Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. 51. J. C., Kasson, J.C.; P. M., van der Spoel, D., Hess, B., Lindahl, E. GROMACS 4.5: A Highthroughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 845–854. 52. Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472.

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 26

53. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Van Buuren, A.R.; Apol, E.; Meulenhoff, P.J.; Tieleman, D.P.; Sijbers, A.L.T.M.; Feenstra, K.A.; Van Drunen, R; Berendsen, H.J.C. Gromacs User Manual version 4.0, 2005. URL http://www.gromacs.org/Documentation/Manual. 54. 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. 55. 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. 56. 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. 57. Domańska, U.; Łachwa, J.; Morawski, P.; Malanowski, S. Phase Equilibria and Volumetric Properties in Binary Mixtures Containing Branched Chain Ethers (methyl 1,1-dimethylethyl ether or, ethyl 1,1-dimethyl ether, or ethyl 1,1-dimethylpropyl ether). J. Chem. Eng. Data. 1999, 44, 974-984. 58. Wu, Z.; Cui Q.; Yethiraj A. A New Coarse-Grained Force Field for Membrane–Peptide Simulations. J. Chem. Theory Comput. 2011, 7, 3793–3802. 59. Zhu, T.; Li, J.; Hawkins, G.; Cramer, C.; Truhlar, D. Density Functional Solvation Model Based on CM2 Atomic Charges. J. Chem. Phys. 1998, 109, 9117. 60. 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. 61. http://www.hypertextbookshop.com/biofilmbook/working_version/artifacts/tables/Module_0 4/Table4-1_DiffCoeffH2O.htm ref (accessed Mar 15, 2016). 62. Potts, R.; Francoeur, M. The Influence of Stratum Corneum Morphology on Water Permeability. J. Invest. Dermatol.1991, 96, 495-499. 63. Blank, I.H.; Moloney, J.; Emslie, A.G.; Simon, I.; Apt, C. The Diffusion of Water Across the Stratum Corneum as a Function of Its Water Content. J. Invest. Dermatol. 1984, 82, 188–194. 64. Š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.

ACS Paragon Plus Environment

Page 25 of 26

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

65. Fu, X. C.; Wang, G. P.; Wang, Y.F.; Liang, W.Q.; Yu Q.S.; Chow, M. S. S. Limitation of Potts and Guy’s Model and a Predictive Algorithm for Skin Permeability Including the Effects of Hydrogen-bond on Diffusivity. Pharmazie. 2004, 59, 282–285. 66. Katritzky, A.; Dobchev, D.; Fara, D.; Hür, E.; Tämm, K.; Kurunczi, L.; Karelson, M.; Varnek, A.; Solov'ev, V. Skin Permeation Rate as a Function of Chemical Structure. J. Med. Chem. 2006, 49, 3305-3314. 67. Michael, H.A.; Filomena, M. Human Skin Permeation and Partition: General Linear FreeEnergy Relationship Analyses. J. Pharma. Sci. 2004, 93, 1508–1523. 68. Lin, X. B.; Li, Y.; Gu, N. Nanoparticle’s Size Effect on its Translocation across a Lipid bilayer: A Molecular Dynamics Simulation. J. Comput. Theor. Nanosci. 2010, 7, 269–276. 69. Steven, L. F.; Angela, V. Simulation of Nanoparticle Permeation through a Lipid Membrane. Biophysical J. 2010, 99, 144-152. 70. Gupta, R.; Rai, B. Penetration of Gold Nanoparticles through Human Skin: Unraveling Its Mechanisms at the Molecular Scale. J.Phys. Chem. B 2016. DOI: 10.1021/acs.jpcb.6b03212 71. Neale, C.; Pomès, R. Sampling Errors in Free Energy Simulations of Small Molecules in Lipid Bilayers. Biochim. Biophys. Acta, Biomembranes 2016.doi:10.1016/j.bbamem.2016.03.006 72. Lv, C.; Aitchison, E.; Wu, D.; Zheng, L.; Cheng, X.; Yang, W. Comparative Exploration Of Hydrogen Sulfide And Water Transmembrane Free Energy Surfaces Via Orthogonal Space Tempering Free Energy Sampling. J. Comput. Chem. 2015, 37, 567-574. 73. Jämbeck, J.; Lyubartsev, A. Exploring The Free Energy Landscape Of Solutes Embedded In Lipid Bilayers. J. Phys. Chem. Lett. 2013, 4, 1781-1787. 74. Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Molec. Graphics. 1996, 14.1, 33-38.

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

TOC

ACS Paragon Plus Environment

Page 26 of 26