In-Silico Skin Model: A Multiscale Simulation Study ... - ACS Publications

Subscriber access provided by UNIV OF NEWCASTLE

Article

In-Silico Skin Model: A Multiscale Simulation Study of Drug Transport Kishore Gajula, Rakesh Gupta, D B Sridhar, and Beena Rai J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00224 • Publication Date (Web): 18 Jul 2017 Downloaded from http://pubs.acs.org on July 19, 2017

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.

Journal of Chemical Information and Modeling 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 19

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

Journal of Chemical Information and Modeling

In-Silico

Skin

Model:

A

Multiscale

Simulation Study of Drug Transport Kishore Gajula, Rakesh Gupta, D B Sridhar and Beena Rai* Physical Sciences 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 Accurate in-silico models are required to predict the release of drug molecules through skin in order to supplement the in-vivo experiments for faster development/testing of drugs. The upper most layer of the skin, stratum corneum (SC), offers main resistance for permeation of actives. Most of the SC’s molecular level models comprise of cholesterol, and phospholipids only, which is far from the reality. In this study we have implemented a multiscale modelling framework to obtain the release profile of three drugs namely Caffeine, Fentanyl and Naphthol through skin SC. We report for the first time diffusion of drugs through a realistic skin molecular model comprised of ceramides, cholesterol and free fatty acid. The diffusion coefficients of drugs in the SC lipid matrix were determined from multiple constrained molecular dynamics simulations. The calculated diffusion coefficients were then used in the macroscopic models to predict the release profiles of drugs through the SC. The obtained release profiles were in good agreement with available experimental data. The partition coefficient exhibits more effect on the release profiles. The reported multiscale modelling framework would provide insight into the delivery mechanisms of the drugs through the skin, and shall act as a guiding tool in performing targeted experiments to come up with a suitable delivery system.

Introduction Human skin protects body from the attack of bacteria and pathogens, provides permeation barrier to many harmful molecules and also maintains the hydration level. The outer layer of

ACS Paragon Plus Environment


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

skin epidermis, known as Stratum Corneum (SC), is mainly responsible for these barrier properties.1 Because of skin’s high surface area, it provides a convenient route for administration of drugs.1,2 Transdermal delivery provides a potential alternative to other regular modes of administration of drugs. Understanding the mechanism of drug permeation through the skin, and, dermal uptake, is beneficial for the development of effective transdermal delivery formulations. Currently, both pharma and cosmetic industries rely on the extensive in-vivo/in-vitro experiments to arrive at the optimized formulations, which are highly time consuming and expensive. In view of this, a realistic in-silico molecular model of human skin coupled with accurate transport models, which could supplement/replace some of the elaborate in-vivo/in-vitro tests with in-silico tests, are required. The in-silico skin model shall couple the atomistic/molecular level mechanisms with the continuum model of the skin SC so as to capture the trans-dermal delivery mechanisms. The molecular level decoding of the SC’s structure has enabled us to develop an appropriate virtual model to accurately mimic its barrier properties.3 The SC is made up of corneocytes and lipid matrix which are arranged in brick and mortar fashion, respectively.1,2 The corneocytes are mostly impermeable and molecules penetrate only through lipid matrix. The drug permeation across the SC happens largely through passive diffusion and is highly selective in nature. It only allows relatively lipophilic compounds to diffuse into its lipid layers. The main physical property that governs species transport in the skin is diffusion coefficient. Thus the lipids of the skin SC should be approximated with realistic composition for accurate calculation of the diffusion coefficients. The lipid matrix is composed of a heterogeneous mixture of long chain ceramides (CER), cholesterol (CHOL) and free fatty acid (FFA) in certain ratios.4,5 Based on the structure of head group and attached fatty acid chain length, 11 different types of ceramides have been found to be present in human SC lipid matrix..6,7 Atomistic molecular dynamics (MD) simulation of heterogeneous mixture of all possible CER, FFA and CHOL is challenging and beyond the current computational capabilities. However, in order to simulate a realistic SC layer8, most abundant ceramide, CER-NS and free fatty acid, FFA were chosen in this study. The diffusion and permeability of small molecules through skin lipid matrix were reported recently.9-11 The diffusion coefficient calculated from MD simulation, can further be used to predict the dermal uptake/cumulative release through the SC using computational fluid dynamics (CFD) techniques.


