Finite Element Analysis for Predicting Skin Pharmacokinetics of Nano Transdermal Drug Delivery System Based on the Multilayer Geometry Model

Introduction

Transdermal drug delivery system (TDDS) is of interest for the drug delivery system that can be subject to a strong first-pass effect, or the skin-targeted drug.¹ Skin pharmacokinetics, the monitor of detecting the fate of TDDS in vivo, can directly reflect the drug concentration change, the intensity of the drug at the target site, the relationship of dose-effect, and drug absorption and metabolism in the skin. However, there is little research on the TDDS skin pharmacokinetics throughout the whole skin.

SC is considered to be the main barrier to percutaneous absorption of drug. So, the dermatopharmacokinetics, referred to as the pharmacokinetics of TDDS in SC, was widely applied to investigate the drug metabolic in skin.² However, the heterogeneous structure of skin can be arranged into stratum corneum (SC), viable epidermis, dermis, and subcutaneous tissue from a macro perspective: SC, the brick-and-mortar model, is the formidable barrier of drug absorption through the skin; viable epidermis, structured with viable cells, inherently presents certain obstacles to TDDS diffusion; while the dermis surely effects TDDS on penetration depth and extent underlying the skin.³ The penetration process of TDDS includes overcoming the SC barrier, penetrating into the viable epidermis and dermis, and then being absorbed by capillaries. So, the dermatopharmacokinetics might not fully explain the skin pharmacokinetics of the drug. In addition, TDDS pharmacokinetics was also determined by measuring the drug concentration over time in blood⁴ or the drug retention in skin tissue after subcutaneous administration,⁵ or a combination of the two.⁶ However, the blood pharmacokinetics did not accurately reflect the drug dynamic distribution in the whole skin and the drug retention in skin could not insight the drug transport properties in different skin layers. To more accurately explore the dynamic distribution of TDDS in various layers of the skin, drug transport properties in skin different layers require a more systematic study, especially for the drugs to target the specific skin layer.

Mathematical models of skin have lead to increasing concern in predicting the permeability of transdermal drug delivery. The predictive mathematical model for skin penetration was developed from the steady-state models (quantitative structure–permeation relationship models, structure-based models, and porous pathway model) to the transient models with time dependence (including basic models, compartment models, complex models, and slow binding/partitioning kinetics in the SC), in which the compartment models, also known as PK models, could be applied to trace the drug fate after penetration into skin.⁷ McCarley et al⁸ described one/two-compartment models representing the SC/viable epidermis to predict the drug absorption into and through the skin. Unfortunately, the one/two-compartment mode could not fully accurately reflect the diffusion process of the drug in the whole skin. (Previously, we also used a two-compartment model to study the transdermal process of the drug. However, the in vivo skin pharmacokinetics results were significantly different from the simulation results (P<0.01) as shown in the Supplementary data). So, a multi-compartment model should be established to simulate the entire skin and drug penetration process. However, drug permeability in different skin layers cannot be obtained by experiment due to the complex skin structure and limited experimental facilities. In recent years, a mathematical method was reported to compute drug penetration processes in time and space.⁹

The FEA provides an approximate numerical solution of a partial differential equation. The primary basis for the FEA is handle domains with polygon meshes (elements) and boundaries, which could discrete the continuous domain. For example, skin could be dispersed into connecting subdomains. As reported, the two-dimensional (2D) FEA model was used to investigate the diffusion of lipophilic solute in SC.¹⁰ Rim et al¹¹ developed a FEA model to mimic the percutaneous diffusion of compounds into two isotropic materials, respecting vehicle and skin. It was also reported that 2D FEA was used to simulate the dynamic water diffusion in SC and macromolecules penetration in the epidermis.^12,13 However, the simulation of the drug diffusion over time in skin was confined to studying that in SC, or the skin described as homogeneous, as well as the simulation results are rarely experimentally verified. To probe drug absorption and metabolism over time and space (three-dimensional, 3D) in whole skin structured with different layer, the multilayer geometry of TDDS, SC, viable epidermis, and dermis should be meshed with different polygons and density in the FEA. And, the simulation of FEA was executed with the input Kps, including the Kp of drug permeability from vehicle into SC, from SC to the viable epidermis, and from viable epidermis to dermis.

