Predicting Transdermal Fentanyl Delivery Using Mechanistic Simulations for Tailored Therapy

Transdermal drug delivery is a key technology for administering drugs. However, most devices are “one-size-fits-all”, even though drug diffusion through the skin varies significantly from person-to-person. For next-generation devices, personalization for optimal drug release would benefit from an augmented insight into the drug release and percutaneous uptake kinetics. Our objective was to quantify the changes in transdermal fentanyl uptake with regards to the patient’s age and the anatomical location where the patch was placed. We also explored to which extent the drug flux from the patch could be altered by miniaturizing the contact surface area of the patch reservoir with the skin. To this end, we used validated mechanistic modeling of fentanyl diffusion, storage, and partitioning in the epidermis to quantify drug release from the patch and the uptake within the skin. A superior spatiotemporal resolution compared to experimental methods enabled in-silico identification of peak concentrations and fluxes, and the amount of stored drug and bioavailability. The patients’ drug uptake showed a 36% difference between different anatomical locations after 72 h, but there was a strong interpatient variability. With aging, the drug uptake from the transdermal patch became slower and less potent. A 70-year-old patient received 26% less drug over the 72-h application period, compared to an 18-year-old patient. Additionally, a novel concept of using micron-sized drug reservoirs was explored in silico. These reservoirs induced a much higher local flux (µg cm^-2 h^-1) than conventional patches. Up to a 200-fold increase in the drug flux was obtained from these small reservoirs. This effect was mainly caused by transverse diffusion in the stratum corneum, which is not relevant for much larger conventional patches. These micron-sized drug reservoirs open new ways to individualize reservoir design and thus transdermal therapy. Such computer-aided engineering tools also have great potential for in-silico design and precise control of drug delivery systems. Here, the validated mechanistic models can serve as a key building block for developing digital twins for transdermal drug delivery systems.

Introduction

Transdermal drug delivery systems (TDDS) are used to deliver non-invasively moderately lipophilic, low-molecular-weight drugs to patients. Delivery through the skin, the largest human organ, enables controlled drug administration to achieve steady blood plasma concentrations without a distinct concentration peak (Perrie and Rades, 2014). Via the skin, drugs that have low bioavailability via oral administration, due to a high first-pass effect, can be effectively delivered. The ideal molecule for transdermal drug delivery has a molecular weight below 500 Da and a log partition coefficient (octanol-water) of 1 to 3 (Wiedersberg and Guy, 2014; Garg and Singh, 2018), so balanced lipophilicity. These molecules can cross the lipophilic stratum corneum but also diffuse through the hydrophilic viable epidermis and dermis, into the aqueous systemic circulation (Wiedersberg and Guy, 2014). Commercial TDDS have been developed for fentanyl, nitroglycerin, estradiol, nicotine, and testosterone, amongst others. In this multi-billion US dollar market (GVR, 2016; Markets, 2018), transdermal fentanyl patches currently provide the highest product sales, as a popular solution for around-the-clock opioid analgesia (Wiedersberg and Guy, 2014).

A key hurdle in these TDDS is that the pathway for delivery—the human skin—is very patient-specific. Absorption kinetics depend on the individual patient’s skin composition and the hygro-thermo-mechanical properties of each of its sublayers. Additionally, the patient’s metabolism, lifestyle, and bio-environment play a role. Studies have quantified the variability in drug uptake with body location (Roy and Flynn, 1990; Larsen et al., 2003; Sandby-Møller et al., 2003), age (Roskos et al., 1989; Sandby-Møller et al., 2003; Boireau-Adamezyk et al., 2014), gender (Sandby-Møller et al., 2003), ethnicity (Leopold and Maibach, 1996), skin hydration and disease state (Singh and Morris, 2011) as well as among patients within the same subject category (Larsen et al., 2003). The interpatient variability in the plasma concentration is caused by the complex combination of the TDDS release, transdermal drug absorption, circulation within the body and drug elimination. Any variation in these steps leads to inter- or intra-individual variability in the effect of the drug. The identified variability is dependent on the specific drug molecule and the composition and size of the target group (Farahmand and Maibach, 2009). The interpatient variability was thereby found to be larger or smaller than with oral delivery, depending on the drug (Farahmand and Maibach, 2009). This patient-induced variability makes precise transdermal dosage challenging. Clinical repercussions are that the therapeutic drug level in the blood is not always reached with certain patients or overdosing occurs when blood concentrations are out of the therapeutic range. Typical side effects are ineffective pain relief, skin irritation, respiratory depression, apnea, or, in some extreme cases, death (Schmid-Grendelmeier et al., 2006; Stanley, 2014).

Instead of “one-size-fits-all”, future TDDS and corresponding control strategies are designed to provide drugs for each patient at the correct rate for a specific body location (Breitbart et al., 2000; Amjadi et al., 2018). Defining a patient-specific therapeutic window was already proposed for certain drugs, including fentanyl (Woodhouse and Mather, 2000), since also the minimum effective concentration differed significantly between patients. However, personalizing transdermal drug delivery requires a clear insight into release and uptake kinetics and the associated biophysical processes that drive the uptake. Current transdermal drug delivery experiments have a rather low spatiotemporal resolution, large inter-sample variability, and are often performed for infinite drug reservoirs (Selzer et al., 2013). A typical example is monitoring the cumulative drug uptake in Franz diffusion cells by analyzing samples with high-performance liquid chromatography (HPLC) at time intervals of several hours. Another example is obtaining steady-state concentration-depth profiles over the skin via tape stripping (Larsen et al., 2003; Rim et al., 2005; Hansen et al., 2008). In several experiments, the skin samples are exposed to water reservoirs for a prolonged time. Water acts as a penetration enhancer for the skin and is often used in clinical dressings to increase drug diffusion (Bond and Barry, 1988). Therefore, the diffusion coefficients determined via Franz diffusion cells will represent an upper limit.

Mathematical modeling is an alternative method to gain complementary insights into the transport processes and to design and optimize next-generation TDDS in silico. The following methods have been used in the literature:

1. Mechanistic models that solve partial differential equations, for example at the macroscale, mesoscale (Naegel et al., 2013) and even cellular level (microscale; Figures 1A–C) with finite elements (Wittum et al., 2017).

2. Molecular dynamics at the subcellular level (e.g., lipid layers) to obtain transport properties and thermodynamic values of the system (e.g., partition coefficients) (Lundborg et al., 2018) that could then be used in a multiscale approach (Rim et al., 2009; Gajula et al., 2017) (Figure 1D). Additionally, alternative modeling strategies have been proposed that account for interactions at the molecular level (Schwöbel and Klamt, 2019) to calculate the skin permeability.

3. Compartmental models, which solve ordinary differential equations and consider each skin layer as a well-mixed compartment (Mitragotri et al., 2011; Selzer et al., 2013; Amarah et al., 2018). Usually, there is no discretization over the skin layer, although sometimes the skin layer is subdivided into a few compartments.

FIGURE 1

Figure 1 Overview of existing mechanistic models [adapted from (A) (Rim et al., 2009), (B) (Selzer et al., 2013), (C) (Naegel et al., 2009)].

Mechanistic models are of particular interest, compared to analytical solutions of drug uptake in the skin. Analytical models of the diffusion process are only valid for simple boundary conditions, which typically do not change in time or space and are derived for simple geometries (Siepmann and Siepmann, 2008). This impedes calculating, for example, a transdermal therapy that changes over time. Also, analytical solutions typically only calculate diffusion, without the presence of any physical binding or metabolization of the drug molecule in the skin. Mechanistic models provide drug concentrations and flow inside the skin and drug reservoir at each point in space and time in three dimensions, to the specific degree of detail that is dependent on the used method. Finite element, finite volume, or finite difference methods are typically used. These methods enable researchers to quantify, among other things, the time lag in drug release and subsequent uptake in the blood, local peak concentrations or fluxes, high-resolution concentration-depth profiles (Naegel et al., 2008), the depletion state of the drug reservoir, and the remaining drug amount stored in the skin. Mechanistic models also enable one to explore a large parametric space of process variables. Considerable work has been performed on mechanistic modeling for TDD (Table 1). These studies, however, do not explicitly focus on intra- or interpatient variability, a key step towards tailoring TDDS for individual patients or categories of patients, e.g., different age groups.

TABLE 1

Table 1 Non-exhaustive overview of mechanistic modeling work of transdermal drug delivery over the past decade (ordered chronologically).

In this study, our objective was to quantify the changes in transdermal fentanyl uptake with the patient’s age and the anatomical location where the patch was placed. We also explored how the drug flux from the patch could be altered as a function of the contact surface area of the patch reservoir with the skin. We applied mechanistic modeling to quantify in-silico differences between the scenarios. A finite-element model for transdermal drug delivery was developed in line with the current state-of-the-art. This mechanistic model was validated, and the sensitivity to the model parameters evaluated. With this model, we took three steps beyond the current state-of-the-art. First, we answered how much more fentanyl is taken up transdermally by a specific patient who applies the same patch now and 50 years in the future (due to skin aging). Second, we quantified how much more effective, or harmful, a transdermal patch is when placed at a different anatomical location. Third, we evaluated how the drug flux can be enhanced by miniaturizing the drug reservoir size, a process that takes advantage of transverse diffusion in the stratum corneum layer.

Materials and Methods

Continuum Model for Transdermal Fentanyl Delivery

Computational System Configuration

A mechanistic continuum model was built to simulate fentanyl release from a transdermal patch (reservoir) and subsequent uptake through healthy human skin. This study targets transdermal delivery of fentanyl, one of the most common drugs delivered transdermally (Wiedersberg and Guy, 2014). This synthetic opioid is approximately 100 times more potent than morphine (Moon and Chun, 2011). Fentanyl transdermal patches are used for patients with severe chronic pain, like cancer patients during treatment or at their end-stage (Muijsers and Wagstaff, 2001). Fentanyl has the appropriate lipophilicity and is sufficiently small to readily penetrate through the skin barrier and then reach the blood circulation. Therefore, conventional, first-generation transdermal patches can be used, where the uptake process through the skin barrier is diffusion-driven (Marier et al., 2006). The model and simulation were built and executed according to best practice guidelines in modeling for medical device design (Casey and Wintergerste, 2000; Morrison and F. D. A Guidance, 2016).

The geometrical setup, along with the boundary conditions, is depicted in Figure 2. The system configuration includes a square-shaped drug reservoir, which contains a finite amount of fentanyl, and the outer part of the human skin, namely the epidermis. TDDS are designed to deliver drugs at a nearly constant rate, similar to other controlled-release dosage systems (Chien and Lin, 2007; Perrie and Rades, 2014). The fentanyl patches, therefore, are commercially labeled with the targeted drug release rate, typically 12–100 µg h^-1 [Table 2, (Kress et al., 2010)] over the 72-h application period. Higher delivery rates are obtained by simply increasing the contact surface area of the patch (range between 4.2–42 cm², Table 2). Conventional transdermal fentanyl therapy consists of estimating empirically the initial dose (so patch size) for the patient, applying the patch transdermally, and replacing it every 72 h (Muijsers and Wagstaff, 2001). Patches are not allowed to be resized manually ad posterior to change the dose, as it may lead to a fast initial release of fentanyl (US Food and Drug Administration, 2005).

FIGURE 2

Figure 2 3D (A) and 1D (B) geometrical models of square-shaped drug reservoir and skin (not to scale; for 3D model, only one-fourth of the system was modeled due to the symmetry).

TABLE 2

Table 2 Delivery rate, active surface contact area, and drug content of selected, commercially available transdermal patches.

For the drug reservoir, a patch length L_pt of 40 mm was chosen. This value implies an active area of 16 cm², lying within the range reported for commercial transdermal patches for fentanyl [4.2-42 cm², Table 2, (Wiedersberg and Guy, 2014)]. The thickness of the patch (d_pt) was chosen to be 50 μm. Thus the volume of the patch reservoir was 80 mm³. These dimensions lead to a realistic initial drug content in the reservoir (mg), based on the initial concentration, as shown in section Boundary and Initial Conditions. The skin’s epidermis (ep) was composed of two layers: the stratum corneum (sc) and the viable epidermis (vep). The reason for this model configuration is that the lipophilic stratum corneum—a brick-mortar structure of lipid bilayers and corneocytes—has markedly different transport properties compared to the hydrophilic viable epidermis (Andrews et al., 2013). The stratum corneum exhibits the primary mechanical barrier function for drug delivery through the skin, together with tight junctions (Kirschner et al., 2010; Andrews et al., 2013; Matsui and Amagai, 2015; Bäsler et al., 2016), which are located inside the viable epidermis. The impact of tight junctions on the drug diffusion was not modeled separately but rather lumped into the transport properties of the epidermis. Tight junctions are rarely included explicitly in existing mechanistic models (Bäsler et al., 2016).