Page 2 of 19

Page 3 of 19

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 macroscopic model of skin SC consists of both corneocytes and lipid matrix in brick and mortar fashion. There are a few studies reported on the macroscopic modeling of skin.12-19 Lindstrom et al.12 used finite difference method to study the drug release from the suspension to skin. Burnette et al.13 studied the diffusion of the drugs through the stratum corneum (SC) using Monte Carlo Random walk method and compared their results with the solutions that of Fick’s 2nd law solved in the homogeneous membrane. Rim et al.14 studied the drug diffusion from finite drug reservoir patch into the skin using finite element method. Rim et al.15 carried out a multiscale simulation for studying transdermal drug delivery. Authors used MD simulations to determine the diffusion coefficient and finite element method to predict the drug flux and compared with the experiments. However in this study, for MD simulations, authors assumed the lipid of the skin to be composed of 1,2-Dimyristoyl-sn-glycero-3phosphorylcholine (DMPC)-cholesterol, which is far from the realistic composition of the skin lipids.15 Tojo et al.16 studied the release of the drug from a device into the skin which consisted of SC and viable skin layers subsequently into body. Kushner et al.17 modelled the skin SC lipid pathway from first principles approach to come up with the modified Fick’s 2nd law. Authors accounted for the branched pathways which exist in the SC and developed the two-tortuosity model.17 Permeation data from the experiments was fitted to the model to predict the diffusion coefficient and partition coefficient. The fitted diffusion coefficient is an approximate one because it doesn’t account for heterogeneity of the SC lipids. Frash et al.18 studied the diffusion of drug in SC lipid pathway using finite element method simulations and compared the steady state flux with that of homogeneous membrane. Authors have also used the finite element data to find the effective diffusion coefficient by fitting to a homogeneous model.18 Hasen et al.19 considered the biphasic SC and deep skin layers while performing the permeation experiments and fitted with the solution of Fick’s 2nd law in order to predict diffusion and partition coefficients. To the best of our knowledge, the only effort towards multiscale modelling of skin has been reported by Rim et al.15 Authors first calculated the microscopic diffusion coefficient of a drug in a lipid bilayer of SC using MD simulations. Authors assumed brick and mortar model while calculating the effective macroscopic diffusion coefficient. Authors calculated effective macroscopic diffusion coefficient using homogenization procedure20 where they account for microstructure of SC. The effective diffusion coefficient was used in FEM simulations to predict the release profile of active. Authors assumed that the bilayer consists of DMPC and cholesterol with 50 mol % cholesterol to mimic the skin lipid matrix, which is a crude



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

Page 4 of 19

approximation. The lipid matrix of skin is actually composed of long chain CER, CHOL and FFA not DMPC. The lateral packing of lipids plays crucial role in determining the bilayer barrier properties. Recently it has been reported that small chain ceramide bilayer were more permeable as compared to long chain ceramide bilayer, both by experiment21 and by MD simulation.11

The

MD

simulation

showed

that

the

packing

of

CHOL

with

Dipalmitoylphosphatidylcholine (DPPC)22 and CER23 were not similar and the bilayer properties were highly different for both the cases. In this study, we report a multiscale modelling framework which would provide insight into the delivery mechanisms of the drugs through skin, and will act as guiding tool in performing targeted experiments (derived from simulations) to come up with a suitable delivery system. We have modelled a realistic skin SC layer at both molecular and macroscopic scale. The diffusivity of drug molecules namely Fentanyl, Caffeine and Naphthol through the SC lipid layer was calculated using multiple atomistic MD simulations. The averaged diffusivity along the bilayer normal, was used as an input for the macroscopic CFD model. The macroscopic model deals with solving Fick’s 2nd law in the specified geometry (brick and mortar) with appropriate boundary conditions. To account for parallel and branched pathways in SC (lipid region), the modified Fick's 2nd law is used, which was developed from the first principles approach.17