Paeonol (PAE, 2ʹ-hydroxy-4ʹ-methoxyphenyl), extracted from the traditional Chinese Medicine Paeonia moutan and Radix Cynanchi Paniculati, is widely used as an anti-inflammatory and anti-allergic.¹⁴ In addition, PAE loaded in the nanoemulsion could improve the properties of drug loading and stability.¹⁵ Therefore, PAE-NEs was selected as a model nano transdermal drug to study the skin pharmacokinetics. In this study, we established a multilayer geometry model, and the Kps of PAE-NEs in SC, viable epidermis, and dermis was obtained by MD simulation, in vitro permeation experiments, and in vivo skin microdialysis probe technology. Besides, the FEA method was applied to predict the fate of PAE-NEs in the skin according to the input parameters. Finally, the feasibility of simulating drug skin pharmacokinetics using FEA was verified by in vivo experiments. The developed combination of multilayer geometry model, numerical simulation, and experiment is a powerful tool to push the development of TDDS skin pharmacokinetics.

Materials and Methods

Materials

Lecithin E200 was purchased from Degussa (Germany). Isopropyl myristate (IPM) was obtained from Sinopyharm Chemical Reagent Co., Ltd. Alkyl poly-glucosides (APG) was purchased from China Research Institute of Daily Chemical Industry. Paeonol, with a purity of not less than 98.0%, was supplied by the National Institutes for Food and Drug Control (Beijing, China). Methanol was obtained from Sigma-Aldrich (HPLC grade). Other reagents were of AR grade.

Animals

Male SD rats of SPF grade, weighing 200±20 g and about 1.5 months of age, were purchased from the Second Military Medical University. All animal protocols complied with the International Ethical Guideline and National Institutes of Health Guidelines on the Care and Use of Laboratory Animals, and with the approval of the Committee on Ethics of Biomedical, Second Military Medical University; the approval number is PREC2017-073.

Preparations of PAE-NEs

PAE-NEs were prepared using emulsification technology. Briefly, The PAE (1%) dissolved in IPM (12%) was mixed with alkylpolyglucosides (APG, surfactant, 12%), lecithin (surfactant, 6%), and 1, 2-propylene glycol (co-surfactant, 9%). Then, the mixture was added with aqueous (60%) dropwise at room temperature under stirring at 300 revs/minute until the system was clear.

Establishment of the Model of the Drug Transport System

The drug transport system was composed of the TDDS and the skin (Figure 1). Due to the heterogeneous skin, the drug permeability in the skin depends on the depth of penetration. The fate of drug exposed to skin was diffusion into SC, penetration or metabolism into the viable epidermis, dermis, and absorbed by the system or binding to tissues.¹⁶ Hence, the multilayer structure model of the intact drug transport system was divided into TDDS and the skin structured with SC, viable epidermis, and dermis, as shown in Figure 1A and B.

Figure 1 The multilayer geometry model of the drug transport system: vehicle, SC, viable epidermis, and dermis. (A) The 3D schematic diagram of vehicle and skin texture. (B) The 2D schematic diagram of drug permeation process from vehicle to subcutis and the highlighted pathways remarked with diffusion coefficients (Dm, Ds, De, and Dd).

The drug permeability profile of the drug transport system as a function of time (t) and depth (x) in the vehicle, skin, and capillary via the Fick’s diffusion equation are as the following equations (Eqns. 1–4).^17,20 And the initial conditions are as follows (Eqns. 5–8):

(1)

(2)

(3)

(4)

where C_m, C_s, C_e, C_d; D_m, D_s, D_e, D_d; and L_m, L_s, L_e, L_d stand for the drug concentration; drug Kp; and thickness of vehicle, SC, epidermis, and dermis, respectively.

(5)

(6)

(7)

(8)

Where C_m0, C_s0, C_e0, and C_d0 stand for the initial drug concentration in the vehicle, SC, viable epidermis, and dermis of the first application, respectively.