The dermis could also play a role in drug transport and drug retention during transdermal delivery (Menczel and Maibach, 1970). Previous works modeled the dermis 1) as an additional depot in which drugs can be stored (Heikkinen et al., 2015); 2) as an additional layer in which only diffusion and no drug removal by the blood flow was modeled (Manitz et al., 1998); 3) by adding a sink condition to the system (Grassi et al., 2011); 4) by modeling simplified capillary loops in the dermis (Calcutt and Anissimov, 2019). An even more realistic model of dermis would require modeling the patient’s blood flow since capillaries and vessels are present in the papillary and reticular dermis, respectively. Modeling blood flow would also require including the associated biochemical processes in the capillaries, an undertaking that was beyond this study’s scope. Therefore, the dermis was not included in the computational domain in the present study.

Two computational domains were constructed: 1) a one-dimensional (1D) model, where the drug reservoir was as wide as the skin, so this simulated unidirectional drug transport (Figure 2B). This case is representative of Franz diffusion cell experiments (Larsen et al., 2003; Rim et al., 2005; Bartosova and Bajgar, 2012) and has a low computational cost; 2) a three-dimensional (3D) model where the skin was wider than the drug reservoir to capture possible transverse diffusion of the drug (x and y directions in Figure 2A), since the drug spreads laterally, the actual delivery surface area becomes larger than the patch surface area (Vieille-Petit et al., 2015). The skin was extended by 1.7 mm [= 20 x (d_sc + d_vep)] on each side of the patch to avoid an impact of the transverse boundary on the simulated drug uptake kinetics at the interface between epidermis and dermis (surface 1 in Figure 2), as determined by a sensitivity analysis. Only one-fourth of the geometry was explicitly modeled by considering the symmetry in the system.

Governing Equations

Only the diffusion of fentanyl was solved, and isothermal conditions were assumed, namely, ambient body temperature. Neither water transport due to skin de-/rehydration, nor the resulting skin shrinkage or swelling was included. To derive the mass conservation equation, we started from the following equation, defined for each material i [kg m^-3] to describe the drug concentration $c_i^{\alpha }$ of substance α:

$\begin{array}{ll}\frac{\partial c_i^{\alpha }}{\partial t}+\nabla ·\left(-D_i^{\alpha }\nabla c_i^{\alpha }\right)=S_S^{\alpha }& (1)\end{array}$

where $D_i^{\alpha }$ is the corresponding diffusion coefficient or diffusivity (m² s^-1), $S_S^{\alpha }$ is a volumetric source term for substance α (kg m^-3s^-1), and t is the time (s). No source term was included in this study since the contributions of the following processes could be neglected (Naegel et al., 2013):

-Clearance of the drug via the blood capillaries and vessels into the body. As only the epidermis was included in the model (section Computational System Configuration), in which no blood capillaries or larger blood vessels are present, this effect was not included.

-Metabolization of the drug molecule within the skin by chemical reactions that may lead to a conversion of the drug into other compounds. Fentanyl is mainly metabolized in the liver and also by the intestinal mucosa for oral delivery, but not within the skin (Kim et al., 2016; Wishart et al., 2018). Therefore, metabolization within the skin is low and was not modeled.

-Adsorption of the drug molecule into the skin, and thus physical binding of the drug molecules. This process is different from the storage of unbound drug molecules in the tissue during transient drug uptake. The bound molecules could sometimes also unbind later and diffuse into the blood circulation. Adsorption of the drug in the skin is one of the processes that reduce the bioavailability of the drug. Bioavailability is the amount of drug administered through the skin that reaches the systemic circulation in an unchanged state. Bioavailability is a key pharmacokinetic characteristic that is determined by physical adsorption and chemical metabolization processes of the drug in the skin. The bioavailability of fentanyl transdermal delivery is very large [e.g., 92% in (McEvoy et al., 2017)]. This means that most of the drug which diffuses from the patch will reach the blood circulation. By considering this fact, we can conclude that one of the processes reducing this bioavailability—adsorption of fentanyl molecules in the skin—is also rather limited. Therefore, as an approximation, no physical adsorption source term was modeled in the skin domain.

A key phenomenon during drug uptake in multi-layer assemblies, such as the skin, is drug partitioning (de Monte et al., 2015). Partitioning implies that when a drug α is brought into contact with two materials (A and B), the drug concentration in these two materials will equilibrate to different values, namely $c_A^{\alpha }$ and $c_B^{\alpha }$. The ratio of these equilibrium concentrations is called the partition coefficient $K_{A/B}^{\alpha }$.

$\begin{array}{ll}{K}_{A/B}^{\alpha }=\frac{c_A^{\alpha }}{c_B^{\alpha }}& (2)\end{array}$

For drug partitioning in liquids, the octanol-water partition coefficient is often determined $(K_{o/w^{\alpha }})$, where values larger than 1 indicate drug lipophilicity and values smaller than one indicate drug hydrophilicity. The log $(K_{o/w^{\alpha }})$ is also often reported, where positive/negative values indicate lipophilicity/hydrophilicity, respectively.

Partitioning leads to a discontinuity in the drug concentration at the interface between the two materials, for example, between the stratum corneum and the viable epidermis. This discontinuity can affect the numerical stability of the simulation. An elegant solution is to substitute the dependent variable drug concertation $(c_i^{\alpha })$ for another variable $({\psi }_i^{\alpha })$ in Eq.(1) in the following way (Rim et al., 2008; Naegel et al., 2013):

$\begin{array}{ll}{c}_i^{\alpha }=K_i^{\alpha }{\psi }_i^{\alpha }& (3)\end{array}$

where ${\psi }_i^{\alpha }$ is termed in this study, the drug potential for every material i and is defined in [kg m^-3]. $K_i^{\alpha }$ is termed the drug capacity of the drug in the material i (-). Similar substitutions were made in other research fields (Datta, 2007; Defraeye et al., 2012; Defraeye and Verboven, 2017). This choice of dependent variable avoids numerical stability issues. This substitution is elaborated in Supplementary Material 1 and results in a single mass conservation equation instead of one for each material [Eq.(1)]:

$\begin{array}{ll}{K}^{\alpha }\frac{\partial {\psi }^{\alpha }}{\partial t}+\nabla ·\left(-D^{\alpha }{K}^{\alpha }\nabla {\psi }^{\alpha }\right)=0& (4)\end{array}$

where $K^{\alpha }=K_i^{\alpha }$ and $D^{\alpha }=D_i^{\alpha }$ for each material i, so these parameters were defined separately for each material. The dependent variable ψ^α is continuous throughout all materials and over all the interfaces (Naegel et al., 2013). The partition coefficient can also be defined as the ratio of drug capacities by combining Eq.(2) and Eq.(3):

$\begin{array}{ll}{K}_{A/B}^{\alpha }=\frac{K_A^{\alpha }}{K_B^{\alpha }}& (5)\end{array}$

Note that at all interfaces between different materials, continuity of fluxes is inherently maintained:

$\begin{array}{ll}\begin{array}{l}\left(-D_A^{\alpha }\nabla c_A^{\alpha }\right)·{\mathbf{n}}_{AB}=-\left(-D_B^{\alpha }\nabla c_B^{\alpha }\right)·{\mathbf{n}}_{BA}\hfill \\ \left(-D_A^{\alpha }{K}_A^{\alpha }\nabla {\psi }^{\alpha }\right)·{\mathbf{n}}_{AB}=-\left(-D_B^{\alpha }{K}_B^{\alpha }\nabla {\psi }^{\alpha }\right)·{\mathbf{n}}_{BA}\hfill \end{array}& (6)\end{array}$

where n is the unit vector normal to the interface.

Material Properties and Transport Characteristics of Skin and Patch

Fentanyl is a synthetic opioid used as a pain medication. It has a low molecular weight (337 Da) and is moderately lipophilic with a log(K_o/w) of 3 to 4 (Rim et al., 2009; Wiedersberg and Guy, 2014; Kim et al., 2016). The material transport properties of the skin components and the drug reservoir are given in Table 3 for fentanyl for the different cases that were simulated (see section Spatial and Temporal Discretization). Values were taken from the literature (Rim et al., 2005; Rim et al., 2009). The different material property data sets that were used (validation study versus parametric study) led to similar uptake kinetics (see Supplementary Material 3). The drug capacities are derived from the partition coefficients via Eq.(5), as the latter are typically available from measurements. To be able to do so, one drug capacity needed to be fixed, and therefore K_pt^α was set (arbitrarily) equal to one. Similar to most other simulation studies (Rim et al., 2009; Naegel et al., 2011; Gajula et al., 2017), the diffusion and partition coefficients were taken as constants and isotropic, meaning that they are independent of the drug concentration because more detailed data was not available.

TABLE 3

Table 3 Material transport properties used in the model for different configurations for fentanyl.

Boundary and Initial Conditions

The drug was assumed to be removed from the computational domain only via the interface with the dermis. In the dermis, it will be transported by diffusion and blood flow. To this end, a constant concentration (and potential), equal to zero, was imposed at the bottom of the epidermis (Figure 2), as in previous studies (Rim et al., 2005; Rim et al., 2008). This condition represents a Dirichlet boundary condition. This assumption is justified by the fact that the drug concentration in the dermis is very low because of a higher drug diffusion coefficient, drug partitioning between dermis and epidermis, and the fact that the blood flow extracts drugs. Zero-flux conditions were imposed at all vertical boundaries. At t=0, the skin was assumed to be drug-free. The initial drug concentration in the patch was set at 80 kg m^-3, according to a previous study (Rim et al., 2005). This dose implies a total initial amount of fentanyl in the reservoir m_pt,ini = 6.4 mg, which corresponds to amounts typically present in commercially available patches (Table 2). Complete contact between the patch and the skin was assumed, without any inclusion of air layers or discontinuities like hairs or skin roughness.

Spatial and Temporal Discretization

The finite element mesh was built based on a grid sensitivity analysis on the 1D and 3D models. Using Richardson extrapolation, the spatial discretization error on the total mass flux to the dermis was estimated to be 0.1% for both 1D and 3D models (Roache, 1994; Franke et al., 2007). The grid consisted of 120 quadrilateral finite elements (1D, elements with a size of about 1 µm) and 107,000 hexahedral finite elements (3D) for the base case. The grid was gradually refined towards the different material interfaces to enhance numerical accuracy and stability, because the largest gradients are found at such interfaces, particularly at the initial stage of the uptake process.

Starting from these initial conditions, the transient simulations calculated a drug uptake process that lasted 240 h (10 days), to capture the drug uptake history until depletion. These simulations applied adaptive time-stepping, with a maximal time step of 600 s (10 min). This time step was chosen to ensure high temporal resolution for the output data and was determined from a sensitivity analysis.

Alternative Configurations

The base case (Table 3) mimicked the drug uptake in the skin, as released from a finite drug reservoir. Additional configurations were simulated with partially different geometries, process conditions, or reservoir/skin material properties:

-Validation of the drug release and uptake with experimental data from a previous study (Rim et al., 2005). A detailed description of the experiments and the corresponding simulations is given in Supplementary Material 2 and Table 3.

-Sensitivity analysis to multiple model parameters (detailed in section Sensitivity to Model Parameters)

-An infinite reservoir with a very high diffusion coefficient that cannot be depleted (detailed in section Infinite Reservoir and Patch Removal After 72 h).

-The case where the patch was removed after 72 h, which was performed by changing the boundary conditions (detailed in section Infinite Reservoir and Patch Removal After 72 h).

-Variability with respect to the anatomical location where the patch was placed on the patient’s body and the age of the patient (detailed in section Patient-Specific Parameters).