2. Method, Model and Initial Structure. 2.1 Molecular Scale Simulation The force field parameters of the CER were based on Berger et al. and GROMOS87 parameters.24 The CH2 and CH3 groups were represented by united carbon atom with zero partial charge and the polar hydrogen atoms were included explicitly. The charges on polar groups were taken from earlier simulation study.9-11 The parameters for FFA and CHOL were taken from the Holtze et al.25 The simple point charge (SPC) model was used for water molecule.26 Gromos 54a7 and 54a6 parameter were used for drug molecules. The force field parameters for each drug molecule are given in the supporting information.


Page 5 of 19

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 lipids and drug molecules used in simulations. Oxygen, hydrogen, nitrogen and carbon atoms are shown in red, white, blue and cyan color, respectively. CH2 and CH3 are represented by single C (cyan) atom type. The simulations were performed in NPT ensemble using the latest GROMACS MD package.27-30 The temperature was controlled at 310 K by Nose-Hoover thermostat with a time constant of 0.5 ps and it was 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 drug molecules were constrained using LINCS algorithm31 while SETTLE algorithm was used for water.32 A time step of 2 fs was used for all simulations. A cut off of 1.2 nm was used for van der Waal and electrostatic interactions and later one was computed using particle mesh Ewald (PME) method. The equilibration time for all the systems was around (~ 10-20 ns), according to the longest relaxation correlation time. The final 25 ns run of constrained simulation was used in the calculation of constrained force. The systems were periodic in all three directions. The individual skin lipid constituents and drug molecules used in this study, are shown in Figure 1. The equimolar bilayer structure, which was equilibrated for ~200 ns, was taken from the earlier simulation.10,23 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. The sampling of heterogeneous membrane interior is very much challenging. The constrained MD simulations34,35 were performed in order to calculate transport properties across the skin



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

lipid bilayer. Recently, we have successfully modelled coarse grained skin lipid bilayer and calculated transport properties across this layer.36-38 To sample the whole phase space of the bilayer, multiple drug molecules were constrained at different xy and z plane. We have followed an approach in which multiple molecules are constrained in single window based on electrostatic and Vander Waal cutoff.10,11,33 The membrane normal z, was chosen as reaction coordinate of the system, where z = 0 nm corresponds to the center of mass (COM) of the lipid molecules (CER+CHOL+FFA). For each system, multiple drug molecules were inserted in the upper water layer at the distance of ~ 4.8 nm from the COM of the bilayer. Overlapping water molecules with constrained drug molecules were removed and system was then energy minimized. Each system was further equilibrated for another 20 ns by keeping drug molecules fixed at the current position. The drug molecules were then pulled towards the center of the bilayer (-z direction) at the rate of 0.1 nm/ps using pull code of GROMACS.

Figure 2. Final configurations of constrained drug molecules at different xy and z position in three different window. CER, CHOL, FFA and water are shown in red, green, orange and cyan color. Constrained molecule is shown in the vdW atomic style.

The time step for these runs 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 further constrained simulation. In each window, drug molecule was constrained at different z and xy position as shown in Figure 2. Total 25 windows were generated using the above procedure for each drug molecule which spanned the whole phase space from the


Page 6 of 19

Page 7 of 19

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


bulk water from upper leaflet to the middle of the bilayer. The stored equidistance configurations were further simulated for 30-35 ns by constraining the distance between the COM of drug molecules and COM of bilayer. But the drug 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. Last 25 ns runs of each simulation were used to compute the thermodynamics and transport properties of drug molecules through lipid layer. To improve the statistical accuracy of sampling, each drug molecule was constrained at four different xy plane for a given z coordinate. The results presented here, are averages of these four simulations. 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 diffusivity D(z) of drug molecules along the bilayer normal is given by,34,35

=

∆ ,∆ ,

∆ , = , −< , >

1 2

Where R is gas constant, T is temperature, F(z,t) is constrained force on solute at a given z and is averaged force over time.

2.2 FEM Simulation The diffusion of actives in intercellular region of SC is interest of this study, because the corneocytes are impermeable. The diffusion is governed by Fick’s 2nd law.

=

3