For the boundary condition, there is no exchange of drug between the drug delivery system and the surrounding (Eqn. 9). Eqn. (10) shows the equilibrium condition at the vehicle/skin surface (SC). Eqns. (11–13) represent the continuity of flux across the vehicle/SC, SC/viable epidermis, and viable epidermis/dermis, respectively. Moreover, Eqn. (14) states the drug absorbed and eliminated by the capillary.

(9)

(10)

(11)

(12)

(13)

(14)

where C_m, C_s, C_e, and C_d mean the drug concentration in vehicle, in SC, in viable epidermis, and in dermis, respectively, K_m means the isolation factor at the interface between the vehicle and skin, and k_cl refers to the clearance rate of drugs by capillaries.

The Kps of Drug in the Multilayer Geometry Skin Model

The parameters of drug Kps in the different layers of the multilayer structure model are the key points for the numerical simulation. In this work, the Kps in different skin layers were conducted by different methods, including the MD simulation, in vitro penetration test, and in vivo skin microdialysis probe technology. The drug Kp was calculated based on Fick’s law, as in Eqn. (15):

(15)

where D (cm²/h) is the drug Kp, J (μg/cm²/h) is the permeation flux in each layer, C (μg/cm³) is the initial interface drug concentration, and x (cm) is the thickness of model layers.

MD Simulation for PAE Kp in SC

The technology of coarse-grained (CG) MD simulation was applied to simulate the process of the drug diffusion in skin SC. In the simulation, the dimyristoyl-sn-glycero-3-phosphocholine (DMPC) was structured to simulate the bi-monolayer system of SC using the software of GROMACS V4.5.1 and van der Waals (VMD).^21,23 Besides, the TIP3P were constructed to simulate the water in the initial system.^24,25

The parameters of CG DMPC and water force field were set according to the theory of Marrink et al.²⁶ The simulation environment is room temperature and atmospheric pressure. The atomic bond length is kept within the equilibrium distance under the LINCS algorithm, and the electrostatic interaction is calculated using Particle-Mesh-Ewald (PME) algorithm. In addition, the descent algorithm and the conjugate gradient algorithms were used to minimize the system energy. Then drug molecules were placed on the upper surface of the phospholipid bilayer to simulate the process of membrane permeation and the dynamics diffusion under 310.15 K for 100 ns. The Kp of the drug molecule through SC was calculated according to the Einstein relation:

(16)

where D represents the drug Kp, r(t) represents the proton coordinates at time t, and represents the mean square displacement (MSD) calculated from the initial time t₀.

The Kp of PAE-NEs in the Viable Epidermis

Franz diffusion cells were employed to investigate the Kp of PAE-NEs in the viable epidermis. Briefly, rats abdominal hair was removed with electrical clippers and the skin harvested without subcutaneous tissue and then the SC was stripped with tape stripping.² The skin sample was treated with 1 M NaBr solution for 4 hours, and the viable epidermis was obtained with a medical cotton swab.²⁷ Then the tissue of the epidermis layer was mounted between the donor and acceptor compartments and treated with 0.5 g of PAE-NEs. The acceptor compartments were filled with receptor medium (normal saline) and stirred under 300 revs/minute at 32±0.5°C.^28,29 At regular intervals (0.5, 1, 2, 4, 6, 8, 10, and 12 hours), 1 mL of receiver medium was withdrawn for analysis and the acceptor compartments refilled with fresh medium. The cumulative permeability was detected and the Kp was fitted using the penetration profile and Eqns. (15), (17), and (18):^29,30

(17)

(18)

where Q_n (μg/cm²) is in vitro accumulation percutaneous permeability, C_n is the drug concentration of extracted samples, V₀ and V_i are the volume of receiver medium and extracted samples, C_i is the drug concentration of the ith extracted samples, S is the effective penetration area, and J (μg/cm²/h) is the slop of PAE-NEs penetration profile.

The Kp of PAE-NEs in the Dermis