-Reservoir contact surface area and size, in terms of its transverse dimensions (x,y in Figure 2; detailed in section Reservoir Contact Surface Area).

Unless specified otherwise, 1D unidirectional models were used (Figure 2). To support this decision, the differences with a 3D model and the impact of the anisotropic transport properties between transverse and longitudinal directions were quantified for the evaluated patch widths in Supplementary Material 3 and section Validation.

Sensitivity to Model Parameters

The sensitivity of the drug uptake process to various model parameters was explored to identify the ones with the largest impact. To this end, the relative sensitivity S_{U, Xj} of a process quantity U(t) (e.g., drug flux) to a change in a model input parameter X_j over time was calculated using partial derivatives:

$\begin{array}{ll}{S}_{U,X_j}=\left|\frac{\frac{\partial U(t)}{U(t)}}{\frac{\partial X_j}{X_j}}\right|=\left|\frac{\frac{U_{X_j+\Delta X_j}(t)-U_{X_j-\Delta X_j}(t)}{U(t)}}{\frac{(X_j+\Delta X_j)-(X_j-\Delta X_j)}{X_j}}\right|=\left|\frac{U_{X_j+\Delta X_j}(t)-U_{X_j-\Delta X_j}(t)}{2\Delta X_j}\frac{X_j}{U(t)}\right|& (7)\end{array}$

where ΔX_j/X_j was set equal to a 1% deviation from the nominal value of X_j. The sensitivity to the following parameters (X_j) was probed: d_sc, d_ep, $D_{sc}^{\alpha }$, $D_{ep}^{\alpha }$, $K_{pt/sc}^{\alpha }$, $c_{pt,ini}^{\alpha }$. The following quantities, U(t), were evaluated at specific time points (12, 24, 48, and 72 h): (1) the uptake flux across the skin into the dermis, so into the blood g_bl,up [kg m^-2 s^-1]; (2) the total amount of drugs taken up via the dermis by the blood flow m_bl,up(t) [kg]. Such a sensitivity analysis is an essential step to designing drug delivery systems and therapies that achieve a constant drug delivery for the patient.

Infinite Reservoir and Patch Removal After 72 h

An infinite reservoir with a very high diffusion coefficient was simulated. As this reservoir cannot be depleted, the gradient over the skin remains constant once a steady state is reached. To this end, the boundary conditions were adjusted. For these simulations, the patch was removed in the system configuration and a constant concentration $(=c_{pt,ini}^{\alpha })$ was imposed at surface 2 (Figure 2). To simulate the removal of the patch after 72 h, a zero-flux condition was imposed at surface 2, and only transport in the epidermis was solved.

Patient-Specific Parameters

The patient-specific uptake of fentanyl through the skin epidermis was quantified in two steps. First, the intra-patient variability was analyzed, namely by investigating the effect of the anatomical application site on the human body of the transdermal drug delivery device. The different body sites that were tested and their corresponding stratum corneum and viable epidermis thickness are presented in Table 4. Simulations were performed with the average thickness (µ_i), but also with thicknesses that corresponded to µ_i -/+ 2σ_i, respectively, where σ_i is the standard deviation of material i. This range corresponds to the 2.5% and 97.5% quantiles of the sublayer thickness, and, thus, about 95% of the thickness range was covered.

TABLE 4

Table 4 The thickness of stratum corneum and viable epidermis at different body sites—mean (standard deviation).

Second, the interpatient variability on the drug uptake kinetics was evaluated, namely by assessing patients from different age categories. The stratum corneum thickness increases with age (Boireau-Adamezyk et al., 2014). In addition to the thickness of the stratum corneum, which is the main barrier for transdermal drug delivery, also the composition of the skin changes, where the lipid content changes. These compositional changes, in turn, affect the diffusion coefficient and the partition coefficient. Previous studies suggest that with increasing age, the diffusion coefficient, and thus the permeability of the skin slightly decreases (Thakur et al., 2009; Boireau-Adamezyk et al., 2014). This decrease in permeability of the skin does not significantly affect the permeation for lipophilic drugs (Roskos et al., 1989). Therefore, in this study the possible changes in the composition of the stratum corneum layer due to age and the resulting diffusion and partition coefficient were not investigated. Note that these effects would even strengthen the reduction of transdermal drug uptake with age. Only the effect of age on the stratum corneum thickness was included, so only the skin layer thicknesses were altered, while the other material properties were assumed to remain constant. Identifying the intra- and interpatient variability of the drug uptake kinetics will help to design drug delivery systems and therapy so that every patient can receive an optimal drug dose.

The following correlation was used to relate the stratum corneum thickness [d_sc (m)] and age [A (a), i.e., years] for the dorsal forearm, as reported previously (Boireau-Adamezyk et al., 2014) (40 female subjects, 10 in each of the following age categories: 18–30, 30–40, 40–55, and 55–70 years):

$\begin{array}{ll}{d}_{SC}=(0.125A+11.8)·{10}^{-6}& (8)\end{array}$

Reservoir Contact Surface Area

In addition to the 40 mm wide drug reservoir, other reservoir sizes were evaluated. The rationale for this decision was to quantify how the contact surface area (patch length or width L_pt) of the patch affects the drug flux to the skin. The following patch widths were evaluated: 40 mm (base case), 4 mm, 400, 40, 4 µm. Since smaller reservoirs are depleted faster, simulations with an infinite reservoir (section Infinite Reservoir and Patch Removal After 72 h) were also conducted.

Numerical Implementation and Simulation

The model was implemented in COMSOL Multiphysics® software (version 5.4, COMSOL AB, Stockholm, Sweden), a finite-element-based commercial software program. This software was verified by the code developers. Therefore, additional code verification was not performed by the authors. Transient diffusive drug transport [Eq.(4)] in the patch and skin during drug release and uptake was solved using the partial differential equations interface (coefficient form). The conservation equation was solved for the dependent variable ψ. Quadratic Lagrange elements were used with a fully-coupled direct solver, which relied on the MUltifrontal Massively Parallel sparse direct Solver (MUMPS) solver scheme. The tolerances for solver settings and convergence were determined by means of sensitivity analysis in such a way that a further increase in the tolerance did not alter the resulting solution.

Metrics

The simulated drug delivery process was analyzed quantitatively by calculating several metrics. With respect to the stored drug amount, we quantified the remaining (residual) amount contained in the patch as a function of time {m_pt,res(t) [kg]}, and the total amount of drugs stored (present) in the skin (viable epidermis and stratum corneum) as a function of time {m_ep,stor(t) [kg]}.

For the transported drug quantity, the drug amount released by the patch {m_pt,rel(t) [kg]} was quantified as a function of time, which is the cumulative integrated flux over surface 2 (Figure 2) over time. The corresponding release flow rate {G_pt,rel [kg s^-1]}, namely the slope of the m_pt,rel(t) curve, was also determined. Furthermore, the drug amount that was taken up by the dermis was also determined. Due to the large bioavailability of fentanyl, we assumed that almost all drugs exit the computational system. Hence, the amount uptaken by the dermis equals the drugs taken up into the blood flow. Using the outgoing drug flow, the drug amount that was taken up by the blood flow, m_bl,up(t) [kg], was determined as a function of time. It is the cumulative integrated flow over surface 1 (Figure 2) over time. The corresponding uptake flow rate {G_bl,up [kg s^-1]} was also determined, which is the slope of the curve. Using this parameter, it was possible to calculate the flux across the skin into the blood {g_bl,up [kg m^-2 s^-1]}.

The maximal fentanyl concentration in the stratum corneum {c_sc,max(t) [kg m^-3]} was also quantified. This quantity is important because concentrations of specific drugs that are too high may induce irritation (Schmid-Grendelmeier et al., 2006). With these quantities of stored and transported amounts of drugs, the mass balance of the drug can be written as:

$\begin{array}{ll}{m}_{pt,ini}=m_{pt,res}(t)+m_{ep,stor}(t)+m_{bl,up}(t)& (9)\end{array}$

Out of the drug uptake profiles, m_bl,up(t), the uptake kinetics were assessed by defining the fractional drug release of the patch (Y_bl,up):

$\begin{array}{ll}{Y}_{bl,up}=\frac{m_{bl,up}(t)-m_{pt,fin}}{m_{pt,ini}-m_{pt,fin}}& (10)\end{array}$

where subscripts ini and fin represent the initial and final drug amount in the patch when the patch is depleted, where m_pt,fin =0. m_bl,up(t) represents the amount of the drug taken up by blood flow. From the definition of Y, the half-uptake-time (HUT, t_1/2) can be calculated as the time required to take up half of the initial drug amount in the patch via the skin. HUT is a useful parameter to characterize and compare the release behavior of the patch, as a single value can be used to characterize the uptake kinetics.

Results and Discussion

Validation

The comparison with previous experimental data (Rim et al., 2005) in Figure 3 identifies how accurate our mechanistic model predictions on drug uptake through the epidermis are for two different initial drug concentrations in the reservoir. For both initial concentrations, the simulations are in good agreement with experimental data, i.e., within the error bars, during the first 27 h with a root mean square deviation (RMSD) of 0.13 and 0.14 µg cm^-2 h^-1 at 60 kg m^-3 and 80 kg m^-3, respectively. This initial period of the drug uptake, characterized by a strong increase in the flux, is thus predicted accurately. This increase is caused when drugs are being transported into the epidermis. As a result, the drug concentration in the epidermis increases, since a part of the drugs is stored there. For the remaining 2 days, there are larger deviations. The experimental decrease in the flux exceeds that of the simulations. Here, a quasi-steady-state condition sets in since drugs do not accumulate anymore in the epidermis. The skin’s capacity to hold drugs is reached, by which the stored amount of drugs remains rather constant (see section Drug Release and Uptake in the Skin). As such, the drugs entering the epidermis primarily diffuse through. However, full steady-state conditions with a constant flux are never reached, since the drug concentration in the patch - a finite reservoir—decreases over time. As a result, the uptake flux slowly decreases over time. This predicted decline has a rather constant slope in the simulations. The experiments predict a much steeper decline in flux than in the simulations. The simulations seem thus to miss capturing one or more key processes.

FIGURE 3

Figure 3 Drug fluxes of fentanyl at surface 1 (g_bl,up), so leaving the epidermis, from experiments (Exp.) and simulations (Sim.) for two initial patch concentrations (60 kg m^-3, 80 kg m^-3) as a function of time. The fentanyl molecule is illustrated as well, where the oxygen is depicted in red, carbon in black, hydrogen in white, and nitrogen in blue.

Specific model improvements that could enhance the accuracy are discussed in section Mechanistic Modeling. However, for this particular discrepancy during the later stage of the drug uptake process, a particular factor is likely to be responsible in part. Due to the high diffusion coefficient of the patch, the concentration in the patch is very uniform and almost equals that at the patch-skin interface. As such, the main concentration gradients occur over the skin, as predicted in-silico. Over time, the concentration in the patch decreases, which leads to a reduced gradient over the skin. The steep decrease in the flux in the experiments could indicate that also a concentration gradient inside the patch could be present. This would be caused by a lower diffusion coefficient in the patch than currently used in the simulations. The path length for the molecules to travel through the skin would progressively increase over time, hence reducing the flux more than based on a concentration gradient over the skin alone.

The diffusion coefficient of the patch $D_{pt}^{\alpha }$ was determined experimentally by performing a release experiment of the patch directly in a receptor medium and fitting an analytical expression to the results. Perhaps this type of experiment in an aqueous medium leads to too high diffusion coefficients compared to when the patch placed on the skin. If the diffusion in the reservoir holding the drug should be restricted more, a better agreement with experiments over the entire uptake process could likely be obtained but this hypothesis should be further explored.

The current discrepancy of the model with the experimental data leads to a lower predicted amount of fentanyl than in reality. As a result, the targeted pain relief, as predicted from the in-silico drug dose, will not be reached in reality. However, the patient will certainly not get a higher drug dose than predicted by the model. The model thereby provides a conservative drug uptake estimate, which is safe for the patient. Note that the differences in these predicted fluxes and the cumulative drug amount that was taken up by blood flow between the 3D (cylindrical, Supplementary Material 2) and unidirectional (1D) transport models were < 0.4%. Therefore, a 1D transport model can be used as a viable alternative for the 3D model.