Where, C = concentration of active, Db = Diffusion coefficient of active in lipid region of SC and t = time. To account for the tortuous path of lipid region in the SC, Kushner et al.17 came up with modified Fick’s 2nd law which accounts for the tortuous path of diffusion.

= !$

"#

%&'( $)*&'+,

- ∇

4

Where, τflux, τvolume = tortuosity factors to account for parallel and branched transport or active in lipid region. Modified Fick’s 2nd law can be solved in FEM framework to get the concentration in the intercellular region of SC. The tortuous factors used in the modified Fick’s 2nd law can be calculated from the geometric parameters of the SC. τflux and τvolume are calculated from the geometric parameters of SC



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

/012345 = / ?@>

6

67869:;

Where, N = number of corneocyte layers, h = thickness of corneocyte, g = width of lipid channel, ω = offset ratio, d = corneocyte width. The modified Fick’s 2nd law is solved in the intercellular region of SC using FEM framework with proper boundary conditions in order to calculate the concentration at each point in SC intercellular area. The amount of active that permeated through the SC can be obtained from integrating the diffusive flux at the end of the SC with respect to time. Using modified Fick’s 2nd law the cumulative release of an active through lipid region of SC is given by17

A =

B

$%&'(

C

D, D

E

DF

G

7

Where, Q(t) = amount of active diffused through the lipid region of SC, τflux = tortuosity factor, Db = diffusion coefficient of active in lipid region of SC and H = porosity.

Figure 3. Geometry of skin SC lipid region considered in the simulation. Applied boundary conditions are also shown.

The geometry of SC intercellular region was constructed in COMSOL multiphysics software as shown in Figure 3. Geometric parameters of SC were obtained from the literature.39 The


Page 9 of 19

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


geometry was meshed using triangular elements with inbuilt meshing tool available in COMSOL software. The module of transport of diluted species in COMSOL was used for solving the modified Fick’s 2nd law in the geometry. The following boundary conditions were used in the simulation (please see Figure 3).

•

Concentration at inlet (donor side) C (t, y=L) = Co mol/m3 = Kvehicle/SC X Cvehicle

•

Concentration at outlet (receiver side) C (t, y=0) = 0 mol/m3

•

Initial concentration in the lipid region of SC C (t=0) = 0 mol/m3

•

Periodic boundary condition at both sides of lipid region of SC

Corneocytes were assumed to be impermeable. The Fentanyl, Caffeine and Naphthol drug molecules were used in the study. The release profile of these actives through epidermis of the skin were taken from the literature.15,17,40 The diffusion coefficient of the actives through skin lipids of the SC were obtained from MD simulations, as described in earlier sections. The partition coefficient of the actives from the formulation into SC was taken from the literature.15,17,40 This partition coefficient decides the amount of the active that enters into the SC from the formulation. The initial concentration of the actives were taken from the literature.15,17,40 The FEM simulation was run for the experimental time as reported in the literature. 15,17,40 The geometric parameters used in the simulation for a particular active were taken from the literature where the experimental release profile was obtained for the direct comparison with the simulation results.15,17,40 An absolute tolerance of 1x10-6 was used in simulations for all actives. The concentration gradient at the end of the SC was obtained from the FEM simulations. The amount of the active permeated through SC was calculated from the computed concentration gradient data.

3. Results and Discussion 3.1 Molecular Dynamic Simulation. Diffusion coefficient D(z) profiles for each drug molecule along the bilayer normal, calculated using equation 1, are plotted in Figure 4. To clearly show the difference in the diffusion values of drug molecules, the error bars have not been incorporated in Figure 4. The diffusion profiles for each drug molecule with error bars have been plotted in Figure S1 (please see supporting information). The shape of D(z) profile of each molecule was similar in nature. The D(z) values in bulk water was found to be greater than that of bilayer interior and it may be because of higher density of the lipid environment inside the bilayer. As the drug molecule moved down towards the bilayer from the water phase, the D(z) value



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

decreased and reached a minimum value and then it remained almost constant till the middle of the bilayer. It is interesting to note that the sudden drop in the D(z) for each drug molecule happened at the same z position. The D(z) values were decreased by an order of magnitude in the bilayer interior as compared to bulk water for each drug molecule.