The Kp of PAE-NEs in the skin dermis was obtained using skin microdialysis probe technology combined with MATLAB software. Rats anesthetized with 25% of urethane (0.4 mL/100 g) were fixed on warm-pad supine to maintain body temperature. Then, a skin microdialysis probe was implanted in the dermis guiding with an 18 G puncture needle under monitoring by 3D computer tomography (CT). Then, the probe implanted in the abdominal skin was balanced for 1 hour with perfusate (PBS solution) at a flow rate of 5 μL/min.²⁹ After the rats were treated with PAE-NEs (area administration=3×4 cm²), the dialysis solution was sampled and analyzed using LC-MS every 20 minutes (100 μL) for 12 hours. The drug concentration in the dermis was calculated according to the established microdialysis methodology, as displayed in the supplementary data. Then, the concentration–time curve of drug in dermis was Fourier series fitted using the software of MATLAB. Further, the concentration-time curve was integrated to obtain the total drug flux. And the Kp of PAE-NEs in the dermis was calculated using Fick’s law (Eqn. 15).

FEA for Simulation of the PAE-NEs Skin Pharmaceutics

Numeralization of the Multilayer Geometry Models

Mapped mesh is used for numeralization of the multilayers of the drug transport system with the thickness of the vehicle, SC, viable epidermis, and dermis at 0.1 cm, 0.02 cm, 0.02 cm, and 0.16 cm, respectively.^31,32 The element grids were divided and optimized using the grid division tool provided by PDE tool software in Matlab 12.0. The initial drug concentration in the simulation was 0.2 g/cm³, and that in the skin was 0 g/cm³. The dimension of the transverse (spread direction) is much larger than that of the longitudinal (depth direction) in entire skin layers, which meet the boundary conditions of the numerical multilayer geometry model. The time step and convergence criteria were set at 1 and 0.0001, respectively.

2D and 3D Dynamic Penetration Processes of PAE-NEs in Skin Using FEA

The multilayer geometry model and the input parameters (Kps of the drug in SC, viable epidermis, and dermis detected by MD simulation, in vitro permeation test, and in vivo microdialysis probe technology, respectively) were used to simulate the transport process of the PAE-NEs in skin. 2D dynamic diffusion processes of PAE-NEs over a period of time within the numeralization multilayers model and the 3D drug concentration in the entire skin over time and space was profiled using the FEA method assisted with the PDE Tool in Matlab12.0.

Validation of FEA Method Used for Predicting Skin Pharmacokinetics of TDDS

The experimental and simulated drug concentration–time curve results were compared to validate the feasibility of the FEA method to predict skin pharmacokinetics of TDDS. For experiment, the drug concentration over time in the dermis layer was acquired using the skin microdialysis probe technology and the depth of the probe in the skin was monitored by a CT. Correspondingly, a single line (the profile of drug concentration–time) of the simulation curve was separated from the 3D dynamic processes of PAE-NEs obtained using FEA, in which the depth was the same as the probe in the dermis.

Statistical Analysis

The data was shown as mean±SD. The PK parameters were analyzed with Kinetica 5.0 software. P-values less than 5% were considered to be significant.

Results and Discussion

The Kps of Drug in the Multilayer Geometry Skin Model

MD Simulation for PAE Kp in SC

Research on the drug permeation process in SC is restricted with the traditional method, whereas MD simulation is a powerful tool to visually investigate the drug transport process over time and skin depth.^33,34 In the MD simulation, coarse-grained (CG) modeling allowing long-time period and larger length scale simulation compared with the traditional atomistic models.^22,35 Considering the main penetration barrier in the SC is located in lipid and the phospholipid molecule of DMPC is considered to be similar to the lipid membrane structure, DMPC was selected to construct the SC phospholipid bilayer as shown in Figure 2A and B.^22,36 Further, the simulation system structured by DMPC, water (Grey), and PAE (blue) is displayed in Figure 2C. As shown in Figure 2D–F, the dynamic diffusion process of PAE in the skin indicated that the drug gradually approached the DMPC as the simulation time increases.