Comparison With Commercial TDDS

The drug uptake from the simulated transdermal patch (base case) is compared with that of commercial TDDS in Figure 4. From our simulations, we obtain a peak uptake rate by the blood of 23 µg h^-1 after 7–10 h (Figure 4). This uptake decreases to 18 µg h^-1 after 72 h for a 16 cm² patch. The delivery rate is rather constant, but there is a slight decline due to the reduction of reservoir concentration, which is the driving force for drug transport. Commercial patches, however, only report the targeted steady-state value that does not reflect this slight decline in flow rate, because these TDDS are designed to deliver the drug at a nearly constant rate. Our results lie between the performance of Durogesic® DTrans® 12 µg h^-1 and 25 µg h^-1 patches, for example, which have a surface area of 5.25–10.5 cm². Thereby, the patch used in this study, with a surface area of 16 cm², delivers a flux that is slightly lower than the Durogesic® DTrans® patches. The reason is a different patch design or material composition. In summary, our mechanistic model produces uptake rates and kinetics in a similar range as commercial products. Combined with the validation study results, this tool is reliable for product design and optimization.

FIGURE 4

Figure 4 Drug flow as a function of time for drugs released by the patch (g_pt,rel), for drugs taken up by the blood flow (g_bl,up) for the base case and commercial fentanyl patches. The results of an infinite reservoir are also shown.

Sensitivity Analysis to Model Input Parameters

The relative sensitivity of the flux across the skin into the blood (g_bl,up) to the input parameters is shown in Figure 5A. These results quantify the parameters that affect the most the predicted drug uptake, and how this sensitivity changes over time. The total amount of drugs taken-up after a certain period is shown in Figure 5B. A sensitivity value S_U,Xj of one implies that the impact on the solution (U_{Xj + ΔXj} - U_Xj) is more significant than the prescribed disturbance (1%, ΔX_j) of the model input parameters [Eq.(7)].

FIGURE 5

Figure 5 Relative sensitivity to different model parameter inputs for the mass flux taken up by the blood [U(t) = g_bl,up] as a function of time (A) and the total amount of drugs taken up after 12-24-48-72 h [U(t) = m_bl,up] (B, sc: stratum corneum, vep: viable epidermis, diff. coeff.: diffusion coefficient, part. coeff.: partition coefficient, conc.: concentration).

Concerning the magnitude of the sensitivities, the data show the largest sensitivity for the model parameters that describe the stratum corneum. This finding is not surprising given that this epidermal sublayer has the most significant resistance to drug transport, and exhibits the skin’s primary barrier function. The partition coefficient has the largest overall impact of all model parameters. The thickness and the diffusion coefficient have similar sensitivities due to their comparable role in the conservation equation [Eq.(1)].

There is a distinct temporal sensitivity to the different model parameters, especially the thickness and diffusion coefficient of the stratum corneum and the initial drug concentration in the reservoir. This temporal sensitivity originates from the transient nature of the uptake process. Initially, loading of the drug in the skin occurs (first 10 h), which involves drug storage and transport through the skin. After this initial period, the sensitivity is more constant over time. The highest temporal sensitivity occurred for the initial drug concentration, with a sensitivity to the mass flux that reduces to only a few percent after the first 10 h (Figure 5A). For designing drug delivery systems and therapy, the initial drug concentration in the patch is an important factor, and also comes into play when replacing the patch.

In summary, uncertainty in model parameters impacts the solution. Still, this effect is smaller than the disturbance itself (i.e., S_U,Xj <1), except during the very early stage of drug uptake (< 24 h). Nevertheless, the solution is particularly sensitive to the model parameters of the stratum corneum, especially the partition coefficient, so these parameters must be known as accurately as possible.

Drug Release and Uptake in the Skin

Release and Uptake Kinetics

The drug release and uptake kinetics from the patch into the skin are presented in Figure 6 for the base case, specifically by analyzing the contributions to diffusion and storage processes (including partitioning). These data display what the temporal delay is between drug release and uptake, and how much drugs are stored in the different skin layers over time. After drug release initiation, it takes roughly 30 min before the drug diffused through the 85 μm thick epidermis and can be started to be taken up by the blood via the capillaries in the dermis. In parallel to diffusive transport, drugs accumulate in the skin during the first 6–7 h. Despite its small thickness, the stratum corneum stores approximately 99% of the drug in the epidermis, because of its much larger capacity due to partitioning compared to the viable epidermis (K_i^α, Table 3). Equilibrium occurs after this initial period because the skin’s capacity to store drugs is reached then. After that, the drug primarily diffuses through the skin, and the stored amount remains rather constant. The stored drug amount in the skin is, however, only a fraction of the total drug amount initially present in the patch, with a maximal value of 0.16 mg, or 2.4% of the initial content. This low value indicates that the drug bioavailability from the model was > 97.6%. Due to the decreasing drug concentration in the patch (a finite reservoir), the uptake flux slowly decreases over time, simply because the concentration gradient decreases.

FIGURE 6

Figure 6 (A) Amount of drug stored in the epidermis [m_ep,stor(t)] and patch [m_pt,res(t)], and cumulative amounts of drugs released by the patch m_pt,rel(t) and taken up by the blood flow m_bl,up(t) as a function of time for the base case. The initial amount of drugs is also shown. (B) Corresponding flux released by the patch into the epidermis (g_pt,rel) and taken up by blood (g_bl,up) as a function of time, where one insert shows the percentage for the maximal value. The results of an infinite reservoir are also shown.

The patch is depleted by half of its initial amount after 190 h (almost 8 days). Typically, fentanyl patches are replaced every 72 h or 3 days (e.g., Table 2). After this 72-h application period, 77% of the initial drug amount is still present in the patch for the base case. However, transdermal patches are designed to deliver the drug at a controlled rate to achieve a constant blood plasma concentration rather than delivering the entire amount of drug and depleting the reservoir (Perrie and Rades, 2014). Therefore, the concentration in the patch cannot reduce too much, as this would imply the driving force for drug transport, i.e. the concentration difference, would become too low. As such, a significant drug concentration should still be present when removing the patch, to guarantee a rather constant drug flow (McEvoy et al., 2017). Also, in our validation study, only 31% of the initial drug amount diffused through the skin after 72 h. Larsen (Larsen et al., 2003) reported absorbed amounts below 20%.

For comparison, the results of an infinite reservoir with high diffusive drug transport are also shown in Figure 6. As expected, a constant flux is reached after an initial uptake period. In other words, there is a linear increase in the drug amount taken up with time. Ideally, TDDS are designed to deliver drugs at a constant rate. The resulting constant fentanyl flux (1.5 μg cm^-2 h^-1) is similar to the range reported for commercial patches [Table 2; (Wiedersberg and Guy, 2014)], namely 1.7–3.0 μg cm^-2 h^-1. The finite reservoir deviates from this constant rate. After 72 h, the flux is only 80% of the maximum (insert in Figure 6B). The infinite reservoir overpredicts the amount of drug uptake by 17% after 72 h compared to the finite reservoir (base case), and this difference is significant.

Additionally, Supplementary Material 3 shows that the differences with a 3D model are limited, and the anisotropic transport properties between transverse and longitudinal direction have a limited impact on the evaluated patch width.

Spatial and Temporal Resolution

We compare the high spatiotemporal resolution data from the simulations against typical data obtained in TDDS experiments. This endeavor identifies additional insights and benefits gained from mechanistic modeling. In experiments, drug uptake is typically measured using a Franz diffusion cell at discrete time intervals (e.g., every 5 h) and subsequent measurements of the samples via HPLC (Larsen et al., 2003; Rim et al., 2005). From the measured drug amount taken up (μg) during the sampling interval (h) and the patch’s surface area (cm²), the time-averaged flux through the entire skin sample can be calculated (μg cm^-2 h^-1). Even if the total amount of drug taken up is correct, such temporal averaging masks local flux peaks in time and the associated maxima in concentrations. This information is essential because concentrations of specific drugs that are too high may induce skin irritation (Hogan and Maibach, 1990; Brockow et al., 2013). The discrepancies induced by such temporal averaging are quantified in Figure 7, as derived from our simulation data. The average flux for different time-averaging intervals is calculated by grouping our simulation results over specific intervals because these can be derived from the fluxes at each point in time. Averaging intervals >5 h substantially underestimate fluxes by approximately 20%. However, this accuracy strongly depends on the process kinetics during that timeframe. Simulations also allow researchers to mimic experimental conditions in a deterministic way, without suffering from biological variability and the resulting uncertainty.

FIGURE 7

Figure 7 Flux taken up by the blood as a function of time for 3 time-averaging intervals, based on simulation data. The horizontal dashed lines indicate the time range over which the averaging took place, i.e., 2.5, 5, or 10 h.

To further illustrate the added spatiotemporal insights on drug transport obtained by simulations, the vertical concentration profiles through the patch and skin are shown in Figure 8. The experimental counterpart to obtain such profiles would be tape stripping (Lademann et al., 2009). This process, however, only extracts the profiles at a single point in time, where typically steady-state conditions are targeted. The obtained 1D spatial resolution by tape stripping is in the micron range, but is strongly dependent on the anatomical site, age, stratum corneum thickness, number of cell layers, and corneocytes, among others (Lademann et al., 2009). Contour plots of concentration and potential from the simulations are shown in Figure 9. Partitioning is indicated by the much higher concentrations of the moderately lipophilic drug in the stratum corneum compared to the viable epidermis. The discontinuities at the interface of the skin-patch layer are visible and contrast the continuous distribution of the dependent variable drug potential ψ^α [Eq. (4)] through all layers. These large concentration jumps challenge the numerical stability of the simulation, which is why the conservation equations were solved against the drug potential instead of drug concentration in this study.

FIGURE 8

Figure 8 Concentration profiles at different points in time as a function of the height of the model (z-direction, scaled with the epidermal thickness) in the vertical centerline (z-axis) of the patch for the base case: (A) entire model (logarithmic scale), (B) reservoir, (C) stratum corneum, (D) viable epidermis.

FIGURE 9

Figure 9 Color contours of (A) drug concentration (c_i^α) and (B) drug potential (ψ^α) in the drug reservoir and epidermis for different points in time. Different scaling is used in (A, B) to improve clarity. Note that the maximal range (80 kg m^-3) is not depicted.

After initially loading the epidermis with the drug, a quasi-steady-state diffusion process develops, with the typical linear concentration distribution over each epidermal layer (Figures 8C, D). The concentration reduction in the patch (Figure 8B), however, still induces a slight temporal shift in the concentration profiles in the epidermis (Figures 8C, D).

Removing Patch After 72 h

The simulations enable us to analyze what occurs within the epidermis when the fentanyl patch is removed after the recommended 72-h period. This analysis is possible because simulations can quantify volume-averaged drug contents as well as surface-averaged fluxes. In Figure 10, the drug uptake amount in the blood, the corresponding flux, and the drug storage in the skin are shown as a function of time. A very small drug amount was stored in the skin (0.12 mg after 72 h). Because only this amount can further diffuse into the blood once the patch is removed, the drug amount that is taken up steeply decreases after removal. However, it still took approximately 24 h before this residual amount diffused from the skin to the blood. This effect is more pronounced for drugs that can be stored in higher amounts in the skin, namely those with a larger drug capacity (K_i^α).

FIGURE 10

Figure 10 (A) Amount of drug stored in the epidermis [m_ep,stor(t)] and the cumulative amount of drugs taken up by the blood flow m_bl,up(t) as a function of time for the case where the patch is removed after 72 h. (B) The corresponding flux that is taken up by blood (g_bl,up) as a function of time.

Patient-Induced Variability for Anatomical Location

The impact of the anatomical location of the drug reservoir on the human body was also investigated. The absorption of drugs through the skin differs in different body locations due to a different thickness of the stratum corneum or the presence of the hair follicles (Feldmann and Maibach, 1967). In Figure 11, the amount of drug uptake is given for three-body sites that differ in terms of epidermal thickness. The HUT (section Metrics) is also indicated for each body location, as calculated based on the curves with average thickness. The HUT for the base case was 190 h.