Figure 4. Local Diffusion coefficient of Naphthol, Fentanyl and Caffeine molecules along the bilayer normal, calculated from the constrained molecular dynamics simulations. The bilayer normal is in z direction and z= 0 correspond to the center of the bilayer. The diffusion coefficient changed along the bilayer normal but was found to be constant in the bilayer interior, so for simplification, the average value along the bilayer normal were calculated and used in FEM simulations. It has been shown earlier that diffusivity of bigger drug and protein molecules33, 36-38 do not vary much in the interior of the skin lipid bilayer but for small molecule it varied significantly.9-11 Since, the drug molecules used here are bigger in size, so the obtained observations are well in agreement with previous computational studies.33,

36-38

The diffusivity average was taken from the z = -2.8 to 2.8. The averaged


Page 10 of 19

Page 11 of 19

diffusivity of Naphthol, Fentanyl and Caffeine in bilayer interior was found to be 2.12x10-10, 3.67x10-10 and 2.37x10-10 m2/s, respectively. These diffusion coefficients were further used in the FEM simulation to obtain the release profile of drug molecules through SC layer. 3.2 FEM Simulation: 3.2.1 Mesh Independent Study: Five different meshes G1, G2, G3, G4 and G5 were considered in the study. The number of triangular elements in the mesh was varied as 14354, 38448, 65368, 96188 and 117910, respectively. FEM simulations were performed using these five meshes for time of one hour. The cumulative release profiles for different meshes have been plotted in Figure 5.

Q(t) gm/cm2

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.8x10

-6

1.6x10

-6

1.4x10

-6

1.2x10

-6

-6

1x10

-6

8x10

-7

6x10

-7

4x10

-7

2x10

-7

9.580x10 -6 9.560x10 -6 9.540x10 -6 9.520x10 -7 9.500x10 -7 9.480x10 -7 9.460x10 -7 9.440x10

0.8 0.802 0.804 0.806 0.808 0.81

G1 G2 G3 G4 G5

0 0

0.2

0.4

0.6

0.8

1

t (hrs) Figure 5. Cumulative release profile through skin SC lipid region from FEM simulations for different meshes G1-G5. The number of triangular elements in the meshes were as 14354(G1), 38448(G2), 65368(G3), 96188(G4) and 117910(G5), respectively. These profiles were used to check the mesh independency of the FEM model. The difference between the cumulative release amounts from different meshes decreased with increase in the number of elements in the mesh. The difference between meshes G4 and G5 was almost



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

negligible. Thus in order to reduce the computational effort without losing much accuracy, the mesh G4 was considered for further simulations. 3.2.2 Model Results: The fundamental equations 4 to 7, from two-tortuosity model17 were used in macroscopic FEM model. The diffusion coefficients obtained from the MD simulation and partition coefficients obtained from the literature were used for solving the modified Fick’s 2nd law of two-tortuosity model. The concentration gradient data across the SC was exported. The cumulative amount of drugs permeated through SC was calculated by numerical integration of concentration gradient data across the SC. The release profile (cumulative amount of drug permeated through SC) of each drug molecule through skin SC is plotted in Figure 6. Release profiles were compared with experimental release profile obtained from the literature.15,17,40

Figure 6. (a) Cumulative release profile of Fentanyl through skin SC lipid region. FEM simulation vs. experiments.15 Figure 6. (b) Cumulative release profile of Caffeine through skin SC lipid region. FEM simulation vs. experiments.40 Figure 6. (c) Cumulative release profile of Naphthol through skin SC lipid region. FEM simulation vs. experiments.17


Page 12 of 19