Figure 2 MD simulation for PAE dynamic penetration in SC. (A) 2D schematic diagram of phospholipid bilayer structure of SC; (B) the 3D structure of phospholipid molecules modeled by van der Waals spherical model. Carbon atoms: dark gray; oxygen atoms: red; nitrogen atoms: blue; hydrogen atoms: light gray; and phosphorus atoms: gold; (C) simulation system, the hydrophobic part of DMPC is modeled with the van der Waals model: Blue, red, and gold marked nitrogen atoms, Oxygen and phosphorus atoms, respectively. While the hydrophilic tail carbon chain is modeled with the current model and the carbon atoms are gray. In addition, both sides of the bimolecular shaded in gray are water molecules; (D–F) the dynamic processes of PAE in SC, the coordinate of (D–F) were (35.758, 29.998, 23.602), (40.145, 25.609, 24.273), and (33.722, 21.450, 31.305), respectively; (G) changes in phospholipid bilayer thickness; (H) the diffusivities trajectories of the drug molecules (n=6).

As shown in Figure 2G, the thickness of the phospholipid bilayer was kept within a reasonable and stable range within 100 ns of the simulation, indicating reasonable parameters, settings, and a stable simulation system. In the analysis of drug diffusion trajectories (n=6, Figure 2H), the stochastic initial conditions (the relative position of the drug and DMPC) were contributed to the incomplete match of drug molecular trajectories. In addition, the mean square displacement (MSD) (Figure 2H) combined with Eqn. (16) was used to fit the Kp of PAE in SC, calculating the Kp was approximately 2.00×10⁻⁶ cm²/h. However, the MD simulation of PAE-NEs penetration in SC or that in the epidermis and dermis of the skin requires further research.

The Kp of PAE-NEs in the Viable Epidermis

The permeation profile of PAE-NEs in the viable epidermis is illustrated in Figure 3. According to Eqns. (17) and (18) (Q_n=103.03 t+66.86, R²=0.996), Eqn. (15) and the viable epidermis thickness (0.20 mm),³⁶ the Kp of PAE-NEs in the viable epidermis was 1.58×10⁻⁵ cm²/h. Besides, the cumulative permeability reached 75.82±2.85% at 12 hours. In general, tape-stripping is a minimally invasive approach to obtain the SC layer for measuring the drug concentration in SC.² In the research, the method of tape-stripping was applied to remove the SC. The in vitro penetration test was a wide alternative method and a useful tool to investigate the TDDS pharmacokinetics.^37,38 So, the in vitro permeability of PAE-NEs in the viable epidermis could be calculated according to the in vitro penetration profiles. The Kp of PAE-NEs in the viable epidermis is slightly different from previous reports, which might be due to the fact that the Kp of drug varies with the drug molecule and the skin species.^13,39

Figure 3 In vitro permeation profiles of PAE-NEs. The solid lines in red represent accumulation percutaneous permeability and the solid lines in blue represent the accumulation permeability (n=5).

The Kp of PAE-NEs in the Dermis

Microdialysis probe technology used to detect the TDDS pharmacokinetics was frequently reported.^29,40,41 In the research, we creatively applied the technology to detect the Kp of drug in the dermis. The depth of the microdialysis probe implanted in the dermis visualized and located by the CT was 1.36 mm. A Fourier series curve of the PAE-NEs concentration–time in the dermis was fitted with MATLAB, as illustrated in Figure 4. According to Fick’s law, the Kp of PAE-NEs was 3.20×10⁻⁵ cm²/h. The results of drug Kp in SC, viable epidermis, and dermis demonstrated that the penetrability of PAE-NEs was gradually increased with skin depth and that in the dermis is about twice that in the viable epidermis, which is consistent with previous reports.^13,39,42

Figure 4 Concentration–time profile of PAE-NEs in microdialysate of skin probe located at 0.136 cm in the skin (solid points) and the Fourier series fitting curve of the drug concentration over time (solid line).

FEA for Simulation of the PAE-NEs Skin Pharmaceutics

Numeralization of the Multilayer Geometry Models

The FEA is a powerful grid-based method for the numerical solution of partial differential equations for dynamic drug penetration within the skin.^13,23 The spatial domain interconnected by node was divided into elementary shapes, which can be triangles, quadrangle, and polygons. The nodes density varies in the domains and each point and node have a unique property.⁴³ The meshing of the multilayer domain, the diffusivities, as well as the internal boundary conditions of the skin multilayer models are shown in Figure 5. The time domain is 10 hours, the time step = 0.001 hours = 0.06 minutes and the convergence was 0.001.