FIGURE 11

Figure 11 Cumulative amount of drug taken up in the blood as a function of time for the forearm (A), shoulder (B), buttock (C) epidermal thickness averaged for all body sites (D). For comparison, the results of the base case are also depicted. The range corresponding to the simulated average thickness values +/- 2 x standard deviation (µ_i -/+ 2σ_i) is depicted in gray to quantify interpatient variability via the 2.5% and 97.5% quantiles. The thickness of the epidermis (d_ep) and half-uptake time (HUT) are indicated as well..

The largest drug uptake was observed on the shoulder, while the smallest was obtained for the forearm. The difference in drug uptake (mg) between these two body sites was 36% (relative to that of the forearm) after 72 h. The results for the average epidermal thickness over all the body sites agree with the base case. Interpatient variability in the stratum corneum and epidermis thickness directly manifested itself in drug uptake rates. For the forearm, the variation in total drug uptake (mg) over the simulated range (µ_i -/+ 2σ_i) is 113%, and this variation is not evenly distributed across the average value. Clinically, this finding implies that specific patients can have blood serum concentrations that can be too high (toxic) or too low to be effective for that specific patient. This interpatient variability for a specific anatomical location induces significant variability on the drug uptake results. This spread makes it challenging to distinguish significant differences between anatomical locations in clinical experiments. In contrast to the present study, specific previous studies reported limited differences in fentanyl uptake between the anatomical location (Roy and Flynn, 1990; Larsen et al., 2003).

Patient-Induced Variability for Age

The impact of patient age on drug uptake kinetics was investigated by calculating variations in the stratum corneum thickness with age via Eq. (8). In Figure 12, the drug uptake and flow rate are given for different age groups [as previously defined (Boireau-Adamezyk et al., 2014)]. There is a 26% difference in the total drug amount taken up after 72 h between an 18- and a 70-year-old patient (relative to the 18-year-old patient). This finding suggests that if a patient applies the same fentanyl patch now and 52 years later, this patient will receive 26% less drug when s/he is older; therefore, the patch will be less effective. When examining flow rates, one could even consider performing therapy using a patch with a higher dose and/or a larger surface area. The current mechanistic model can even be used to define an age-specific dosage systematically.

FIGURE 12

Figure 12 (A) Cumulative amount of drug taken up in the blood as a function of time for different age categories. The results of the average age of each group are shown, as well as the minimal and maximal ages. (B) Drug flow taken up by the blood as a function of time for different age groups and commercial fentanyl patches. For both A, B the gray band indicates the entire range of ages that is evaluated (18–70 years).

Furthermore, the time before the minimal effective concentration is reached in the blood differs with age. As an example, it took 23 h to uptake 0.5 mg for an 18-year-old patient, whereas it took 9 h longer when the patient was 70 (Figure 12A). Our simulations showed that, with aging, the patch delivered drugs more slowly and less potently. Mechanistic simulations enable researchers to quantify this difference deterministically and theoretically, without introducing additional statistical uncertainty concerning interpatient variability.

Impact of Contact Surface Area of the Reservoir

We explored how the reservoir width, and thus the contact surface area affected the released flux. For normal TDDS, the reservoir is much wider than the epidermal thickness, which is the longitudinal transport pathway for the drugs. This phenomenon leads to unidirectional transport. Hence, the 3D edge effect induced by transverse diffusion at the edges of the patch is negligible. This examination aimed to identify whether a reservoir with a smaller contact surface area released drugs faster and at a higher rate than a standard patch. Since transport will occur in 3D in the case of patches with a smaller contact area, this transverse diffusion could induce higher fluxes. In Figure 13, the released flux (surface 2) is shown as a function of time for all reservoir surface areas for finite and infinite reservoirs. In Figure 14, the released flux at equilibrium (for an infinite reservoir) is shown as a function of reservoir width (L_pt), where the reservoir width is made dimensionless with the epidermal thickness (d_ep = 85 μm).

FIGURE 13

Figure 13 Flux released by the patch into the skin (g_pt,rel, surface 2) as a function of time for different reservoir sizes (L_pt), so contact surface areas, for (A) finite reservoir, (B) infinite reservoir.

FIGURE 14

Figure 14 Flux released by the patch into epidermis after 72 h (g_pt,rel) for the infinite reservoir as a function of the reservoir size (L_pt), so contact surface area, scaled by the epidermal thickness (d_ep). In the schematics, dimensions of the patch and skin thickness are to scale, but the skin width is not to scale.

The flux leaving the reservoir is dependent on the size so contact surface area of the drug reservoir in contact with the skin, for both finite and infinite reservoirs. For the finite reservoirs, the depletion of the smaller reservoirs causes the flux to decrease sharply over time. This depletion can, however, be mitigated by increasing the thickness of the reservoir (d_pt) or by connecting all small reservoirs to a large bulk reservoir. The infinite reservoirs evolve to a steady-state, a condition that is more convenient for comparison of the reservoir contact surface area. Once the size (L_pt) enters the submillimeter range, or the patch size goes below approximately 10 x d_ep, the flux increases to more than double of that of a conventional patch (base case). This phenomenon occurs because the drug can diffuse in three dimensions instead of predominantly one direction for the large reservoirs.

Furthermore, the transverse diffusion coefficient of the stratum corneum is a few orders of magnitude larger than the longitudinal one [(Rim et al., 2009); Table 3]. Thereby, for smaller contact surface areas, longitudinal and transverse diffusion occur, a process that induces a higher flux at the contact interface. This edge effect is illustrated in the contour plots presented in Figure 15. For a reservoir size (L_pt) of 40 μm, the release rate increases by a factor of 20 at equilibrium for an infinite reservoir compared to the base case. For L_pt = 4 μm, the increase was 200-fold. Note that these factors decrease to 3 and 18, respectively, when the transverse diffusion coefficient would just equal the longitudinal one (results not shown). This data implies that transverse diffusion is a key parameter for the large observed differences, partially due to the higher transverse diffusion coefficient of the stratum corneum layer, especially when the reservoir is not very wide compared to the skin thickness. Note that for small patches, drug storage in the skin will delay uptake into the blood, because it takes longer to saturate the skin sublayers due to their lower total flow rate (µg h^-1). It would take longer to reach an uptake equilibrium (g_bl,up), namely approximately 35 h for L_pt = 4 μm versus 15 h for L_pt = 40,000 μm. The amount stored in the skin is, in all cases equal to, or smaller than the base case. Note that the smallest reservoirs are of the same size as the corneocytes (Figure 1) but still much larger than the lipid bilayer thickness [approximately 10¹ nm (Das and Olmsted, 2016)]. As such, the drug concentration contours could depend to some extent on where the reservoir is precisely placed, relative to the corneocyte or lipid bilayer at the surface. This can be visualized with a mesoscale model (Figure 1). However, in the current lumped approach with anisotropic diffusion in the SC layer, we receive an average drug uptake profile. This is justified if we assume that for a complete patch, the reservoirs are randomly located on the skin by which the average uptake we simulate is still representative. Note that the trend we see for these very small reservoirs is already present for the reservoirs that are much larger than the corneocytes (Figure 14). This implies that smaller reservoirs progressively take more benefit more of the transverse diffusion, compared to larger reservoirs, to increase the drug uptake flux.

FIGURE 15

Figure 15 Color contours of drug concentration over the skin for different sizes of the reservoir for simulations with an infinite reservoir after 72 h, so when a steady-state is reached. Note that the maximal range is not depicted (80 kg m^-3), and only one-fourth of the patch-skin system shown due to symmetry.

Outlook

Mechanistic Modeling

The current state of the art in mechanistic modeling of TTDS was summarized in Table 1. Such mechanistic modeling provides several advantages compared to the analytical solution of the diffusion-driven drug uptake process (Rim et al., 2005; Khanday and Rafiq, 2016). Such analytical solutions enable to calculate drug dose taken up for several drugs, based on their diffusive properties of the skin. Mechanistic modeling is, however, required to target more complex situations, for example, when considering the skin as a multi-layer structure (stratum corneum, viable epidermis) or when the patch is replaced so if the boundary conditions change over time. Compared to the current study and state of the art (Table 1), further advancements should be pursued to enhance the realism and accuracy of transdermal drug delivery predictions in terms of (1) the modeled transport processes, (2) the targeted drug delivery system, (3) numerical modeling, and (4) the model parameters. These future targets for model development are detailed below.

Modeled Transport Processes

Concerning transport processes, the physical adsorption of molecules or chemical binding should be included to enhance accuracy. For fentanyl, bioavailability through the skin is very high [e.g., 92% in (McEvoy et al., 2017)] but not equal to 100%. Thus, some of the fentanyl does not reach the systemic circulation, due to absorption in the epidermis or by chemical changes. Including these mechanisms in mechanistic models is not yet a standard practice (Yamaguchi et al., 2007; Pontrelli and De Monte, 2014). The current mechanistic model for transdermal drug delivery is built up for first-generation systems (Prausnitz and Langer, 2008). Thereby, in addition to fentanyl, it can also be calibrated for lipophilic drugs with a small molecular weight, for example, ibuprofen (Tombs et al., 2018), buprenorphine, rotigotine, and rivastigmine (Pastore et al., 2015). For next-generation systems (Bartosova and Bajgar, 2012; Lee et al., 2018) that enable the delivery of larger molecules like insulin, additional processes would need to be included in the current mechanistic model. These systems aim to enhance skin permeability by increasing the driving force for drug uptake via chemical permeation enhancers, iontophoresis, or by disrupting the stratum corneum using microneedles or thermal ablation (Bartosova and Bajgar, 2012).

Furthermore, when including the dermis in the model (Naegel et al., 2011; Selzer et al., 2013; Selzer et al., 2015), which is implicitly done when modeling skin thicknesses in the millimeter range, blood flow (in addition to diffusion) must be modeled. Modeling the dermis without drug extraction by blood flow will underpredict the drug uptake rate. The impact of including the dermis in the modeled system configuration is illustrated in Supplementary Material 4 (Figure S2), where the impact of different dermis thicknesses on diffusive drug transport is illustrated, so without modeling blood flow. Due to the relatively large dermis thickness and volume, modeling only diffusive transport overpredicts the amount of stored drug and transverse spreading found in specific studies (Selzer et al., 2015). The impact of tight junctions (Bäsler et al., 2016) on diffusion was not explicitly modeled but should be considered in the future.

Additionally, skin swelling or shrinkage, for example, driven by changes in skin hydration, could be included (Dąbrowska et al., 2016). For an infinite reservoir under steady-state conditions, where Fickian diffusion would predict a constant flux, swelling will introduce a time-dependency into the flux (Perrie and Rades, 2014). The swelling or shrinkage of the patch could also be considered (Rim et al., 2005). Finally, diffusion and partition coefficients that are a function of drug concentration (rather than constant values) should be used, especially if there are large variations for the drug of interest. For example, when increasing the Azone concentration from 0 to 9 wt%, the partition coefficient changed from -4.0 to 2.1 (Lundborg et al., 2018). This alteration significantly changes its equilibrium distribution. Including all of the aforementioned physical processes in the model is not always required as they don’t all play critical roles for each drug molecule. Their importance should be assessed on a case by case basis.

Drug Delivery System

Concerning the modeled delivery system, finite reservoirs should always be preferred over infinite reservoirs, which are currently still commonplace (Table 1). For finite reservoirs, the gradient, therefore the driving force, will decrease, a phenomenon that makes delivery at a constant rate more challenging (Figure 6). For this reason, controlled drug delivery systems were designed to alleviate the decreasing concentration gradient [e.g., multilayer Deponit® system (Chien and Lin, 2007; Perrie and Rades, 2014)]. Future models should also include the dermis. This inclusion is essential to evaluate the drug share taken up by the blood versus the amount of the drug that diffuses into and is stored in the thick dermis. This factor will affect total bioavailability and uptake kinetics. In this study, the storage in the dermis was assumed small compare to the amount of drug taken up via the dermis in the blood flow, due to the large bioavailability of fentanyl. If the dermis is included in the computational model, it is essential to include the blood flow in capillaries and vessels and to adjust this blood flow as a function of patient age and activity level, among others (Simmons et al., 2011). Another reason to include the dermis is that the current boundary condition imposed at the viable epidermis (surface 1), namely a zero concentration, has a certain limitation. This condition is valid in experimental setups with Franz diffusion cells, but it can be disputed for a real human tissue. Here, the concentration profile at the viable epidermis will result from the trade-off between diffusion and uptake by the blood flow. Finally, the mechanistic model should be linked to a pharmacokinetic model that relates the uptake amount to the metabolization process in the body to obtain the final blood plasma concentration.