Page 13 of 19

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 geometric parameters for all the above simulations were obtained from the literature.39 Kushner.et al17 reported that the value of ω (offset ratio) is 3-8 for native human skin samples. In case of the active fentanyl the symmetric arrangement of the corneocytes was considered as mentioned in the article15 from where the experimental release profile was compared and all other geometric parameters were taken from the literature.39 Dias.et.al40 performed experiments to study permeation of caffeine through human epidermis. No information about geometric parameters was available in the article. In case of fentanyl, the experimentally measured partition coefficient (0.14)15 was used in the simulation. The experimental fentanyl flux15 reported was integrated with respect to time in order to get the cumulative experimental release profile. The obtained simulation release profile of fentanyl was well in agreement with the experiments15, order of magnitude of release profile was comparable to that of experiments. Rim et al15 reported simulated flux for the time of 40 hrs, whereas the same authors reported the experimental flux for the time of 70 hrs. The reported simulated flux was integrated with respect to time to get the cumulative release profile and the same has been compared with our simulation release profile along with cumulative experimental release profile. Although both results cannot be compared completely as simulated time in both cases are different, the order of magnitude of our simulation result is comparable to the reported simulated release profile as well. The plots of reported experimental, simulation results and the current FEM simulations have been presented in Figure S4 for the comparison (please see the supporting information). In case of caffeine there are multiple partition coefficients reported in the literature.19,41 From the Potts and Guy correlation,41 partition coefficient was calculated to be 0.8918. Hansen et al.19 reported the value of partition coefficient as 2.15 and 4.51 for lipid/donor and SC/donor, respectively. These partition coefficients are fitted values of the experimental data to analytical expressions of concentration as depth of SC. A sensitivity analysis was performed using the different ω values and partition coefficients. The release profile obtained from FEM simulations for value of partition coefficient 2.15 and ω =3 compared well to the experimental release profile. In case of Napthol also there are multiple values of partition coefficients reported in the literature17 which are fitted values of simulation data to a known model. A sensitivity analysis was performed using the different ω values and partition coefficients. The release profile obtained from FEM simulations for value of partition coefficient 267 and ω =8 was closer to the experiment release profile. The deviation of



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

simulation results as compared to experimental results arise due to the value of the partition coefficient used in the simulation. The partition coefficient values reported in the literature varied with different literature sources.15,17,19 The partition coefficient between vehicle and SC bilayer can also be determined from well-known Potts and Guy correlation.41 The partition coefficient signifies the amount of active partitioning into SC layer from the formulation. Since the partition coefficient values used in our simulations are obtained from literature which in turn are fitted values from a homogeneous model, a shift of the release profile is observed which could arise due to error in the estimation of partition coefficient. The shift in the profiles of experimental vs. FEM simulation plotted in all simulations may also be due to the assumption in our simulations where the permeation through lipid region of the SC is only considered. In the experiments whole epidermis is used for obtaining permeation data. The details of sensitivity analysis are given in the supporting information.

4. Conclusions In this study, we have performed multiscale simulation (MD and FEM) of skin SC layer. The obtained release profile of different drug molecules through SC model shows good qualitative agreement with experiments. A realistic model of skin’s SC layer having equimolar mixture of ceramide, cholesterol and free fatty acid, was used. Several constrained MD simulations were performed to calculate the diffusion of three drug molecules namely Caffeine, Fentanyl and Naphthol. The diffusivity of the drugs changed with the bilayer depth. In the water layer, the diffusivity was ~5 times higher than that of bilayer interior. The diffusivity remains constant inside the bilayer interior. The averaged diffusivity from the water-bilayer interface to bilayer center was calculated and used as an input for the FEM simulation. The FEM simulation was performed on the whole SC layer comprising of lipid matrix and corneocytes. The corneocytes were assumed to be impermeable. The modified Fick’s second law equation17 was solved to obtain the cumulative amount of drug diffused through SC layer. It is important to note that integrated multi-scale model (molecular and macroscopic models) presented here does show a good qualitative match between experiments and simulations. However, there are certain limitations which need to be overcome to make it more robust and quantitative. It has been shown recently that the CER and FFA distribution changes the skin permeability and morphology significantly.20,42 Thus the molecular model must account for variability in skin lipid compositions and contain major CER and FFA present in SC. On the


Page 14 of 19

Page 15 of 19

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


other hand the macroscopic model suffers from the lack of reliable partition coefficients of actives between the vehicle/formulation and the SC. The model is dependent on the experimentally derived partition coefficients which are very few. It is heartening to note that in spite of these limitations, the presented multiscale modelling framework would provide useful insight into the delivery mechanisms of the drugs through the skin, and might act as a guiding tool in performing targeted experiments to come up with a suitable delivery system. The efforts are on at our research labs to refine the models further to address some of these shortcomings.

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 for his constant encouragement and support during this project.

