Introduction

Throughout the past three decades, the transdermal patch has become a proven technology that offers significant clinical benefits over other dosage forms. Compared to the more conventional oral dosing, systemic delivery of drugs by percutaneous permeation offers several advantages such as avoidance of pre-systemic metabolism in the gastrointestinal tract and liver, a controlled drug release rate that can provide relatively constant concentrations of drug in the circulation for several days, and improved patient compliance. Although transdermal patches have been marketed for about 25 years, only around 10 different compounds, all of them with molecular weights <500 amu, are currently delivered transdermally. A major obstacle to the wider application of this delivery route is the low skin permeability of many drugs, which prevents administration of effective doses from patches of acceptable physical dimensions.

The transdermal drug delivery system consists of the delivery device (the patch) and the various layers of the skin. The goal is to deliver the drug from the patch through the intervening skin layers to the microvasculature underneath. The parts of the skin that are of main interest in transdermal delivery are the uppermost stratum corneum and underlying layers of the viable epidermis. The stratum corneum has been identified as providing the primary barrier against permeation of compounds through the skin. It consists of two phases with very different characteristics—tightly packed keratinized cells (corneocytes) surrounded by the intercellular lipid matrix. The implications of the biphasic structure of the stratum corneum on its transport properties have been explored previously.30,33 For most compounds that are soluble in oil, the primary transport pathway seems to be the intercellular lipid matrix.23

One way of increasing the permeation of drugs in the skin is by including permeation enhancers in the transdermal delivery formulation. Permeation enhancers facilitate permeation of penetrants through the skin by temporarily diminishing the barrier properties of the skin.32,33 Ideally, these materials should be pharmacologically inert, nontoxic, and nonirritating, and the skin should immediately regain its barrier properties on its removal. Various classes of permeation enhancers and their mechanisms of action have been reviewed by Sinha and Kaur.38 The mechanisms by which enhancers affect skin permeability have been described as: disruption of the intercellular bilayer lipid structure; interaction with the intracellular proteins of the stratum corneum; and improvement of partitioning of a drug into the stratum corneum.4

For development and optimization of effective transdermal formulations, it is important to understand the mechanism of drug permeation from the delivery device across the skin. During the last two decades, many researchers have studied percutaneous absorption of drugs both experimentally and theoretically. Both transport through the device and through the skin are generally considered to follow Fickian diffusion through a simple homogeneous membrane. The studies on transport of compounds through the skin may be divided into two large categories: those that deal with the skin as a whole11,25,26,34 and those that concentrate on the stratum corneum as a separate layer.11,18,27,30 Here, we take the second approach, by explicitly considering the biphasic microstructure of the stratum corneum.

Among studies that have explored the implications of the biphasic structure of a medium on its transport properties, several have estimated the lipid pathway tortuosity and the effective diffusivity of the stratum corneum,23,26,30,41 while others have focused on diffusion across artificial barrier membranes.13,16 Several were 2-D studies of the cross-sections of the stratum corneum (or membrane), estimating the macroscopic diffusivity in the trans-membrane direction from geometric arguments,23 Monte-Carlo calculations16 or by measuring the pathlength directly from micrographs.41 For the 3-D structure, the method of asymptotic expansions was used to find the effective diffusivity for both impermeable and permeable corneocytes.35,36

The thinness and fragility of the stratum corneum layer make it difficult to measure its permeability directly. In addition, experimentally measured permeabilities offer scant information regarding the different contributions to the observed values, namely the interaction of the drug with its environment in the skin, and the heterogeneous structure of the stratum corneum. Therefore, the effect of designing a drug with a different chemical moiety, or the effects of adding permeation enhancers that act on the lipid and corneocytes of the stratum corneum differently cannot be readily inferred. On the other hand, the heterogeneity of the microstructure is a challenge for numerical determination of the permeability and necessitates a multiscale approach. The method should be able to account in a consistent manner for both the diffusion of substances in their molecular-level environment and the microstructure of the stratum, and be capable of considering enhancement of the permeability by a second (or third) component. Toward this goal, a hierarchical multiscale framework of modeling the transdermal diffusion of molecules is presented in this study. We start by using molecular dynamics (MD) simulations to calculate the microscopic diffusion coefficient in the lipid bilayers of the SC. MD simulations of the drug fentanyl in a lipid/cholesterol system that is close in character to the SC lipid are performed, and the molecular diffusivities of fentanyl are calculated from the resulting trajectories. Oleic acid, a chemical enhancer, is added to the model where it affects the diffusion characteristics of the fentanyl molecules.

In the next step, a homogenization procedure is performed over an idealized unit cell of the heterogeneous SC, resulting in effective diffusion parameters. The method of asymptotic expansions is used to obtain the macroscopic diffusion tensor for the stratum corneum. The effective diffusivity tensors are calculated for the case of coupled multi-component diffusion of drugs diffusing along with permeation enhancers. These effective parameters are the macroscopic diffusion coefficients for the homogeneous medium that is “equivalent” to the heterogeneous SC.

In the final step, finite element modeling of the macroscopic diffusion problem together with the effective diffusion coefficients results in the quantification of drug flux delivered through the skin.

The next three sections describe each of the micro-, meso-, and macroscale problems. The resulting drug flux through the skin shows very reasonable agreement to experimental data.

Microscopic Diffusion of Fentanyl in a Lipid Bilayer

The model drug is fentanyl (C22H28N2O) (Fig. 1). It is a synthetic opioid 75–100 times more potent than that of morphine. It was the first analgesic available for transdermal administration, with a low molecular weight of 336.5 amu and octanol-water partition coefficient of about 800. In addition, studies of the bioavailability of fentanyl delivered transdermally found that 92% of the administered dose reached the systemic circulation, with no significant cutaneous metabolism or degradation.42 The pharmacological and pharmacokinetic properties of transdermal fentanyl along with its therapeutic efficacy have been extensively reviewed.19,21