Numerical Modeling

Concerning numerical modeling, there were large discontinuities in concentration over the skin layers due to partitioning. To improve numerical stability and accuracy, it is advised to solve for a dependent variable other than concentration, as was performed in the present study using drug potential.

Model Parameters

Concerning the model parameters, the diffusion and partition coefficients are rarely measured explicitly before modeling using a separate in vitro experiment. Instead, data from the literature are utilized, often even from other drugs with similar molecular weight and lipophilicity (Rim et al., 2009). Alternatively, data are also fitted (Rim et al., 2005) or corrected (Naegel et al., 2008) to match experimental data. Obtaining a good agreement, in this case, is not surprising and can mask missing physical processes within the model. Therefore, one cannot claim the model is validated, but rather it should be considered calibrated. This procedure to obtain model parameters is not necessarily discouraged, but one cannot use the same dataset for model fitting and experimental comparison. Selzer et al. (2015) recently fit model parameters to in vitro experiments and compared the calibrated model to a set of in vivo experiments. This approach is certainly viable, but the resulting model parameters still led to differences in the experiments.

Patient Variability and Personalization

The impact of aging on transdermal drug delivery was considered by changing the stratum corneum thickness. For a 70-year-old patient, the drug taken up by the blood was lower, and the maximum peak flow of the drug occurred later in time than for an 18-year-old patient. This in-silico data analysis quantified the dose delivered for differently aged patients for a specific drug concentration in the patch. Moreover, the slight time shift in peak drug flow is also helpful for defining the time window where patients of different ages may be more susceptible to developing side effects. This information enables more precise monitoring and prevention of serious side effects. These data could be used to tailor devices for specific age groups or to develop devices that can monitor the drug flux and tailor it depending on the person’s age. This would be a step forward compared to current conventional transdermal fentanyl therapy. Where the initial dose (so patch size) is being estimated based on previous daily doses of oral morphine for the patient (McPherson, 2010), with applying the patch transdermally and replacing it every 72 h (Muijsers and Wagstaff, 2001). These features would allow a tailored treatment and a constant delivery rate below defined thresholds (Lee et al., 2018). Additional age-related changes in the stratum corneum structure, such as a decrease with lipid content and its composition (Rogers et al., 1996), lipid peroxidation (Kammeyer and Luiten, 2015; Yadav et al., 2019), or reduced hydration level (Farage et al., 2007; Boireau-Adamezyk et al., 2014), were not accounted for in this study. There are several challenges with designing such tailored devices or therapy for transdermal drug delivery and implementing then into the clinics. The most straightforward solution that could be implemented the most swiftly would be to use conventional transdermal therapy based on existing fentanyl patches. The most optimal therapy with respect to the initial concentration in the patch, the amount of time it should be applied (currently 72 h), and the location where it should be applied could be determined per age category. This means that a clinician could use the mechanistic model to decide which patch to use (e.g., Durogesic® DTrans® 12 µg h^-1 or 25 µg h^-1), how long to use it (72 h or shorter/longer) and where to place it on the body.

Before undertaking clinical trials in vivo, a detailed in vitro study is still required. The model and its transport processes were validated already using Franz diffusion cells (section Validation). Nevertheless, the findings with respect to age and anatomical location need to be confirmed experimentally to further consolidate the findings out of the present study. Such an in vitro study is also an essential next step to further widen the possible use of the numerical model to help designing fentanyl therapy.

For now, we only discussed the “individualization” of therapy based on patients with different age categories. In reality, patients within this individual category will also have a certain variability in skin build-up, so drug uptake. This interpatient variability is an additional point to tackle when completely individualizing therapy for every single patient. The mechanistic model can in principle, calculate the individualized drug uptake for every single patient, but this would require the right input parameters, such as skin thickness, which can be challenging to estimate for each patient in a clinical context.

Reservoir Design and Contact Area

From our results, multiple smaller isolated reservoirs with the same total exposed surface area will likely release a drug at a higher rate than a single reservoir. This finding suggests that drug uptake could be enhanced only by patch design, namely by taking advantage of the transverse drug transport in the epidermis and particularly in the stratum corneum, without penetrating the stratum corneum layer with microneedles or ablating it. This finding could help in designing and individualizing TDDS by customizing the reservoir size and, thus, contact surface area to the patient. Identifying and quantifying this effect is only possible with simulations because accurately measuring these local fluxes over small reservoirs (e.g., 2 µm) would be very challenging experimentally. Although the flux increases with smaller contact surfaces, the total amount of drug delivered (mg) will be lower, and these smaller reservoirs will also deplete faster (see Figure 13A). The latter is manifested by a rapidly changing magnitude of the flux over time for finite reservoirs. As such, for smaller reservoirs, it will be more challenging to deliver drugs at a constant rate, which is a key target for TDDS. However, the commercial application of such small reservoirs implies that multiple reservoirs could be integrated into a single patch with a larger thickness or a buffer reservoir to mitigate fast depletion. This design may mitigate the faster depletion and decreased rate, but this concept should be evaluated in more detail in the future to render these findings more conclusive.

Based on the predicted fluxes in Figure 14, we calculate how many small reservoirs are required to achieve the same flow rate of the drug through the skin in order to achieve the same therapeutic effect as a normal patch (L_pt = 40,000 μm). For reservoirs of 400 μm placed on a patch of 16 cm², 4.7 x 10³ reservoirs are required. These reservoirs would only take up 47% of the total surface area of the patch. For reservoirs of 4 μm placed on a patch of 16 cm², 7.4 x 10⁵ reservoirs are required. These reservoirs would take up less than 1% of the total surface area of the patch. One could therefore also choose to design patches much smaller than the conventional patches, when using these smaller reservoirs. However, more detailed analysis and simulations are required to determine the possible interaction between individual reservoirs that are too close in each other’s vicinity, to determine the optimal spacing.

Note that in this study, complete contact between skin and drug reservoir was assumed. This phenomenon implies that the material in the drug reservoir (e.g., a gel) is sufficiently deformable to ensure such perfect contact upon application. If the contact is not perfect, the skin roughness and the hydration level will also play a role in drug diffusion and delivery and thus must be explicitly considered (Dąbrowska et al., 2016).

Future Use of Mechanistic Modeling in TDDS

A decisive factor in the future use of mechanistic modeling for TDDS design is to increase their efficiency. A key bottleneck in the workflow is typically the low availability of accurate and suitable model parameters and combining in vivo data from different studies, a factor that lowers the accuracy of the model solution (Selzer et al., 2015; Keurentjes and Maibach, 2019). Currently, it is possible to obtain a satisfactory agreement between mechanistic model results and experimental data, but differences remain during the uptake process (Rim et al., 2005; Selzer et al., 2013; Gajula et al., 2017). A better absolute agreement could be pursued by obtaining more accurate model parameters from detailed experimental measurements, for example, the a priori in vitro parameter determination that was applied in one study (Selzer et al., 2015), but even then some discrepancies with between simulations and experiments remained. However, if this time-intensive step needs to be performed for every new set up or drug of interest, in silico TDDS design will not be used beyond academic studies. This endeavor would require more time to obtain the model parameters for simulations rather than to perform stand-alone experiments. In this case, simulations would likely be of primary interest when a single drug is considered, and a very large parametric space is explored. An alternative for obtaining the model parameters from experiments would be inverse modeling.

Instead of pursuing a high precision in the predicted uptake kinetics (Selzer et al., 2015), computer-aided engineering (CAE) in TDDS could, however, easily be used to probe for relative differences among systems, devices, drugs, or patients. For this purpose, model parameter data can be simply based on the literature, a feature that would make the entire modeling suite much swifter and more attractive. With this perspective, CAE has been used to push innovation in many other engineering fields, from modeling blood flow in vessels to the integrated design and construction of civil structures (Sanderse and Weippl, 2018).

Finally, mechanistic modeling could become an essential component of fourth-generation controlled, feedback-induced TDDS. These TDDS enable multiple drug release rates and use feedback control based on measured biomarkers via wearable sensors to regulate drug delivery (Lee et al., 2018). Currently, such wearable technologies are being developed to individualize treatment (Lee et al., 2018). Wearable sensors measure biomarkers and thus, monitor the patient’s physiological condition, which is used to trigger the release of medication. Consequently, actuated transdermal drug delivery patches steer the drug release to provide the correct dose. The sensors and actuators are integrated in a closed feedback loop for better individualizing therapy to deliver the optimal rate at the right time to the correct body location. These TDDS rely on the measured patient’s response, i.e. the change of certain biomarkers, for control. There can be, however, a significant time lag between the drug release and the drugs reaching the systemic circulation. Mechanistic simulations could provide key complementary quantitative insight on the drug-release and percutaneous absorption kinetics for each patient, and help to better control fourth-generation TDDS.

In this context, a next step in these fourth-generation TDDS could be the use of digital twins. A digital twin is a virtual representation of the TDDS, which is linked to the real-world patient by sensor data of certain biomarkers. Digital twins can be used here for predictive modeling of the drug release and uptake in the human body, and can quantify the time lag, and account for it during control. Here, mechanistic modeling is an essential building block for digital twins of human organs.

Conclusions

Validated mechanistic modeling was used to gain new insights into transdermal fentanyl uptake. First, we quantified the changes in transdermal fentanyl uptake with the patient’s age and the anatomical location where the patch was placed. We also evaluated how much the drug flux can be enhanced by miniaturizing the drug reservoir size. Additionally, we obtained quantitative insights into the release and uptake kinetics of fentanyl transdermal drug delivery by analyzing drug diffusion, storage, and partitioning. The main findings are summarized below.

-Differences in drug uptake amount between anatomical locations of the drug reservoir on the human body of 36% after 72 h were found, but there was also strong interpatient variability.

-With aging, the transdermal drug delivery patch worked slower and less potently. An 18-year-old patient received 26% more drugs over the 72 h application period than a 70-year-old patient.

-Our proposed novel concept of using micron-sized drug reservoirs induced a much higher flux (µg cm^-2 h^-1) than conventional patches. Due to enhancing transverse diffusion in the stratum corneum layer, simply by changing patch design, a 200-fold increase in the drug flux was possible for a micron-sized patch. The identification of fluxes at the micron scale in a straightforward and accurate way was only possible in silico. With this concept, the reservoir surface area can be tuned and individualized for a specific patient or patient with a certain age category.

-For commercial patches, where the size is much larger than the epidermal thickness, drug transport was mainly unidirectional, and 1D models can be used reliably, assuming perfect contact between skin and reservoir.

-We showed in silico that sampling intervals in experiments > 5 h significantly underestimated peak drug fluxes, for example, by 20% when considering a 10-h interval.

-The role of transverse diffusion, particularly in the stratum corneum, was strongly dependent on the patch size and did not play a critical role in conventional patches.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

TD and RR conceptualized the study and acquired funding. TD did the project administration. TD performed the investigation, developed the methodology, performed the validation, and executed the simulations with key input from LD, who did exploratory simulation work during her MSc thesis. TD and AT supervised LD and FB. TD wrote the original draft of the paper and did the visualization, with key input from RM and FB. RR, FB, LD, RM, and AT performed critical review and editing.

Funding

This work was supported by the Novartis Research Foundation (grant “Virtual twinning for intelligent, personalized transdermal drug delivery”). The authors declare that this study received funding from Novartis Research Foundation. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication. This manuscript has been released as a pre-print at bioRxiv (Defraeye et al., 2020).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2020.585393/full#supplementary-material

References

Amarah, A. A., Petlin, D. G., Grice, J. E., Hadgraft, J., Roberts, M. S., Anissimov, Y. G. (2018). Compartmental modeling of skin transport. Eur. J. Pharm. Biopharm. 130, 336–344. doi: 10.1016/j.ejpb.2018.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Amjadi, M., Sheykhansari, S., Nelson, B. J., Sitti, M. (2018). Recent Advances in Wearable Transdermal Delivery Systems. Adv. Mater. 30, 1–19. doi: 10.1002/adma.201704530