Author Contributions R.G. and B.R. conceived the study. K.G. performed CFD simulation. R.G and D.S. performed MD simulations. All authors contributed to the interpretation and discussion of the results.

Supporting Information S1. Force field parameters for drug molecules; S2. Diffusion Profile with Error Bars and S3. Sensitivity Analysis

References 1. Elias, P. M. Epidermal Lipids, Barrier Dunction, Desquamation. J. Invest. Dermatol. 1983, 80, 44s-49s. 2. Pirot, F.; Kalia, Y.; Stinchcomb, A.; Keating, G.; Bunge, A.; Guy, R. Characterization Of The Permeability Barrier Of Human Skin In Vivo. Proc. Natl. Acad. Sci. 1997, 94, 15621567.

3. Iwai, I.; Han, H.; Hollander, L.; Svensson, S.; Öfverstedt, L. G.; Anwar, J.; Brewer, J.; Bloksgaard, M.; Laloeuf, A.; Nosek, D.; Masich, S. The Human Skin Barrier Is



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 19

Organized As Stacked Bilayers Of Fully Extended Ceramides With Cholesterol Molecules

Associated

With

The

Ceramide

Sphingoid

Moiety. J.

Invest.

Dermatol. 2012, 132, 2215-2225.

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

5. Bartosova, L.Bajgar, J. Transdermal Drug Delivery In Vitro Using Diffusion Cells. Curr Med. Chem. 2012, 19, 4671-4677

6. Weerheim, A.; Ponec, M. Determination of Stratum Corneum Lipid Profile by Tape Stripping in Combination with High-Performance Thin-Layer Chromatography. Arch. Dermatol. Res. 2001, 293, 191-199.

7. Das, C.; Noro, M.; Olmsted, P. Lamellar and Inverse Micellar Structures of Skin Lipids: Effect of Templating. Phys. Rev. Lett. 2013, 111, No. 148101.

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

9. Das, C.; Olmsted, P. D.; Noro M. G. Water Permeation through Stratum Corneum Lipid Bilayers from Atomistic Simulations. Soft Matter. 2009, 5, 4549-4555. 10. Gupta, R.; Sridhar, D. B; Rai, B. Molecular Dynamics Simulation Study of Permeation of Molecules through Skin Lipid Bilayer. J. Phys. Chem. B. 2016, 120, 8987-8996. 11. Gupta, R.; Sridhar, D. B.; Rai, B. Molecular Dynamics Simulation of Skin Lipids: Effect of Ceramide Chain Lengths on Bilayer Properties. J. Phys. Chem. B. 2016, 120, 1253612546. 12. Tom Lindstrom, F.; Ayres, J. W. Diffusion Model for Drug Release from Suspensions II: Release to a Perfect Sink. J. Pharm. Sci. 1977, 66, 662-668. 13. Burnette, R. R. A Monte-Carlo Model for The Passive Diffusion of Drugs Through the Stratum Corneum. Int. J. Pharm. 1984, 22, 89-97. 14. Rim, J. E.; Pinsky, P. M.; van Osdol, W. W. Finite Element Modeling of Coupled Diffusion with Partitioning in Transdermal Drug Delivery. Ann. Biomed. Eng. 2005, 33, 1422-1438. 15. Rim, J. E.; Pinsky, P. M.; van Osdol, W. W. Multiscale Modeling Framework of Transdermal Drug Delivery. Ann. Biomed. Eng. 2009, 37, 1217-1229. 16. Kakuji, T. Mathematical Modeling of Transdermal Drug Delivery. J. Chem. Eng. Jpn.

1987, 20, 300-308.


Page 17 of 19

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


17. Kushner, J.; Deen, W.; Blankschtein, D.; Langer, R. First‐principles, Structure‐Based Transdermal Transport Model to Evaluate Lipid Partition and Diffusion Coefficients of Hydrophobic Permeants Solely from Stratum Corneum Permeation Experiments. J. Pharm. Sci. 2007, 96, 3236-3251. 18. Frederick, F. H.; Barbero A. M.