Figure 1
figure 1

Chemical structure of fentanyl

Full size image

The stratum corneum lipid composition is highly complex; in our approach, the stratum corneum lipids are modeled by a bilayer composed of dimyristoyl-phosphatidylcholine (dmpc)-cholesterol with 50 mol% cholesterol. It has been shown experimentally that the diffusion coefficients in the dmpc-cholesterol lipids correspond well to those in the stratum corneum lipids.22 Figure 2 shows the dmpc-cholesterol bilayer with a fentanyl molecule in the center of the bilayer.

Figure 2
figure 2

Molecular dynamics unit-cell of dmpc-cholesterol bilayer with a fentanyl molecule in the center of the bilayer. The bilayer is solvated by water molecules. For clarity, the fentanyl molecule is depicted with larger atoms. Colors of atoms are: H (white), C (green), O (red), N (blue)

Full size image

Permeation enhancers are an important part of many transdermal drug delivery formulations. The enhancers diffuse along with the drug and increase the permeability of the drug by temporarily altering the skin barrier (primarily the stratum corneum). One example is oleic acid, which increases drug partitioning and drug permeation in the skin. Oleic acid molecules were added to the fentanyl-dmpc-cholesterol systems in order to study the effect of oleic acid on the lipid bilayer and the diffusivity of fentanyl.

The diffusion coefficient of fentanyl molecules within the lipid membranes can be obtained from the slope of the mean-square displacement (MSD) of fentanyl vs. time:

$$ D = \mathop { \lim }\limits_{t \to \infty } \frac{1}{2s}\frac{d}{dt}\left\langle {\left( {r\left( {t + t_{0} } \right) - r\left( {t_{0} } \right)} \right)^{2} } \right\rangle,$$
(1)

where s is the number of spatial dimensions, r(t) is the center-of-mass (COM) coordinates at time t.39 The average 〈 〉 is taken over all initial times t 0.

The transport of small molecules through lipid bilayer membranes is central to many biological processes, such as passive drug and metabolite uptake. Since the early 1990s molecular dynamics studies of diffusion of small molecules in lipid bilayers have been extensively pursued.2,5,6,29 These studies have shown that small solutes diffuse at different rates depending on their location in the bilayer, with high diffusion rates found in the center of the bilayer, where the membrane density is low, while much lower diffusion rates were found in the region of the headgroups. In order to get a more detailed description of the position-dependence of the diffusion process, a series of simulations of dipalmitoylphosphatidylcholine (dppc) and dmpc lipids systems were performed with a single fentanyl molecule and different numbers of oleic acid molecules, where the fentanyl molecule is constrained in the z-direction (normal to the bilayer) to certain depths in the bilayer, while still being free to move in the xy plane.

The lateral diffusivity of fentanyl as a function of the depth of the bilayer can be found from Eq. (1). The mean-square displacements were calculated for the center-of-mass position of each fentanyl molecule and averaged over time. In order to correct for the relative bilayer leaflet motions, the center-of-mass coordinates of each fentanyl molecule was corrected by subtracting the center-of-mass coordinates of the lipid leaflet. The time origin was shifted every 10 ps for better statistics.

Homogenization of the Biphasic Structure of the Stratum Corneum—Multicomponent Coupled Diffusion

The main geometric features of the stratum corneum are the corneocytes and the lipid bilayers. The diffusion of drug through the stratum corneum involves not only the molecular diffusion of the drug in the corneocytes and the lipid bilayers, but also depends on the details of the three dimensional arrangement of the two phases. The length scale of the heterogeneity is on the order of the size of the corneocytes—about 40 μm wide and 1 μm thick—and is much smaller than the typical diameter of several centimeters of the stratum corneum of interest in transdermal drug delivery. The addition of permeation enhancers results in a coupled multi-component diffusion system. In this section, we expand the development of our previous work35,36 by considering the homogenization of a system of equations describing coupled two-component diffusion, where the diffusivity of one component (drug) is a function of the concentration of the other component (enhancer).

Consider the domain Ω of the stratum corneum such that \( \Upomega = \Upomega^{\text{cor}} \cup \Upomega^{\text{lip}},\) where Ωcor and Ωlip are the corneocyte and lipid regions, respectively (see Fig. 3). The problem is governed by Fickian diffusion in Ω and expressed by

$$ \frac{{\partial c_{a} \left( {{\varvec{x}},t} \right)}}{\partial t} - \frac{\partial }{{\partial x_{i} }}\left( {D_{ij}^{ab} \left( {{\varvec{x}},c_{1} ,c_{2} } \right)\frac{{\partial c_{b} \left( {\varvec{x},t} \right)}}{{\partial x_{j} }}} \right) = 0,\quad {\text{in }}\Upomega $$
(2)

for the components a = 1, 2, with appropriate boundary conditions. Summation is implied over the spatial indices i, j ranging from 1 to the number of spatial dimensions (3 for 3-D or 2 for 2-D), and the component index b where b = 1, 2, and will be the case in the following as well, unless indicated otherwise. The following development is for the 3-D case, but it is a simple matter to modify it for the 2-D problem. The diffusivities D ab ij are now functions both of position and concentration, and assume different functional forms depending on whether the position x falls in the corneocytes or in the lipid matrix,