CrossRef Full Text | Google Scholar

Andrews, S. N., Jeong, E., Prausnitz, M. R. (2013). Transdermal delivery of molecules is limited by full epidermis, not just stratum corneum. Pharm. Res. 30, 1099–1109. doi: 10.1007/s11095-012-0946-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartosova, L., Bajgar, J. (2012). Transdermal Drug Delivery In Vitro Using Diffusion Cells. Curr. Med. Chem. 19, 4671–4677. doi: 10.2174/092986712803306358

PubMed Abstract | CrossRef Full Text | Google Scholar

Bäsler, K., Bergmann, S., Heisig, M., Naegel, A., Zorn-Kruppa, M., Brandner, J. M. (2016). The role of tight junctions in skin barrier function and dermal absorption. J. Control. Release 242, 105–118. doi: 10.1016/j.jconrel.2016.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Boireau-Adamezyk, E., Baillet-Guffroy, A., Stamatas, G. N. (2014). Age-dependent changes in stratum corneum barrier function. Ski. Res. Technol. 20, 409–415. doi: 10.1111/srt.12132

CrossRef Full Text | Google Scholar

Bond, J. R., Barry, B. W. (1988). Limitations of Hairless Mouse Skin as a Model for In Vitro Permeation Studies Through Human Skin: Hydration Damage. J. Invest. Dermatol. 90, 486–489. doi: 10.1111/1523-1747.ep12460958

PubMed Abstract | CrossRef Full Text | Google Scholar

Breitbart, W., Chandler, S., Eagel, B., Ellison, N., Enck, R. E., Lefkowitz, M., et al. (2000). An alternative algorithm for dosing transdermal fentanyl for cancer-related pain. Oncol. (Williston Park). 14, 695–705.

Google Scholar

Brockow, K., Garvey, L. H., Aberer, W., Tanaskovic-Markovic, M., Barbaud, A., Bilo, M. B., et al. (2013). Skin test concentrations for systemically administered drugs - An ENDA/EAACI Drug Allergy Interest Group position paper. Allergy Eur. J. Allergy Clin. Immunol. 68, 702–712. doi: 10.1111/all.12142

CrossRef Full Text | Google Scholar

Calcutt, J. J., Anissimov, Y. G. (2019). Physiologically based mathematical modelling of solute transport within the epidermis and dermis. Int. J. Pharm. 569, 118547. doi: 10.1016/j.ijpharm.2019.118547

PubMed Abstract | CrossRef Full Text | Google Scholar

Casey, M., Wintergerste, T. (2000). Special Interst Group on “Quality and Trust in Industrial CFD” Best Practice Guidelines. 1st ed. (ERCOFTAC).

Google Scholar

Chen, L., Han, L., Saib, O., Lian, G. (2015). In Silico Prediction of Percutaneous Absorption and Disposition Kinetics of Chemicals. Pharmaceutical 32, 1779–1793. doi: 10.1007/s11095-014-1575-0

CrossRef Full Text | Google Scholar

Chien, Y. W., Lin, S. (2007). “Drug Delivery : Controlled Release,” in Encyclopedia of Pharmaceutical Technology Informa Healthcare, 3rd ed. Ed. Swarbrick, J. (New York), 2092.

Google Scholar

Dąbrowska, A. K., Spano, F., Derler, S., Adlhart, C., Spencer, N. D., Rossi, R. M., et al. (2016). In vivo confirmation of hydration-induced changes in human-skin thickness, roughness and interaction with the environment. Biointerphases 11, 031015. doi: 10.1116/1.4962547

PubMed Abstract | CrossRef Full Text | Google Scholar

Dąbrowska, A. K., Spano, F., Derler, S., Adlhart, C., Spencer, N. D., Rossi, R. M. (2018). The relationship between skin function, barrier properties, and body-dependent factors. Ski. Res. Technol. 24, 165–174. doi: 10.1111/srt.12424

CrossRef Full Text | Google Scholar

Das, C., Olmsted, P. D. (2016). The physics of stratum corneum lipid membranes. Philos. Trans. R. Soc A Math. Phys. Eng. Sci. 374, 20150126. doi: 10.1098/rsta.2015.0126

CrossRef Full Text | Google Scholar

Datta, A. K. (2007). Porous media approaches to studying simultaneous heat and mass transfer in food processes. II: Property data and representative results. J. Food Eng. 80, 96–110. doi: 10.1016/j.jfoodeng.2006.05.012

CrossRef Full Text | Google Scholar

de Monte, F., Pontrelli, G., Becker, S. M. (2015). Transdermal Drug Delivery and Percutaneous Absorption: Mathematical Modeling Perspectives (Elsevier Inc.).

Google Scholar

Defraeye, T., Verboven, P. (2017). Convective drying of fruit: Role and impact of moisture transport properties in modelling. J. Food Eng. 193, 95–107. doi: 10.1016/j.jfoodeng.2016.08.013

CrossRef Full Text | Google Scholar

Defraeye, T., Blocken, B., Carmeliet, J. (2012). Analysis of convective heat and mass transfer coefficients for convective drying of a porous flat plate by conjugate modelling. Int. J. Heat Mass Transf. 55, 1–3112–124. doi: 10.1016/j.ijheatmasstransfer.2011.08.047

CrossRef Full Text | Google Scholar

Defraeye, T., Bahrami, F., Ding, L., Malini, R., II, Terrier, A., Rossi, R. M. (2020). “Predicting transdermal fentanyl delivery using mechanistic simulations for tailored therapy,” bioRxiv p, 2020.06.16.154195. doi: 10.1101/2020.06.16.154195

CrossRef Full Text | Google Scholar

EMC (2019). electronic Medicines Compendium. Available at: https://www.medicines.org.uk/emc (Accessed Accessed: 03-Jan-2019).

Google Scholar

Farage, M. A., Miller, K. W., Elsner, P., Maibach, H., II (2007). Structural characteristics of the aging skin: A review. Cutan. Ocul. Toxicol. 26, 343–357. doi: 10.1080/15569520701622951

PubMed Abstract | CrossRef Full Text | Google Scholar

Farahmand, S., Maibach, H., II (2009). Transdermal drug pharmacokinetics in man: Interindividual variability and partial prediction. Int. J. Pharm. 367, 1–15. doi: 10.1016/j.ijpharm.2008.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Feldmann, R. J., Maibach, H., II (1967). Regional variation in percutaneous penetration of 14C cortisol in man. J. Invest. Dermatol. 48, 181–183. doi: 10.1038/jid.1967.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Franke, B., Hellsten, A., Schlünzen, H., Carissimo (2007). Best practice guideline for the CFD simulation of flows in the urban environment (Hamburg).

Google Scholar

Gajula, K., Gupta, R., Sridhar, D. B., Rai, B. (2017). In-Silico Skin Model: A Multiscale Simulation Study of Drug Transport. J. Chem. Inf. Model. 57, 2027–2034. doi: 10.1021/acs.jcim.7b00224

PubMed Abstract | CrossRef Full Text | Google Scholar

Garg, A., Singh, S. (2018). “Transdermal drug delivery systems,” in In-Vitro and In-Vivo Tools in Drug Delivery Research for Optimum Clinical Outcomes, 1st ed., vol. 3. Eds. Misra, A., Shahiwala, A. (Boca Raton: CRC Press), 51–77.

Google Scholar

Grassi, M., Lamberti, G., Cascone, S., Grassi, G. (2011). Mathematical modeling of simultaneous drug release and in vivo absorption. Int. J. Pharm. 418, 130–141. doi: 10.1016/j.ijpharm.2010.12.044

PubMed Abstract | CrossRef Full Text | Google Scholar

GVR (2016). Transdermal Drug Delivery System Market Analysis By Technology (Electroporation, Radiofrequency, Iontophoresis, Microporation, Thermal, Mechanical Arrays, Ultrasound), By Application (Central Nervous System (CNS), Pain Management, Cardiovascular, Hormone).

Google Scholar

Hansen, S., Henning, A., Naegel, A., Heisig, M., Wittum, G., Neumann, D., et al. (2008). 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. 68, 352–367. doi: 10.1016/j.ejpb.2007.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Heikkinen, E. M., Voipio, H. M., Laaksonen, S., Haapala, L., Räsänen, J., Acharya, G., et al. (2015). Fentanyl Pharmacokinetics in Pregnant Sheep after Intravenous and Transdermal Administration to the Ewe. Basic Clin. Pharmacol. Toxicol. 117, 156–163. doi: 10.1111/bcpt.12382

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, D. J., Maibach, H., II (1990). Adverse dermatologic reactions to transdermal drug delivery systems. J. Am. Acad. Dermatol. 22, 811–814. doi: 10.1016/0190-9622(90)70113-V

PubMed Abstract | CrossRef Full Text | Google Scholar

Kammeyer, A., Luiten, R. M. (2015). Oxidation events and skin aging. Ageing Res. Rev. 21, 16–29. doi: 10.1016/j.arr.2015.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Keurentjes, A. J., Maibach, H., II (2019). Percutaneous penetration of drugs applied in transdermal delivery systems: an in vivo based approach for evaluating computer generated penetration models. Regul. Toxicol. Pharmacol. 108, 104428. doi: 10.1016/j.yrtph.2019.104428

PubMed Abstract | CrossRef Full Text | Google Scholar

Khanday, M. A., Rafiq, A. (2016). Numerical Estimation of Drug Diffusion At Dermal Regions of Human Body in Transdermal Drug Delivery System. J. Mech. Med. Biol. 16, 1650022. doi: 10.1142/S0219519416500226

CrossRef Full Text | Google Scholar

Kim, S., Thiessen, P. A., Bolton, E. E., Chen, J., Fu, G., Gindulyte, A., et al. (2016). PubChem Substance and Compound databases. Nucleic Acids Res. 44, 1202–1213. doi: 10.1093/nar/gkv951

CrossRef Full Text | Google Scholar

Kirschner, N., Houdek, P., Fromm, M., Moll, I., Brandner, J. M. (2010). Tight junctions form a barrier in human epidermis. Eur. J. Cell Biol. 89, 839–842. doi: 10.1016/j.ejcb.2010.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreienbuehl, A., Naegel, A., Ruprecht, D., Speck, R., Wittum, G., Krause, R. (2015). Numerical simulation of skin transport using Parareal. Comput. Vis. Sci. 17, 99–108. doi: 10.1007/s00791-015-0246-y

CrossRef Full Text | Google Scholar

Kress, H. G., Boss, H., Delvin, T., Lahu, G., Lophaven, S., Marx, M., et al. (2010). Transdermal fentanyl matrix patches Matrifen® and Durogesic® DTrans® are bioequivalent. Eur. J. Pharm. Biopharm. 75, 225–231. doi: 10.1016/j.ejpb.2010.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Lademann, J., Jacobi, U., Surber, C., Weigmann, H. J., Fluhr, J. W. (2009). The tape stripping procedure - evaluation of some critical parameters. Eur. J. Pharm. Biopharm. 72, 317–323. doi: 10.1016/j.ejpb.2008.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, R. H., Nielsen, F., Sørensen, J. A., Nielsen, J. B. (2003). Dermal Penetration of Fentanyl : Inter- and Intraindividual Variations. Pharmacol. Toxicol. 93, 244–248.

Google Scholar

Lee, H., Song, C., Baik, S., Kim, D., Hyeon, T., Kim, D. H. (2018). Device-assisted transdermal drug delivery. Adv. Drug Deliv. Rev. 127, 35–45. doi: 10.1016/j.addr.2017.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Leopold, C. S., Maibach, H., II (1996). journal of Effect of lipophilic vehicles on in vivo skin penetration of methyl nicotinate in different races. Int. J. Pharm. 5173, 1–7. doi: 10.1016/0378-5173(96)04562-0

CrossRef Full Text | Google Scholar