Steady‐State Flux and Lag Time in The Stratum

Corneum Lipid Pathway: Results from Finite Element Models. J. Pharm. Sci. 2003, 92, 2196-2207. 19. Hansen, S.; Henning, A.; Naegel, A.; Heisig, M.; Wittum, G.; Neumann, D; Kostka, K. H.; Zbytovska, J.; Lehr, C. M.; Schaefer, U. F. In-Silico Model of Skin Penetration Based on Experimentally Determined Input Parameters. Part I: Experimental Determination of Partition and Diffusion Coefficients. Eur. J. Pharm. Biopharm. 2008, 68, 352-367. 20. Rim, J. E.; Pinsky, P. M.; van Osdol, W. W. Using the Method of Homogenization to Calculate the Effective Diffusivity of the Stratum Corneum. J. Membr. Sci. 2007, 293, 174-182. 21. 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. 22. Hofsäß, C.; Lindahl, E.; Edholm, O. Molecular Dynamics Simulations of Phospholipid Bilayers with Cholesterol. Biophys. J. 2003, 84, 2192-2206. 23. 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. 24. 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 25. Höltje, M.; Förster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Höltje, H.D. Molecular Dynamics Simulations of Stratum Corneum Lipid Models: Fatty Acids and Cholesterol. Biochim. Biophys. Acta, Biomembr. 2001, 1511, 156-167. 26. Berendsen, H.; Postma, J.; Gunsteren, W.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermol. Forces. 1981, 331−342. 27. 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.



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

28. 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. 29. 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. 30. J. C., Kasson, J.C.; P. M., van der Spoel, D., Hess, B., Lindahl, E. GROMACS 4.5: A High- throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 845–854. 31. Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. 32. 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. 33. Gupta, R.; Rai, B. Penetration of Gold Nanoparticles through Human Skin: Unraveling Its Mechanisms at the Molecular Scale. J. Phys. Chem. B. 2016, 120, 7133-7142. 34. Marrink, S. J.; Berendsen, H. J. Simulation of Water Transport through a Lipid Membrane. J. Phys. Chem. 1994, 98, 4155-4168. 35. Marrink S. J.; Berendsen, H. J. Permeation Process of Small Molecules across Lipid Membranes Studied by Molecular Dynamics Simulations. J. Phys. Chem. 1996, 100, 16729-16738. 36. Gupta, R.; Kashyap, N.; Rai, B. Transdermal Cellular Membrane Penetration of Proteins with Gold Nanoparticles: A Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2017, 19, 7537-7545. 37. Gupta, R.; Rai, B. Molecular Dynamics Simulation Study of Translocation of Fullerene C60 through Skin Bilayer: Effect of Concentration over Barrier Properties. Nanoscale.

2017, 9, 4114-4127 38. Gupta, R.; Rai, B. Effect of Size and Surface Charge of Gold Nanoparticles on their Skin Permeability: A Molecular Dynamics Study. Sci. Rep. 2017, 7, No. 45292. 39. Johnson, M. E.; Blankschtein, D.; Langer, R. Evaluation of Solute Permeation through the Stratum Corneum: Lateral Bilayer Diffusion as the Primary Transport Mechanism. J. Pharm. Sci. 1997, 86, 1162-1172. 40. Dias, M.; Farinha, A.; Faustino, E.; Hadgraft, J.; Pais, J.; Toscano, C. Topical Delivery of Caffeine from Some Commercial Formulations. Int. J. Pharm. 1999, 182, 41-47.


Page 18 of 19

Page 19 of 19

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


41. Potts, R. O.; Guy, R. H. Predicting Skin Permeability. Pharm. Res. 1992, 9, 663-669. 42. Uchiyama, M.; Oguri, M.; Mojumdar, E. H.; Gooris, G. S.; Bouwstra, J. A. Free fatty acids chain length distribution affects the permeability of skin lipid model membranes. Biochim. Biophys. Acta, Biomembr. 2016, 1858, 2050-2059.

TOC


In-Silico Skin Model: A Multiscale Simulation Study ... - ACS Publications

Recommend Documents