Figure 5 Finite elementary mesh of the finite element solutions (A) and the multilayer geometry of drug transport system (B), Dm, Ds, and De are the Kp of drug from vehicle to SC, from SC to viable epidermis, and from viable epidermis to dermis, respectively.

2D and 3D Dynamic Penetration Processes of PAE-NEs in Skin Using FEA

The 2D dynamic penetration process of PAE-NEs at different simulation times in the drug transport system is shown in Figure 6A–D. The results indicated that the drug penetration is time-dependent: initially, there is no penetration of the drug at the simulation time=0. And as the simulation time prolongs, the drug penetration proceeds from the vehicle to the SC, viable epidermis, and dermis. Besides, as depicted in Figure 6E, the 3D schematic of the drug concentration of PAE-NEs over time and skin depth based on the four-layer model simulates the entire dynamic process from drug administration to metabolism in the skin.

Figure 6 The drug concentration of PAE-NEs over time and skin depth. (A–D) The 2D drug penetration process and the drug concentration at the FEA simulation time of 0, 100, 200, and 300 min, respectively; (E) 3D schematic of the drug concentration–time–skin depth of FEA simulation.

Validation of FEA Method Used for Predicting Skin Pharmacokinetics of TDDS

In the validation experiment, the depth of the liner microdialysis probe implanted was 1.68 mm. The drug concentration–time profile is shown as a green dotted line in Figure 7. The comparison between the drug concentration–time curves at the skin depth of 1.68 mm measured by in vivo skin microdialysis probe technology and the numerical simulation profiles extracted from FEA simulation profile of drug concentration-time-skin at a depth of 1.68 cm were plotted in the same coordinate system as shown in Figure 7: the green dotted line and blue dotted line stand for the experimental results and the numerical simulation results, respectively. The comparative results showed that the overall skin pharmacokinetics could best be fitted by the FEA. The small error between the experimental and simulation results indicated that the FEA method combined with a multilayer geometry model is an accurate method for investigating the skin pharmacokinetics of PAE-NE.

Figure 7 The comparison between the experimental data obtained by skin pharmacokinetics and FEA simulation results. The green dotted line stands for the experimental results, while the blue dotted line refers to the numerical simulation results.

Study Limitations

In the research of the FEA method used for simulating skin pharmacokinetics of TDDS based on the drug Kps in a multilayer geometry model, the limitation is whether it is possible to accurately determine the Kps of the drug in the SC, viable epidermis, and dermis. The Kp of the PAE molecule in SC simulated with MD simulation was an approximation. And the Kp of preparations (PAE-NEs) using the MD simulation model needs further research with the development of computing power. Besides, the permeability of PAE-NEs in the viable epidermis using exfoliated skin might be different from that in vivo. Indeed, the skin inevitably undergoes elastic deformation in the preparing process, which may affect transcutaneous penetration.^13,40 As for PAE-NEs permeability in dermis, the flux J was fitted with Fourier series fitting, which cannot be completely matched with the experimental concentration–time curves. It is important to realize that the more complex models and technologies should be developed to obtain more accurate input parameters for predicting skin pharmacokinetics with FEA.

Conclusions

This was the first study to combine FEA simulation methods based on the mathematical multilayer geometry of drug transport system to investigate skin pharmacokinetics of TDDS. We explored the interaction between PAE and SC on a molecule level using coarse-grained MD simulation to simulate the Kp in SC. Also, an in vitro penetration experiment of PAE-NEs in extraneous skin was carried out to gain the Kp in the viable epidermis. In addition, The Kp of PAE-NEs in the dermis was obtained by the in vivo skin microdialysis probe technology combined with Fourier series fitting. Excitedly, PAE-NEs skin pharmacokinetics profiles simulated by the FEA method with the input parameters showed great consistency with the actual in vivo performance. It can be concluded that the method of FEA based on a multilayer geometry model is a promising strategy to predict the TDDS skin pharmacokinetics.