Lundborg, M., Wennberg, C. L., Narangifard, A., Lindahl, E., Norlén, L. (2018). Predicting drug permeability through skin using molecular dynamics simulation. J. Control. Release 283, 269–279. doi: 10.1016/j.jconrel.2018.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Manitz, R., Lucht, W., Strehmel, K., Weiner, R., Neubert, R. (1998). On mathematical modeling of dermal and transdermal drug delivery. J. Pharm. Sci. 87, 873–879. doi: 10.1021/js970329r

PubMed Abstract | CrossRef Full Text | Google Scholar

Marier, J. F., Lor, M., Potvin, D., DiMarco, M., Morelli, G., Sædder, E. A. (2006). Pharmacokinetics, tolerability, and performance of a novel matrix transdermal delivery system of fentanyl relative to the commercially available reservoir formulation in healthy subjects. J. Clin. Pharmacol. 46, 642–653. doi: 10.1177/0091270006286901

PubMed Abstract | CrossRef Full Text | Google Scholar

Markets, M. (2018). Transdermal Drug Delivery System Market by Type (Patches and Semisolid formulations), Applications (Pain Management, Central Nervous System Disorders, Hormonal Applications, Cardiovascular Diseases), End User, and Region - Global Forecast to 2023.

Google Scholar

Matsui, T., Amagai, M. (2015). Dissecting the formation, structure and barrier function of the stratum corneum. Int. Immunol. 27, 269–280. doi: 10.1093/intimm/dxv013

PubMed Abstract | CrossRef Full Text | Google Scholar

McEvoy, G. K., Litvak, K., Welsh, O. H., Jr. (2017). American Society of Health-System Pharmacists: Drug Information 2017 (Bethesda, Md: American Hospital Formulary Service).

Google Scholar

McPherson, M. L. (2010). Transdermal and parenteral fentanyl dosage calculations and conversions. Bethesda Am. Soc Heal. Pharm.

Google Scholar

Menczel, E., Maibach, H., II (1970). In vitro human percutaneous penetration of benzyl alcohol and testosterone: epidermal-dermal retention. J. Invest. Dermatol. 54, 386–394. doi: 10.1111/1523-1747.ep12259143

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitragotri, S., Anissimov, Y. G., Bunge, A. L., Frasch, H. F., Guy, R. H., Hadgraft, J., et al. (2011). Mathematical models of skin permeability : An overview. Int. J. Pharm. 418, 115–129. doi: 10.1016/j.ijpharm.2011.02.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Moon, J. M., Chun, B. J. (2011). Fentanyl intoxication caused by abuse of transdermal fentanyl. J. Emerg. Med. 40, 37–40. doi: 10.1016/j.jemermed.2007.10.075

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrison, T. M., F. D. A. Guidance. (2016). Reporting of Computational Modeling Studies in Medical Device Submissions - Guidance for Industry and Food and Drug Administration Staff. (US Department of Health and Human Services Food and Drug Administration), 21.

Google Scholar

Muijsers, R. B. R., Wagstaff, A. J. (2001). Transdermal fentanyl: An updated review of its pharmacological properties and therapeutic efficacy in chronic cancer pain control. Drugs 61, 2289–2307. doi: 10.2165/00003495-200161150-00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Naegel, A., Hansen, S., Neumann, D., Lehr, C. M., Schaefer, U. F., Wittum, G., et al. (2008). In-silico model of skin penetration based on experimentally determined input parameters. Part II: Mathematical modelling of in-vitro diffusion experiments. Identification of critical input parameters. Eur. J. Pharm. Biopharm. 68, 368–379. doi: 10.1016/j.ejpb.2007.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Naegel, A., Heisig, M., Wittum, G. (2009). A comparison of two- and three-dimensional models for the simulation of the permeability of human stratum corneum. Eur. J. Pharm. Biopharm. 72, 332–338. doi: 10.1016/j.ejpb.2008.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Naegel, A., Hahn, T., Schaefer, U. F., Lehr, C. M., Heisig, M., Wittum, G. (2011). Finite dose skin penetration: A comparison of concentration-depth profiles from experiment and simulation. Comput. Vis. Sci. 14, 327–339. doi: 10.1007/s00791-012-0186-8

CrossRef Full Text | Google Scholar

Naegel, A., Heisig, M., Wittum, G. (2013). Detailed modeling of skin penetration-An overview. Adv. Drug Deliv. Rev. 65, 191–207. doi: 10.1016/j.addr.2012.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastore, M. N., Kalia, Y. N., Horstmann, M., Roberts, M. S. (2015). Transdermal patches: History, development and pharmacology. Br. J. Pharmacol. 172, 2179–2209. doi: 10.1111/bph.13059

PubMed Abstract | CrossRef Full Text | Google Scholar

Perrie, Y., Rades, T. (2014). Pharmaceutics - Drug Delivery and Targeting. 2nd ed. (London: Pharmaceutical Press).

Google Scholar

Pontrelli, G., De Monte, F. (2014). Mathematical Biosciences A two-phase two-layer model for transdermal drug delivery and percutaneous absorption. Math. Biosci. 257, 96–103. doi: 10.1016/j.mbs.2014.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Prausnitz, M. R. M., Langer, R. (2008). Transdermal drug delivery. Nat. Biotechnol. 26, 1261–1268. doi: 10.1038/nbt.1504

PubMed Abstract | CrossRef Full Text | Google Scholar

Rim, J. E., Pinsky, P. M., Van Osdol, W. W. (2005). Finite element modeling of coupled diffusion with partitioning in transdermal drug delivery. Ann. Biomed. Eng. 33, 1422–1438. doi: 10.1007/s10439-005-5788-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Rim, J. E., Pinsky, P. M., van Osdol, W. W. (2008). Using the method of homogenization to calculate the effective diffusivity of the stratum corneum with permeable corneocytes. J. Biomech. 41, 788–796. doi: 10.1016/j.jbiomech.2007.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Rim, J. E., Pinsky, P. M., Van Osdol, W. W. (2009). Multiscale modeling framework of transdermal drug delivery. Ann. Biomed. Eng. 37, 1217–1229. doi: 10.1007/s10439-009-9678-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Roache, P. J. (1994). Perspective: A Method for Uniform Reporting of Grid Refinement Studies. J. Fluids Eng. 116, 405–413. doi: 10.1115/1.2910291

CrossRef Full Text | Google Scholar

Rogers, J., HArding, C., Mayo, A., Banks, J., Rawlings, A. (1996). Stratum corneum lipids: the effect of aging and seasons. Arch. Dermatol. Res. 288, 765–770. doi: 10.1007/BF02505294

PubMed Abstract | CrossRef Full Text | Google Scholar

Roskos, K. V., Maibach, H., II, Guy, R. H. (1989). The Effect of Aging on Percutaneous Absorption in Man. J. Pharmacokinet. Biopharm. 17, 617–630. doi: 10.1007/BF01062121

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, S. D., Flynn, G. L. (1990). Transdermal Delivery of Narcotic Analgesics: pH, Anatomical, and Subject Influences on Cutaneous Permeability of Fentanyl and Sufentanil. Pharm. Res. 7 (8), 842–847. doi: 10.1023/A:1015912932416

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandby-Møller, J., Poulsen, T., Wulf, H. C. (2003). Epidermal Thickness at Different Body Sites: Relationship to Age, Gender, Pigmentation, Blood Content, Skin Type and Smoking Habits. Acta Derm. Venereol. 83, 410–413. doi: 10.1080/00015550310015419

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanderse, B., Weippl, E. (2018). Digital twins. ERCIM News 114, 1–56.

Google Scholar

Schmid-Grendelmeier, P., Pokorny, R., Gasser, U. E., Richarz, U. (2006). A comparison of the skin irritation potential of transdermal fentanyl versus transdermal buprenorphine in middle-aged to elderly healthy volunteers. Curr. Med. Res. Opin. 22, 501–509. doi: 10.1185/030079906X89829

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwöbel, J. A. H., Klamt, A. (2019). Mechanistic skin penetration model by the COSMO perm method : Routes of permeation, vehicle effects and skin variations in the healthy and compromised skin. Comput. Toxicol. 11, 50–64. doi: 10.1016/j.comtox.2019.02.004

CrossRef Full Text | Google Scholar

Selzer, D., Abdel-Mottaleb, M. M., Hahn, T., Schaefer, U. F., Neumann, D. (2013). Finite dose skin mass balance including the lateral part: Comparison between experiment, pharmacokinetic modeling and diffusion models. J. Control. Release 165, 119–128. doi: 10.1016/j.jconrel.2012.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Selzer, D., Neumann, D., Neumann, H., Kostka, K., Lehr, C., Schaefer, U. F. (2015). European Journal of Pharmaceutics and Biopharmaceutics A strategy for in-silico prediction of skin absorption in man. Eur. J. Pharm. Biopharm. 95, 68–76. doi: 10.1016/j.ejpb.2015.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Siepmann, J., Siepmann, F. (2008). Mathematical modeling of drug delivery. Int. J. Pharm. 364, 328–343. doi: 10.1016/j.ijpharm.2008.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Simmons, G. H., Wong, B. J., Holowatz, L. A., Kenney, W. L. (2011). Changes in the control of skin blood flow with exercise training: where do cutaneous vascular adaptations fit in? Exp. Physiol. 96, 822–828. doi: 10.1113/expphysiol.2010.056176

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, I., Morris, A. (2011). Performance of transdermal therapeutic systems: Effects of biological factors. Int. J. Pharm. Investig. 1, 4. doi: 10.4103/2230-973X.76721

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanley, T. H. (2014). The Fentanyl Story. J. Pain 15, 1215–1226. doi: 10.1016/j.jpain.2014.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Thakur, R., Batheja, P., Kaushik, D., Michniak, B. (2009). “Structural and Biochemical Changes in Aging Skin and Their Impact on Skin Permeability Barrier,” in Skin Aging Handbook (William Andrew Inc.), pp. 55–90.

Google Scholar

Tombs, E. L., Nikolaou, V., Nurumbetov, G., Haddleton, D. M. (2018). Transdermal Delivery of Ibuprofen Utilizing a Novel Solvent-Free Pressure-sensitive Adhesive (PSA): TEPI® Technology. J. Pharm. Innov. 13, 48–57. doi: 10.1007/s12247-017-9305-x

PubMed Abstract | CrossRef Full Text | Google Scholar

US Food and Drug Administration (2005).

Google Scholar

Vieille-Petit, A., Blickenstaff, N., Coman, G., Maibach, H. (2015). Metrics and clinical relevance of percutaneous penetration and lateral spreading. Skin Pharmacol. Physiol. 28, 57–64. doi: 10.1159/000363148

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiedersberg, S., Guy, R. H. (2014). Transdermal drug delivery: 30 + years of war and still fighting! J. Control. Release 190, 150–156. doi: 10.1016/j.jconrel.2014.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Wishart, D. S., Feunang, Y. D., Guo, A. C., Lo, E. J., Marcu, A., Grant, J. R., et al. (2018). DrugBank 5.0 : a major update to the DrugBank database for 2018. Nucleic Acids Res. 46, 1074–1082. doi: 10.1093/nar/gkx1037

CrossRef Full Text | Google Scholar

Wittum, R., Naegel, A., Heisig, M., Wittum, G. (2017). Mathematical modelling of the viable epidermis: impact of cell shape and vertical arrangement. Math. Mech. Solids pp1–14, 1046–1059. doi: 10.1177/1081286517743297

CrossRef Full Text | Google Scholar

Woodhouse, A., Mather, L. E. (2000). The minimum effective concentration of opioids: a revisitation with patient controlled analgesia fentanyl. Reg. Anesth. Pain Med. 25, 259–267. doi: 10.1097/00115550-200005000-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, D. K., Kumar, S., Choi, E. H., Chaudhary, S., Kim, M. H. (2019). Molecular dynamic simulations of oxidized skin lipid bilayer and permeability of reactive oxygen species. Sci. Rep. 9, 1–10. doi: 10.1038/s41598-019-40913-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, K., Mitsui, T., Aso, Y., Sugibayashi, K. (2007). Analysis of in vitro skin permeation of 22-oxacalcitriol from ointments based on a two- or three-layer diffusion model considering diffusivity in a vehicle. Int. J. Pharm. 336, 310–318. doi: 10.1016/j.ijpharm.2006.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Glossary.