$$ D_{ij}^{ab} \left( {{\varvec{x}},c_{1} ,c_{2} } \right) = \left\{ {\begin{array}{*{20}c} {\left( {D_{ij}^{ab} } \right)^{\text{cor}} \left( {c_{1} ,c_{2} } \right)} & {{\varvec{x}} \in \Upomega^{\text{cor}} } \\ {\left( {D_{ij}^{ab} } \right)^{\text{lip}} \left( {c_{1} ,c_{2} } \right)} & {{\varvec{x}} \in \Upomega^{\text{lip}} } \\ \end{array} .} \right.$$
(3)
Figure 3
figure 3

(a) Cross-sectional view of the stratum corneum and the unit cell Y. The grey areas denote the corneocytes, while the white indicate the intercellular lipid matrix; (b) a micrograph of a cross-section of the human stratum corneum, with a unit cell corresponding to that in (a). Reproduced from Talreha et al. 41

Full size image

We assume the distribution of the corneocytes in the stratum corneum is periodic, characterized by a small parameter ε, so that the diffusivities in Eq. (2) depend on ε, and are of the form

$$ \left( {D_{ij}^{ab} } \right)^{\varepsilon } \left( {{\varvec{x}},c_{1} ,c_{2} } \right) = D_{ij}^{ab} \left( {\frac{{\varvec{x}}}{\varepsilon },c_{1} ,c_{2} } \right), $$
(4)

where the D ab ij ’s are periodic functions of a reference period R. The size of the periodic unit cell of the structure is denoted by ε > 0. The periodic model of the stratum corneum structure with a unit cell Y is shown in Fig. 3. For simplicity, we will follow the standard assumption of the corneocytes being impermeable. This has been well established in the case of lipophilic drugs, for which the main transport pathway is the intercellular lipid matrix.23 In other words, D cor ij  = 0, and Eq. (2) becomes

$$ \frac{{\partial c_{a} \left( {{\varvec{x}},t} \right)}}{\partial t} - \frac{\partial }{{\partial x_{i} }}\left( {\left( {D_{ij}^{ab} } \right)^{\varepsilon } \left( {{\varvec{x}},c_{1} ,c_{2} } \right)\frac{{\partial c_{b} \left( {{\varvec{x}},t} \right)}}{{\partial x_{j} }}} \right) = 0,\quad {\text{in }}Y \cap \Upomega^{\text{lip}} $$
(5a)
$$ - \left( {D_{ij}^{ab} } \right)^{\varepsilon } \left( {{\varvec{x}},c_{1} ,c_{2} } \right)\frac{{\partial c_{b} \left( {{\varvec{x}},t} \right)}}{{\partial x_{j} }}n_{i} = 0,\quad {\text{on }}Y \cap \Upgamma $$
(5b)

for a = 1, 2, and Γ denotes the boundary between the corneocytes and the lipid phase. n i is the ith component of the unit normal vector on Γ. However, in the case of hydrophilic compounds an extension to permeable corneocytes could be easily made, parallel to the development for the single component diffusion.36

Writing the solution as \( c_{a} \left( {{\varvec{x}},{\varvec{y}},t} \right) = c_{a}^{(0)} \left( {{\varvec{x}},{\varvec{y}},t} \right) + \varepsilon c_{a}^{(1)} \left( {{\varvec{x}},{\varvec{y}},t} \right) + \varepsilon^{2} c_{a}^{(2)} \left( {{\varvec{x}},{\varvec{y}},t} \right) + \cdots,\) and expanding D ab ij (y, c 1, c 2) using Taylor’s expansion, with (y, c (0)1 , c (0)2 ) as the center, the homogenized diffusion equation can be written as,

$$ r\frac{{\partial c_{a}^{(0)} \left( {{\varvec{x}},t} \right)}}{\partial t} - \frac{\partial }{{\partial x_{i} }}\left( {\hat{D}_{ij}^{ab} \left( {c_{1}^{(0)} ,c_{2}^{(0)} } \right)\frac{{\partial c_{b}^{(0)} \left( {{\varvec{x}},t} \right)}}{{\partial x_{j} }}} \right) = 0 $$
(6)

with effective coefficients

$$ \hat{D}_{ij}^{ab} \left( {c_{1}^{(0)} ,c_{2}^{(0)} } \right) = \frac{1}{{\int\nolimits_{Y} {dy} }}\int\limits_{{Y \cap \Upomega^{\text{lip}} }} {D_{ik}^{ab} \left( {{\varvec{y}},c_{1}^{(0)} ,c_{2}^{(0)} } \right)\left( {\delta_{jk} + \frac{{\partial \chi_{bj} \left( {\varvec{y}} \right)}}{{\partial y_{k} }}} \right)dy.} $$
(7)

for a = 1, 2, and where \( r = {{\int\nolimits_{{Y \cap \Upomega^{\text{lip}} }} {dy} } \mathord{\left/ {\vphantom {{\int\nolimits_{{Y \cap \Upomega^{\text{lip}} }} {dy} } {\int\nolimits_{Y} {dy} }}} \right. \kern-\nulldelimiterspace} {\int\nolimits_{Y} {dy} }} \) is the volume fraction available for diffusion, or the porosity of the stratum corneum (barrier membrane), and the functions χ bk ’s are Y-periodic functions with b = 1, 2 and k = 1, 2, 3 that satisfy the “cell problem”

$$ - \frac{\partial }{{\partial y_{i} }}\left( {D_{ij}^{ab} \left( {{\varvec{y}},c_{1}^{(0)} ,c_{2}^{(0)} } \right)\left( {\delta_{jk} + \frac{{\partial \chi_{bk} \left( {\varvec{y}} \right)}}{{\partial y_{j} }}} \right)} \right) = 0\quad {\text{in }}Y \cap \Upomega^{\text{lip}} $$
(8a)
$$ - D_{ij}^{ab} \left( {{\varvec{y}},c_{1}^{(0)} ,c_{2}^{(0)} } \right)\left( {\delta_{jk} + \frac{{\partial \chi_{bk} \left( {\varvec{y}} \right)}}{{\partial y_{j} }}} \right)n_{i} = 0\quad {\text{on }}Y \cap \Upgamma $$
(8b)
$$ \int\limits_{{Y \cap \Upomega^{\text{lip}} }} {\chi_{ak} dy} = 0, $$
(8c)

where δ jk is the Kronecker delta function.

We postulate a simple linear relationship for the diffusivities, where the diffusivities of the drug (component 1) are linear functions of the enhancer concentration (component 2), while the diffusivity of the enhancer itself is constant. This assumption is valid for the range of concentrations being considered, as shown by the MD results (see the section “Results”). The arrangement of lipid bilayers within the stratum corneum is not understood in full detail. In the face of this uncertainty we take the simplest approach and assume the lipid phase to be isotropic. Therefore D ab ij  = D ab δ ij , and all the cross-diffusivities are taken to be zero.

$$\begin{aligned} D_{ii}^{11} \left( {{\varvec{y}},c_{1} ,c_{2} } \right) &= \left( {D^{11} } \right)^{0} + \left( {D^{11} } \right)^{1} c_{2} ,\quad i = 1,\;2,\;3.\\ D_{ii}^{22} \left( {{\varvec{y}},c_{1} ,c_{2} } \right) &= \left( {D^{22} } \right)^{0} ,\quad i = 1,\;2,\;3.\\D_{ij}^{ab} \left( {{\varvec{y}},c_{1} ,c_{2} } \right) &= 0,\quad i \ne j\quad {\text{or }}a \ne b \end{aligned}$$
(9)

where (D 11)0, (D 11)1, (D 22)0 are constants. Both the spatial and interspecies cross-diffusivities are taken to be zero.

From Eq. (7), note that the linear relationship in the microscopic diffusivities are preserved in the effective diffusivities,

$$ \hat{D}_{ij}^{11} \left( {c_{2}^{(0)} } \right) = \left( {\left( {D^{11} } \right)^{0} + \left( {D^{11} } \right)^{1} c_{2}^{(0)} } \right)\delta_{ik} S_{jk}^{1} ,\quad \hat{D}_{ij}^{22} = \left( {D^{22} } \right)^{0} \delta_{ik} S_{jk}^{2},$$
(10)

where

$$ S_{{jk}} ^{a} = {{\int\limits_{{Y \cap \Omega ^{{lip}} }} {\left( {\delta _{{jk}} + \frac{{\partial \chi _{{aj}} \left( {\varvec{y}} \right)}}{{\partial y_{k} }}} \right)dy} } \mathord{\left/ {\vphantom {{\int\limits_{{Y \cap \Omega ^{{lip}} }} {\left( {\delta _{{jk}} + \frac{{\partial \chi _{{aj}} \left( \varvec{y} \right)}}{{\partial y_{k} }}} \right)dy} } {\int\limits_{Y} {dy} }}} \right. \kern-\nulldelimiterspace} {\int\limits_{Y} {dy.} }} $$

In essence, the effective diffusion coefficients are obtained by multiplying the molecular diffusion coefficients (D aa)0 in the lipids by the factor S a jk .

The Transdermal Diffusion of Fentanyl with Oleic Acid as Permeation Enhancer

The macroscopic model system consists of a thin layer of adhesive patch that is loaded with the drug and enhancer, the stratum corneum layer, and the underlying viable epidermis layer. The viable epidermis is considered as a single homogeneous phase. In fact, a considerable amount of evidence indicates that it is a permeable phase that functions as a viscous aqueous layer to most permeants.17 A sink boundary condition is applied to the lower surface of the viable epidermis, modeling the perfect uptake of the drug and enhancers by the microcapillaries. The diffusion domain Ω is then composed of three cylindrical layers, the dermal patch (subdomain ΩA), the stratum corneum (subdomain ΩB), and the viable epidermis (subdomain ΩC). The patch is a single-layer adhesive in which the drug fentanyl (component 1, with concentration c 1) and enhancer oleic acid (component 2, with concentration c 2) are uniformly dissolved initially. A side view of the system along with its dimensions is shown in Fig. 4. The radius of the skin is set to be large enough so that the flux of the components across the sides of the skin domain is negligible. Also, we assume that no drug or permeation enhancer transport occurs through the top or the sides of the dermal patch, therefore simple zero flux boundary conditions are specified at these boundaries.

Figure 4
figure 4

The macroscopic diffusion domain and its dimensions. The domain consists of the transdermal patch (A), the stratum corneum (B), and the epidermis (C) layers

Full size image

The system of equations describing the diffusion in the domain Ωα, α = A, B, C, is given as

$$ n\frac{{\partial c_{l}^{\alpha } }}{\partial t} - \frac{\partial }{{\partial x_{i} }}\left( {\left( {D_{ij}^{lm} } \right)^{\alpha }\; \frac{{\partial c_{m}^{\alpha } }}\,{{\partial x_{j} }}} \right) = 0\quad {\text{in }}\Upomega^{\alpha}\,\times\, ]0,T[, $$
(11a)
$$ c_{l} = 0\quad {\text{on }}\Upgamma_{g}\,\times\, ]0,T[, $$
(11b)
$$ - \left( {D_{ij}^{lm} } \right)^{\alpha }\; \frac{{\partial c_{m}^{\alpha } }}{{\partial x_{i} }}n_{i}^{\alpha } = 0\quad {\text{on }}\Upgamma_{h}^{\alpha }\,\times\, ]0,T[, $$
(11c)

where the indices l, m = 1, 2, indicate the concentration component and i = 1, … ,n sd, are the spatial indices for the number of spatial dimensions n sd. \( \left( {D_{ij}^{lm} } \right)^{\alpha}\) is the ijth component of the diffusivity tensor that relates a flux of component l occurring due to a concentration gradient of component m in material α. Γ g is the lower surface of the viable epidermis with the sink boundary condition and Γ α h denotes the side and top surfaces of the domain where no flux conditions are applied. n α i indicates the ith (spatial) component of the vector normal to the boundary.

At the interfaces between the subdomains Γ I , there is an interface flux of each chemical species which must satisfy the following conditions that impose the continuity of normal interface flux,

$$ - \left( {D_{ij}^{lm} } \right)^{\alpha }\; \frac{{\partial c_{m}^{\alpha } }}{{\partial x_{j} }}n_{i}^{\alpha } = \lambda_{l} \quad {\text{on }}\Upgamma_{I}^{\alpha } \times ]0,T[ $$
(12a)
$$ - \left( {D_{ij}^{lm} } \right)^{\beta }\; \frac{{\partial c_{m}^{\beta } }}{{\partial x_{j} }}n_{i}^{\beta } = -\lambda_{l} \quad {\text{on }}\Upgamma_{I}^{\beta } \times ]0,T[, $$
(12b)

where α and β are adjacent subdomains that share the interface Γ I : the superscript on Γ I serves as a reminder of the domain to which Γ I is considered to belong for each equation. λ l is the normal interface flux of component l, and the opposite signs on λ l indicates the fact that the flux flowing out of one subdomain must flow into the other to maintain mass conservation. Summation is implied over the chemical species index m and the spatial indices i and j.

In addition to the preceding sets of equations, a condition describing the partitioning of each species across the interface is needed. The partitioning results in a discontinuity in the concentration of each chemical species across the interface Γ I . We write

$$ \left. {\frac{{c_{m}^{\alpha } }}{{K_{m}^{\alpha } }}} \right|_{{\Upgamma_{I}^{\alpha } }} = \left. {\frac{{c_{m}^{\beta } }}{{K_{m}^{\beta } }}} \right|_{{\Upgamma_{I}^{\beta}}}, $$
(13)

for m = 1, 2, where K α m describes the equilibrium distribution of component m in material α, relative to a common reference material. Since only the ratio of \( {{K_{m}^{\alpha } } \mathord{\left/ {\vphantom {{K_{m}^{\alpha } } {K_{m}^{\beta } }}} \right. \kern-\nulldelimiterspace} {K_{m}^{\beta } }} \) is relevant, the reference material can be taken to be any arbitrary material.

The formulation of the finite element equations and matrices of Eqs. (11)–(13) is essentially the same as in Rim et al.,35 except that the diffusivities are no longer isotropic, and there is a retardation factor n in the first term. The presence of n poses no difficulty as it is simply incorporated into the mass matrix.

The parameters needed to solve the problem are the retardation factors, diffusivities, and partition coefficients of the drug and enhancer in the three subdomains. The patch and the viable epidermis are considered to be isotropic, and since the lateral diffusivities in the stratum corneum are equal, the problem can be formulated as an axisymmetric one.

Computational Methods

Molecular Dynamics

The united-atom dmpc force field was taken from the dppc force field of Berger et al. 9 for the most part, with six charge groups for the headgroup.3 The modifications are described below. For the water, we used the simple point-charge (SPC) model.7 The united-atom force field for cholesterol is taken from the work of Höltje et al. 20 and is based on the GROMACS ffgmx force field. Oleic acid has a single hydrocarbon chain; the double bonded hydrocarbon parameters along with the parameters for the –COOH group were taken from the ffgmx force field.

For the fentanyl molecule, the minimum energy structure was calculated with Gaussian98 (Gaussian Inc., Wallingford, CT, USA), using the 6-31G(d) basis. The minimum energy structure is very similar to the crystal structure of the citrate–toluene solvate of fentanyl (CSD refcode: PEPCIT10) with the only significant difference being the dihedral angle of the bond between the benzene ring and the CH3 group. The bond lengths, bond angles, and dihedral angles of the fentanyl minimum energy structure and those of the crystal structure are for the most part very close to each other, and also to the corresponding parameter values of the ffgmx force field. For those bond lengths, bond angles, or the dihedral angles of the fentanyl structure that differed from the equilibrium values of the ffgmx force field by >5%, and the parameter values for the bond lengths, bond angles, or the dihedral angles that do not exist in the ffgmx force field, the equilibrium values were taken to be those from the minimum energy structure of fentanyl.

The force constants for the bond stretching, angle bending, and dihedral torsion interactions were taken from the ffgmx force field as were the Lennard–Jones parameters. For those interactions that were not defined in the ffgmx force field, we took the force constants of the interaction that most closely corresponded to the atom types involved in the interaction.

The partial atomic charges of fentanyl were calculated using the Gaussian98 program. The CHELPG charges of fentanyl were derived with the 6-31G(d) basis and based on the minimum energy configuration of the molecule. The charges for the united-atom approximation were then derived by combining the charges of hydrogen atoms into the attached carbon atom. The charges on the benzene rings were further adjusted by redistributing the charges so that the symmetry of the benzene ring is preserved.

Two sets of simulations were performed, with and without the enhancer oleic acid. The bilayers consist of 32 dmpc and 32 cholesterol molecules arranged in a bilayer form and fully hydrated with water molecules. The bilayers were equilibrated for 4 ns. For the simulations with oleic acid, 12 and 24 molecules of oleic acid were inserted at random positions in the dmpc-cholesterol systems, in equal numbers in each leaflet, and equilibrated for 4 ns each. A fentanyl molecule was inserted at the center of the dmpc-cholesterol and dmpc-cholesterol-oleic acid bilayers, and pulled in the +z and −z directions using the pull code of Gromacs 3.2.18,28 with a pull rate of 0.001 nm/ps, relative to the center of the bilayer. For every 0.2 nm, the fentanyl molecule is constrained in the z-coordinate, while the lateral displacements are collected for four simulations of 1 ns length each, for a total of 4 ns. The constraint on the fentanyl molecule is applied through resetting the z-coordinate of the center-of-mass position of the fentanyl molecule at every time step to the reference depth relative to the center-of-mass of the lipid molecules.

For the electrostatic interactions, the particle mesh Ewald (PME) method was used. The Lennard–Jones interactions were cutoff at a distance of 1 nm. The constant temperature and pressure of 310 K and 1 atm were maintained through the Berendsen weak coupling schemes, with coupling times of 0.1 and 1 ps, respectively. All bond lengths of the lipid molecules were kept constant using the LINCS routine, and the water geometry was maintained with the SETTLE algorithm. The time step for integration was taken to be 2 fs. The simulations were performed with the molecular dynamics package GROMACS version 3.2.18,28 on the Silicon Graphics Orgin 3800 shared memory supercomputer, courtesy of the Center for Biomedical Computing at Stanford University.

Homogenization

As described in the section “Homogenization of the biphasic structure of the stratum corneum—multicomponent coupled diffusion,” the effective diffusivities of the stratum corneum are found by a method of homogenization which uses the molecular diffusivities, extracted according to the methods of the section “Molecular dynamics,” in the “cell problem,” equation (8). The cell problem is solved using the finite element method as follows.

The meshes were created using MSC Patran (MSC Software, Santa Ana, CA, USA). The 2-D meshes consist of 450 to 490 rectangular 4-node elements with 4 integration points per element for numerical quadrature. The 3-D meshes have 8680 to 8950 8-node brick elements with 8 integration points per element. Separate meshes had to be created for each geometric configuration, and the element sizes were kept as consistent as possible. A finite element code that we developed was used, with a sparse gaussian elimination solver SuperLU_MT.14 To test the convergence of the meshes, a 2-D mesh was created with 1168 4-node elements. The relative difference in the results was on the order of 10−4, indicating convergence.

Finite Element Method

The mesh of the domain consists of 1230 4-node rectangular elements with 4 integration points per element, and 1442 nodes. The mesh was again created with MSC Patran. The time integration and iterative algorithms are as described in Rim et al. 35

Results

Molecular Diffusivity of Fentanyl in a Lipid Bilayer

The lateral diffusivities of fentanyl in dmpc-cholesterol and dmpc-cholesterol-oleic acid bilayers as a function of position in the bilayer are shown in Fig. 5. The lateral diffusion profiles are characteristic of solute diffusion in lipids—faster diffusion in the center of the bilayer and much slower diffusion in the dense headgroup region. The average lateral diffusivity of fentanyl in the dmpc-cholesterol bilayer is (1.2 ± 0.4) × 10−7 cm2/s. The average was taken over the inner 2 nm of the bilayer, since for lateral movement it is most likely that the molecule will stay in the tail region of the bilayer. When oleic acid is included in the system, the diffusivity of fentanyl is increased—mostly in the center of the bilayer. The average diffusivity of fentanyl was found to be (2.17 ± 0.6) × 10−7 and (3.97 ± 1.1) × 10−7 cm2/s for the systems with 12 and 24 oleic acid molecules, respectively. The oleic acid concentrations of the two systems calculated from the volume of the bilayers are 0.065 and 0.13 g/cm3. From these the following linear relationship between the oleic acid concentration and the average lateral diffusivity of fentanyl may be deduced,

$$ D_{\text{lat}} = 2.03 \times 10^{ - 6} c_{\text{oleic}} + 1.2 \times 10^{ - 7} = 1.2 \times 10^{ - 7} \left( {16.7c_{\text{oleic}} + 1} \right), $$

with R 2 = 0.965.

Figure 5
figure 5

Lateral diffusivities of fentanyl in dmpc-cholesterol (filled) and dmpc-cholesterol-oleic acid (open) bilayers as functions of depth

Full size image

The lateral diffusion of oleic acid molecules on the other hand was found to be independent of the oleic acid concentration. The oleic acid molecules have an average lateral diffusivity of (2.65 ± 0.9) × 10−8 cm2/s in the dmpc-cholesterol bilayer.

Effective Diffusivity of Fentanyl in a Homogeneous Stratum Corneum

The corneocytes are modeled as square tiles, with edge length 40 μm and thickness 0.8 μm. The horizontal and vertical gaps between the corneocytes which are occupied by the lipid matrix are taken to be 75 nm. The layers of corneocytes cells are stacked so that successive layers are displaced by exactly half the corneocyte edge length in the horizontal directions. These geometric parameters were taken from Johnson et al. 23 for air-dried (partially hydrated) stratum corneum.

For the above geometry, solving Eqs. (8a)–(8c) for the functions χ ak and using those function values to calculate the effective diffusivities result in

$$ \hat{D}_{11}^{11} \left( {c_{2}^{(0)} } \right) = \hat{D}_{22}^{11} \left( {c_{2}^{(0)} } \right) = 0.0876\left( {\left( {D^{11} } \right)^{0} + \left( {D^{11} } \right)^{1} c_{2}^{(0)} } \right),$$
$$ \hat{D}_{33}^{11} \left( {c_{2}^{(0)} } \right) = 0.0006\left( {\left( {D^{11} } \right)^{0} + \left( {D^{11} } \right)^{1} c_{2}^{(0)} } \right),$$
$$ \hat{D}_{11}^{22} = \hat{D}_{22}^{22} = 0.0876\left( {D^{22} } \right)^{0},$$
$$ \hat{D}_{33}^{22} = 0.0006\left( {D^{22} } \right)^{0},$$

and all the other diffusivity terms are zero.

Finite Element Modeling of Transdermal Diffusion of Fentanyl

The retardation factor differs from unity only in the stratum corneum, where it is equal to the porosity of 0.0891. The diffusivity of fentanyl in the patch was taken from experimental studies to be 1.0 × 10−9 cm2/s.34 Preliminary calculations indicated that the diffusivity of oleic acid in the patch does not affect the results significantly, and therefore it was assumed to be equal to the diffusivity of fentanyl. The lateral and transverse diffusivities of fentanyl without any oleic acid in the stratum corneum were found through molecular dynamics simulations and homogenization to be 1.05 × 10−8 and 7.2 × 10−11 cm2/s, respectively. The lateral and transverse diffusivities of oleic acid in the stratum corneum were similarly found to be 2.32 × 10−9 and 1.59 × 10−11 cm2/s. Oleic acid was found to increase the diffusivities of fentanyl as a linear function of oleic acid concentration, with constants of proportionality of 1.78 × 10−7 and 1.22 × 10−9 cm5/(g s), for the lateral and transverse diffusivities. The diffusivity of oleic acid was found to be unaffected by the concentration. Data on diffusivities of compounds in the viable epidermis are unfortunately very rare, and are not available for fentanyl and oleic acid. Chandrasekaran et al. 11 measured the diffusivity of scopolamine (C17H21NO4) in delipidized epidermis, and found a value of 2.0 × 10−7 cm2/s. Since scopolamine has a molecular weight of 303.4, comparable to those of fentanyl (336.5) and oleic acid (282.5), we chose to adopt the diffusivity of scopolamine for the diffusivities of fentanyl and oleic acid in the viable epidermis layer. The lateral and transverse diffusivities of fentanyl and oleic acid are tabulated in Table 1. The values of transverse diffusivities for both fentanyl and oleic acid are well within the range of experimental diffusivities for the stratum corneum of 10−9 to 10−12 cm2/s.37

Table 1 Diffusivities of fentanyl and oleic acid in the different layers of the diffusion domain
Full size table

The partitioning of drug and enhancer between the patch and the stratum corneum, and between the stratum corneum and the viable epidermis, are important determinants of the transdermal permeation process. The partitioning coefficient of fentanyl from the patch into the stratum corneum has been estimated experimentally to be 0.14.34 Both fentanyl and oleic acid are lipophilic molecules, and the same value is assumed for the partitioning of oleic acid into the stratum corneum. For partitioning between the stratum corneum and the viable epidermis, taking into account the aqueous nature of the viable epidermis, the partition coefficients of both fentanyl and oleic acid are approximated by the empirical formula given by Johnson et al.,23

$$ K_{\text{SC/aq}} = K_{\text{o/w}}^{0.76},$$

where K SC/aq is the partition coefficient between the stratum corneum and an aqueous phase, and K o/w is the octanol-water partition coefficient. From the values of the octanol-water partition coefficients for fentanyl and oleic acid, which are estimated to be 860 and 5 × 107, respectively, the partition coefficients between the stratum corneum and the epidermis are 0.006 and 1.41 × 10−6 for fentanyl and oleic acid, respectively.

The initial concentration of fentanyl in the patch was 0.09 g/cm3. The macroscopic diffusion problem was solved for three different initial concentrations of oleic acid in the patch: 0, 0.13, and 0.25 g/cm3. The concentration of fentanyl and oleic acid in the patch, the stratum corneum and the epidermis layers at time t = 10 h are shown as a function of depth (in the center of the domain) in Fig. 6. Note that the increase in fentanyl diffusivity in the stratum corneum due to the enhancer results in the concentration profiles crossing in Fig. 6b. The concentrations at the interfaces between subdomains are discontinuous due to partitioning.

Figure 6
figure 6

Concentration profiles of fentanyl in the (a) patch; (b) stratum corneum; and (c) viable epidermis as a function of depth from the top surface of the patch at t = 10 h. (diamonds = no oleic acid, squares = initial oleic acid concentration in patch = 0.13 g/cm3, triangles = initial oleic acid concentration in patch = 0.25 g/cm3)

Full size image

The flux of fentanyl crossing the lower boundary (Γ g ) of the domain is shown in Fig. 7a as functions of time. The flux of fentanyl when the stratum corneum layer is omitted is shown in Fig. 7b. The experimental fentanyl flux from a transdermal patch through human cadaver skin, courtesy of the ALZA Corporation, is shown in Fig. 7c.34 Details on the experimental setup is included in the supplementary material.

Figure 7
figure 7

(a) Fentanyl flux through the lower boundary of the viable epidermis for different initial concentration of oleic acid in the patch; (b) fentanyl flux through the lower boundary of the viable epidermis when the stratum corneum is removed; (c) experimental fentanyl flux through human skin; for details see the supplementary information

Full size image

Discussion

The lateral diffusivities of small compounds have been measured by fluorescence recovery after photobleaching (FRAP) to be in the range of 6 × 10−8 to 1.6 × 10−7 cm2/s in pure dmpc bilayers and 2–5.6 × 10−8 cm2/s in dmpc-cholesterol bilayers at 300 K.1,22 The average lateral diffusivity of fentanyl in the model lipid system of (1.2 ± 0.4) × 10−7 cm2/s is on the same order as the experimental data, the difference may be attributed to some or all of the following: differences in the liphophilicity and shape factor of the molecules considered and the concentration of the molecules in the bilayer.

The effect of oleic acid on the bilayer is to increase both its cross-sectional area and the disorder of the lipid tails. These may largely be due to the kink in the unsaturated hydrocarbon tail of oleic acid that needs to be accommodated, consistent with the observation that unsaturated fatty acids are more effective in enhancing transdermal permeation of drugs than saturated fatty acids.12 These changes in the bilayer structure act to increase the diffusivity of the fentanyl molecule, proportional to the amount of oleic acid present in the bilayer. The oleic acid itself has a much smaller diffusivity, in agreement with the observation that oleic acid, once introduced, remains in the skin for a prolonged period of time. The oleic acid concentration does not seem to affect its own diffusivity, though this may be due to the small concentrations of oleic acid present in the simulations. The lateral diffusivity of oleic acid in the bilayer is on the same order as that of the lipid molecules, which implies that oleic acid, which is a fatty acid, is easily incorporated into the lipid bilayer and effectively becomes a component of the bilayer itself.

In this work, we assumed that the main mechanism for permeation through the lipid bilayers in the stratum corneum is lateral diffusion. This is contrary to considerations of drug molecules crossing the plasma membrane, where the trans-bilayer diffusion is the primary transport mechanism. One reason for this is that lateral diffusion is faster than transbilayer diffusion by a factor of up to nine orders of magnitude.23 In addition, although the organization of lipid bilayers in the stratum corneum is not well known, a reasonable assumption would be that the lipid bilayers are continuous throughout the stratum corneum, with defects and interconnections between the individual bilayers that allow solute to traverse the entire stratum corneum by lateral diffusion alone.23,43 In fact, electron microscope (EM) studies of the stratum corneum have shown that adjacent bilayers often merge together, such that solutes can diffuse from one bilayer to another by lateral diffusion.15,40 Note that in the limit of a large enough number of interconnections between the bilayers, the lipid can justifiably be considered isotropic as in (10).

For molecules for which the corneocytes are impermeable, the effective diffusivity in the stratum corneum is a measure of the tortuosity of the diffusion path. The large difference of two orders of magnitude between the in-plane and transverse effective diffusivities reflects the different tortuosities of the diffusion path in each direction. The horizontally elongated shape of the corneocytes acts to hinder the diffusion mainly in the trans-stratum corneum direction, and therefore effectively contributes to the permeation barrier properties of the skin.

The stratum corneum with its high resistance to diffusion acts as the rate-regulating membrane in the transdermal system. The small partitioning coefficient of fentanyl between the stratum corneum and viable epidermis also plays an important role in slowing down the rate of permeation. The small diffusivity of fentanyl in the stratum corneum acts to keep the flux at a fairly constant level for longer periods, which is desirable in transdermal delivery systems. The role of the stratum corneum as the rate-regulating layer can be seen more clearly in Fig. 7b. Without the stratum corneum layer, the flux of fentanyl shows an initial sharp peak followed by a rapid decrease. Although the stratum corneum can easily be removed, the result of removing the permeation barrier is not simply that of increased drug flux which may be beneficial, but also a sharp variation in the flux with time that is undesirable.

The peak fentanyl flux of about 1.3 μg/(cm2 s) through the diffusion domain is smaller than, but comparable to the experimental value of about 2.5 μg/(cm2 s) as in Fig. 7c.34 The larger experimental values can be explained by the fact that the skin samples were stored fully hydrated before the experiments were conducted. The hydration level of the skin and subsequent swelling of the corneocytes was shown to result in significant increases in the effective diffusivity through the stratum corneum.35 In particular, experimental studies both in vitro and in vivo suggest that hydration can increase flux up to twofold. In view of this fact, the agreement between the calculations and the experiments seems quite promising, especially considering that many approximations and assumptions were made over the different length scales. Another factor to be addressed is the discrepancy in the flux-time profile between the numerical and experimental data. This suggests that the parameters values are not correct, in view of the fact that the calculated flux peaks earlier and remains much steadier than the experimental data indicates that the diffusion coefficient might be too large, and partition coefficient too low. In addition, some of the parameters such as the partitioning coefficients would be time dependent due to a change in environment. Systematic parametric studies would cast light upon these questions; further exploration of these issues are beyond the scope of the current work.

The addition of oleic acid increases the flux of the fentanyl appreciably. In addition, the small diffusion coefficient of oleic acid in the stratum corneum leads to the change in shape of the flux curves with the addition of oleic acid in Fig. 7a. The flux reaches its maximum value at a later time when oleic acid is added to the system. This is due to the slow diffusion of the oleic acid in the stratum corneum; the enhancement effect increases with time as the oleic acid lingers and builds up in the stratum corneum. This changes the “steadiness” of drug delivery: the length of time where the flux is within 10% of the peak value decreases from 40 to 30 to 27 h as the oleic acid concentration is increased from 0 to 0.25 to 0.5 g/cm3. The addition of oleic acid increases the rate of delivery, but at the cost of decreasing the effective time period of each patch.

A main drawback of the method is the dependence on experimental results, which may not be readily available for other molecules. Another limitation is the computational cost of the molecular dynamics simulations—however, the advance of coarse-grained membrane lipid models10 means that this can be significantly reduced.

Some uncertainties regarding the effects of permeation enhancers have not been addressed here. Comparisons of stratum corneum extracted lipids with dmpc and dmpc-cholesterol bilayer systems have suggested that liquid crystalline phase dmpc-cholesterol bilayers would be a reasonable model of stratum corneum lipids for diffusion coefficients of solute molecules22,24,31 and is further corroborated by the agreement for the calculated and experimental fentanyl flux. However, the effects of permeation enhancers on dmpc/cholesterol bilayer systems might be different from their effects on stratum corneum lipids, and experimental verification of the degree to which oleic acid increases fentanyl flux is needed.

The two major effects of adding a permeation enhancer into a transdermal delivery formulation is an increase in both the diffusion coefficient of the drug in the skin and the solubility (and therefore the partition coefficient) of the drug in the skin. Synergy between the two effects will have a multiplicative effect on the transdermal flux of the drug. However, few enhancers have been studied in sufficient detail to distinguish between these effects, and in this work, we have only considered the effect on the diffusion coefficient. In our earlier work,34 we presented a framework for posing and better understanding the question of the relative magnitudes of these two effects, by modeling the increase in diffusivity and partitioning independently as functions of the enhancer concentration, this may be easily incorporated into the current framework to explore the effects of partition coefficient more thoroughly.

Conclusions

This study presents a framework for studying the multi-component diffusion of drugs coupled with enhancers through the skin by considering the microstructure of the stratum corneum (SC). The different length scales involved in the transdermal diffusion process are addressed in detail and the contribution of each length scale is integrated into the macroscopic permeation barrier of the stratum corneum. This is a significant improvement from the previous pharmacokinetic models or purely macroscopic continuum studies of the transdermal diffusion problem. This framework is applicable to general problems of transport through heterogeneous composite media.

Areas for further investigation would include improved and more efficient lipid models, better qualification of the effects of different classes of permeation enhancers on model lipid systems, study of the corneocytes, and a more realistic geometry of corneocytes. Experimental information on partitioning and diffusion properties of the lipid and corneocytes are also very much in need.