Inverse Mechanistic Modeling of Transdermal Drug Delivery for Fast Identification of Optimal Model Parameters

Transdermal drug delivery systems are a key technology to administer drugs with a high first-pass effect in a non-invasive and controlled way. Physics-based modeling and simulation are on their way to become a cornerstone in the engineering of these healthcare devices since it provides a unique complementarity to experimental data and additional insights. Simulations enable to virtually probe the drug transport inside the skin at each point in time and space. However, the tedious experimental or numerical determination of material properties currently forms a bottleneck in the modeling workflow. We show that multiparameter inverse modeling to determine the drug diffusion and partition coefficients is a fast and reliable alternative. We demonstrate this strategy for transdermal delivery of fentanyl. We found that inverse modeling reduced the normalized root mean square deviation of the measured drug uptake flux from 26 to 9%, when compared to the experimental measurement of all skin properties. We found that this improved agreement with experiments was only possible if the diffusion in the reservoir holding the drug was smaller than the experimentally measured diffusion coefficients suggested. For indirect inverse modeling, which systematically explores the entire parametric space, 30,000 simulations were required. By relying on direct inverse modeling, we reduced the number of simulations to be performed to only 300, so a factor 100 difference. The modeling approach’s added value is that it can be calibrated once in-silico for all model parameters simultaneously by solely relying on a single measurement of the drug uptake flux evolution over time. We showed that this calibrated model could accurately be used to simulate transdermal patches with other drug doses. We showed that inverse modeling is a fast way to build up an accurate mechanistic model for drug delivery. This strategy opens the door to clinically ready therapy that is tailored to patients.

Introduction

Mechanistic, physics-based modeling and simulation is increasingly used in healthcare for device engineering and design (FDA, 2016; Caccavo et al., 2017; Shirazi et al., 2019). A typical example is the transdermal drug delivery systems (TDDS). Physics-based modeling was used to help design and optimize advanced TDDS with microneedles (Ronnander et al., 2018; Chavoshi et al., 2019), thermal skin ablation or iontophoresis (Ferreira et al., 2017; Filipovic et al., 2017), chemical enhancers (Pontrelli and De Monte, 2014), and sonophoresis (Polat et al., 2012). Also, first-generation TDDS are analyzed in-silico to increase our insight into how these systems interact with the patient’s skin (Naegel et al., 2013). Here, drug uptake differences between different substances (Chen et al., 2015) or age categories of patients (Defraeye et al., 2020) were targeted. Also, the impact of the composition of the stratum corneum and viable epidermis on the drug uptake was quantified (Naegel et al., 2009; Wittum et al., 2017).

A key advantage of such first-principles-based modeling is the high spatial and temporal resolution they offer. Transient changes in concentration profiles can be identified within different skin layers during the uptake process (Naegel et al., 2011; Defraeye et al., 2020). This spatiotemporal resolution enables us to virtually probe the drug transport inside and out of the skin in time and space. Also, the cumulative amount of drugs that are taken up by blood flow is available continuously, instead of only every 1–5 h (Defraeye et al., 2020). Another advantage is that such in-silico tests can be performed swiftly and non-invasively, without any side effects for the patient. Physics-based modeling opens new ways of tailoring TDD therapy for patients or patient groups (Defraeye et al., 2020). This deterministic approach is not suffering from statistical variability in the data due to biological variability or measurement uncertainty. As such, only one in-silico experiment is required per case. Thereby, also very small differences between several cases can be accurately identified. Physics-based modeling and simulation provide valuable complementary information to experiments. Besides, they provide useful insights into the relative contribution of different drivers of drug uptake for the processes that can be captured mechanistically. Such physics-based models also provide a higher spatial resolution compared to compartmental models. These compartmental models consider each skin layer as a well-mixed compartment (Mitragotri et al., 2011; Selzer et al., 2013; Amarah et al., 2018) with no discretization over a particular skin layer.

However, the reliability and accuracy of these physics-based models for TDD need to be verified and validated. FDA recently published guidelines for medical device design in this respect (FDA, 2016; U. S. Food and Drug Administration, 2020). A decisive aspect of determining model accuracy is the used model parameters. The properties of the skin, in particular, often play a crucial role, including drug diffusivities and partition coefficients (Naegel et al., 2008).

There are four common ways of obtaining these properties of the skin. First, they are determined directly from ex-vivo measurements (Rim et al., 2005; Hansen et al., 2008). Experimentally determining the diffusion and partition coefficients for a specific drug is resource- and time-consuming, especially when different skin types and anatomical locations on the human body need to be considered. Such an experimental campaign often poses a significant constraint on undertaking a modeling study. Thereby, they are rarely measured explicitly before modeling using a separate in vitro experiment. Second, skin transport properties for various molecules can be simulated with physics-based models (Naegel et al., 2009; Mitragotri et al., 2011). Here simulations at different scales are performed and coupled in a multiscale way (Rim et al., 2009; Gajula et al., 2017). Even molecular dynamics simulations are applied (Lundborg et al., 2018). Simulating, however, also consumes a lot of time and resources, making it is usually not more efficient than measuring the parameters separately. As a third way, modelers just source these parameters from literature data during model setup instead (Barbero and Frasch, 2017; Madhihah et al., 2018; Chakravarty and Dalal, 2019; Defraeye et al., 2020). Due to a lack of appropriate data, often, data of other drugs with similar molecular weight and lipophilicity are taken (Rim et al., 2009). Also, parameters are obtained from multiple sources. As these studies used different setups and boundary conditions (e.g., skin type, temperature, drug reservoirs), the absolute accuracy of the model results is often questioned. This is often the only way for modelers to simulate without a concurrent experimental or simulation campaign to determine material properties. The fourth way is using the quantitative structure-activity relationship (QSAR) modeling method to predict the biological activity of the new drug based on its structural properties (Kwon et al., 2019). This method can provide the stimation of different values for properties of subtances that were not studied experimentally (Toropova and Toropov, 2017), however, the predicted value is limited to the constraints of the implemented dataset (Barratt, 1995).

There is a fifth alternative currently underexplored. All skin material properties could be fitted simultaneously to obtain the best agreement with measured data of the drug uptake kinetics. Calibrating a mathematical model by comparing its predicted response (i.e., drug uptake) with experimental observations is called inverse modeling (Akkaram et al., 2007). This strategy ensures that the in-silico system is calibrated so that it responds as close as possible in the same way as the in vitro (or in vivo) system. A key advantage is that only limited experiments on the drug uptake kinetics need to be performed, instead of performing separate experiments for each property on a different apparatus. A step in this direction was made (Selzer et al., 2015). The diffusion and partition coefficient for the transport of flufenamic acid through the stratum corneum layer were determined by fitting the simulations to the experiments. The in vitro concentration-depth profiles were used to calibrate these material properties. A nonlinear least-squares approach was used to determine these parameters both for every point in time or averaged over all measured times. However, this study did not target the drug uptake amount through the skin by the patient and did not include the transport properties of the patch as possible unknowns that could be fitted.

This study aims at developing and testing a fast and straightforward way to calibrate a physics-based model for transdermal fentanyl delivery. Fentanyl patches for around-the-clock opioid analgesia are currently amongst the most popular transdermal delivery devices (Wiedersberg and Guy, 2014). We target to have an excellent agreement of our model for drug uptake kinetics using a commercial fentanyl patch currently used in the clinics. To this end, a sizable parametric space of the skin epidermis’ material properties and the patch is simulated for transdermal fentanyl uptake through a Franz diffusion cell. A comparison with experimental data on the drug flux out of the epidermis is made. This enables to elucidate in silico the sensitivity of the drug uptake kinetics to each model parameter in the transdermal patch and the skin. Within this parametric space, we determine the best set of properties with a better agreement with experiments during a 72 h drug uptake period, compared to relying on experimentally determined skin properties. A faster alternative procedure to multiparameter inverse modeling by the parametric exploration is also explored, using automated least-squares optimization. We use existing methods for performing this parametric analysis, but aim to show the added value of such methods to gain new insights and to help bring simulation-based insights to the clinics. The validity of the obtained parameter set for simulating other drug patches with a different concentration is also tested. This independent verification proves the generalized use of the model.

Materials and Methods

Continuum Model for Transdermal Fentanyl Delivery

An extensive mechanistic model was set up to simulate fentanyl release from a transdermal patch (reservoir) and subsequent uptake through the human epidermis, which includes the stratum corneum and the viable epidermis (Defraeye et al., 2020). The model and the corresponding simulation were built and executed according to best practice guidelines in modeling, among others, those for medical device design (Casey and Wintergerste, 2000; FDA, 2016). The model was used to simulate the drug uptake experiment described in (Rim et al., 2005). By comparison with the experimental data, we can quantify how accurately the drug uptake kinetics are predicted with our mechanistic model.

Experimental Setup and Data

The experiment involved fentanyl uptake from a cylindrical drug reservoir (radius 9 mm, thickness 50.8 µm) through a cylindrical skin sample (human cadaver epidermis: viable epidermis, and stratum corneum, radius 9.25 mm, thickness 50.8 µm, Figure 1A) into an aqueous solution. A change of the skin thickness due to swelling was not reported in the experiment. The transdermal patch used in the experimental study consisted of a layer of acrylate polymer (Rim et al., 2005). Acrylate polymers are widely used as pressure-sensitive adhesives in TDDs. The acrylate polymers have a good compability with a wide range of drugs. Besides that, they have flexiblility in tailoring the polymers’ properties (Zhan et al., 2015). The drug was embedded in the acrylate polymer, which served as the donor reservoir. Both the reservoir and epidermis were fitted into a Franz diffusion cell and kept at 33°C. The receptor medium in the cell was a phosphate buffer at a pH of 5.65, in which fentanyl had a solubility of 2.5 mg ml⁻¹. The experiments took 72 h. At several points in time over this 3-day period, an aliquot of the receptor medium was removed for this purpose. Its concentration was analyzed via high-performance liquid chromatography (HPLC). Out of the data of the change in concentration and the time, the drug flux was determined in the experiments at discrete points in time. As such, the average flux over a specified time period of multiple hours was obtained. Several initial drug concentrations in the patch (c_pt,ini^α) were evaluated for fentanyl (substance α). We used the data on the drug uptake fluxes that left the epidermis for c_pt,ini^α = 60 kg m⁻³, and c_pt,ini^α = 80 kg m⁻³.

FIGURE 1

FIGURE 1. The geometry of the experimental setup (A) and of the 1D mechanistic model (B) of a cylindrical drug reservoir and epidermis (not to scale).

Computational System Configuration

The geometrical setup is depicted in Figure 1B, along with the boundary conditions. The system configuration is built up, similar to the experiment (Rim et al., 2005). The configuration includes a cylindrical drug reservoir and the outer part of the human skin, namely the epidermis, which includes the stratum corneum and the viable epidermis. The reservoir contains a finite amount of fentanyl. The geometrical specifications and transport properties used in the simulations are indicated in Table 1.

TABLE 1

TABLE 1. Material transport properties are used in the model for the base case for fentanyl from (Rim et al., 2005).

Given the very similar radius of the patch and the skin sample in the experiment, the transport processes can be considered one-dimensional. This led to minimal differences when comparing the 3D and 1D model representation (Defraeye et al., 2020). Differences in predicted fluxes and the cumulative drug amount that was taken up by the aqueous buffer solution were <0.4%. Thus, a 1D model was sufficiently accurate and improved the computational economy, which was beneficial because a large parametric space was explored. Such a 1D model is representative of Franz diffusion cell experiments (Larsen et al., 2003; Rim et al., 2005; Bartosova and Bajgar, 2012). The skin is as wide as the drug reservoir, so unidirectional (longitudinal) drug transport is simulated without any transverse component.

A patch radius R_pt of 9 mm was used for the drug reservoir, similar to in the experiments. This value implies an active area of 2.54 cm², which is in the same order of magnitude as reported for commercial transdermal patches for fentanyl (4.2–42 cm²) (Wiedersberg and Guy, 2014). The thickness of the patch (d_pt) was 50.8 μm. The volume of the patch reservoir was 1.29 × 10^–8 m³. The skin’s epidermis (ep) was modeled as a single layer. This approach lumps the markedly different diffusion and partitioning processes through the lipophilic stratum corneum and the hydrophilic viable epidermis. This strategy does not affect the uptake kinetics for such a one-dimensional setup without extensive lateral transport. As such, a realistic simulated uptake drug flux across the skin into the Franz diffusion cell g_bl,up(t) [kg m⁻² s⁻¹] is obtained. It is the performance metric used to compare with experimental data. However, this one-layer approach does not accurately reflect the spatial concentration gradients within the skin, as shown in (Defraeye et al., 2020), where much higher concentrations are found within the stratum corneum.

Governing Equations

The governing equations for simulating transient transport of fentanyl are detailed. Only drug diffusion was solved, and isothermal conditions were assumed, close to human body temperature. Water transport due to skin de-/rehydration and the resulting skin shrinkage or swelling were not modeled. The following mass conservation equation was derived for fentanyl (substance α) in (Defraeye et al., 2020). This equation is defined for each material i [kg m⁻³] (patch and epidermis) to calculate the drug potential ψ^α:

$K_i^{\alpha }\frac{\partial {\psi }^{\alpha }}{\partial t}+\nabla \cdot \left(-D_i^{\alpha }{K}_i^{\alpha }\nabla {\psi }^{\alpha }\right)=S_s^{\alpha },(1)$

D_i^α is the diffusion coefficient or diffusivity [m² s⁻¹], K_i^α is the drug capacity of the drug in the material i [−], S_s^α is a volumetric source term for substance α [kg m⁻³ s⁻¹], and t is the time [s].

The conservation equation was derived in drug potential (ψ^α) instead of drug concentration (c_i^α = K_i^αψ^α). The use of ψ^α instead of c_i^α leads to a single dependent variable throughout the entire domain, instead of one in each material. This choice mainly avoids numerical stability issues at the interface of the patch and epidermis (Defraeye et al., 2020). Such instabilities can occur due to the discontinuity in the drug concentration that arises at the interface caused by partitioning.

Partitioning implies that when a drug substance α is brought into contact with patch and epidermis (pt and ep), the drug will evolve to a different equilibrium concentration in each of these materials, namely c_pt^α and c_ep^α. The ratio of these equilibrium concentrations is referred to as the partition coefficient K_pt/ep^α.

$K_{pt/ep}^{\alpha }=\frac{c_{pt}^{\alpha }}{c_{ep}^{\alpha }}=\frac{K_{pt}^{\alpha }{\psi }^{\alpha }}{K_{ep}^{\alpha }{\psi }^{\alpha }},(2)$

The octanol-water partition coefficient is frequently determined (K_o/w^α) for drug partitioning in liquids, where values larger than one indicate drug lipophilicity. Values smaller than one indicate drug hydrophilicity. The log (K_o/w^α) is often reported, where positive/negative values indicate lipophilicity/hydrophilicity, respectively. In this study, the partition coefficient between the patch and the epidermis was explicitly determined from experiments, as detailed in Material Properties and Transport Characteristics of Skin and Patch.

No source term was included in this study. The reason is that the contributions of the following processes (Naegel et al., 2013) could be neglected for fentanyl, as motivated in (Defraeye et al., 2020): 1) metabolization of the drug molecule within the epidermis by a chemical reaction, a process that leads to a conversion of the drug into other compounds; 2) adsorption of the drug molecule into the epidermis, and thus physical binding of the drug molecules.

Material Properties and Transport Characteristics of Skin and Patch

Fentanyl is a synthetic opioid that is used as a pain medication. It has a low molecular weight (337 Da) and is moderately lipophilic with a log (K_o/w) of three–four (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 1 for fentanyl for the base case, as taken from (Rim et al., 2005). These parameters were estimated in (Rim et al., 2005) as follows:

• The diffusion coefficient of fentanyl in the patch (D_pt^α) was determined experimentally. This was done by performing a drug release experiment of the patch into a receptor medium, after which an analytical expression was fitted to the results.

• The drug capacity of fentanyl in the epidermis (K_ep^α = 0.14) was derived from sorption experiments, measured over five skin donors. The patch’s drug capacity (K_pt^α) was taken equal to one since it is taken as the reference domain. Note that in (Rim et al., 2005), these drug capacities are actually termed partition coefficients, which is not the same definition used in the present study. The actual partition coefficient between patch and epidermis K_pt/ep^α = K_pt^α/K_ep^α equals 1/0.14.

• The value for the baseline diffusivity of the drug in the epidermis (D_ep^α) was determined by fitting the experimentally measured drug flux in the Franz diffusion cell over a range of initial drug concentrations in the patch for the initial 30 h (Rim et al., 2005).

These different procedures illustrate the complexity of obtaining material transport properties for transdermal drug delivery. Three different experimental setups were required, where even model fitting was required for one parameter (D_ep^α). In addition, each experimental technique had specific uncertainty.

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. This assumption implies that they are independent of the drug concentration. This choice is often unavoidable because more detailed data is rarely available. However, diffusion and partition coefficients have been shown to be a function of drug concentration rather than constant values (Lundborg et al., 2018).

Boundary and Initial Conditions

The drug was assumed to exit the computational domain via the interface with the aqueous solution in the Franz diffusion cell. Therefore, a constant concentration (and potential), equal to zero, was set at the epidermis bottom (Figure 1). This approximates the very low concentration found in the large buffer solution. This condition represents a Dirichlet boundary condition. Zero-flux conditions were imposed at all vertical boundaries. At t = 0 s, the skin was assumed to be drug-free. The initial concentration of drugs in the patch was set at 80 kg m⁻³ or 60 kg m^−3, according to a previous study (Rim et al., 2005). The patch dimensions combined with these initial concentrations lead to an initial amount of fentanyl in the reservoir (m_pt,ini) of 1.03 and 0.776 mg, respectively. This initial drug content corresponds to amounts typically present in commercially available patches (Defraeye et al., 2020). Perfect contact between the patch and the epidermis was assumed. No inclusion of air layers or discontinuities like hairs or skin roughness were accounted for.

Spatial and Temporal Discretization

The grid was built based on a grid sensitivity analysis. The spatial discretization error was quantified based on the total mass flux to the Franz diffusion cell. This error was estimated to be 0.1%, as determined by Richardson extrapolation (Roache, 1994; Franke et al., 2007). The grid consisted of 112 quadrilateral finite elements for the base case (1D, elements with a size of about 1 µm). The grid was gradually refined toward the patch-skin interface to enhance numerical accuracy and stability. The reason is that the largest gradients are found at such interfaces, particularly at the uptake process initial stage.

The transient simulations quantified a drug uptake process that lasted 72 h (3 days), starting from these initial conditions. This time window also agrees with that of conventional transdermal fentanyl therapy. Here the required dose for the patient, so patch size is estimated empirically by the clinician. Afterward, the patch is applied transdermally and is replaced every 72 h (Muijsers and Wagstaff, 2001).

The simulations applied adaptive time-stepping. The maximal time step was 600 s (10 min). This time step ensured high temporal resolution for the output data and was determined from a sensitivity analysis.

Numerical Implementation and Simulation

The model was implemented in COMSOL Multiphysics® software (version 5.5, COMSOL AB, Stockholm, Sweden). COMSOL is commercial finite-element-based software. This software was verified by the code developers. Therefore, additional code verification was not performed by the authors. Transient diffusive drug transport (Eq. 1) in the patch and epidermis was solved using the partial differential equations interface (coefficient form). The conservation equation was solved for the dependent variable ψ. The optimization interface was used to automatically find the optimal model-parameter values by direct inverse modeling (Halliday, 2020). Here, the least-squares method was used. This method minimizes the objective function, namely, the root mean square of the differences between experimental data and simulations. Different starting points were applied to avoid landing in a local minimum. In addition, the results were compared to the results of the full parametric space as well to ensure no local minimum was identified.

Quadratic Lagrange elements were used together with a fully coupled direct solver, which relied on the MUMPS solver scheme (MUltifrontal Massively Parallel sparse direct Solver). For the direct inverse modeling, a derivative-free solver was used, namely the Bound Optimization by Quadratic Approximation (BOBYQA). The tolerances for solver settings and convergence were determined through sensitivity analysis so that a further increase in the tolerances did not alter the resulting solution. For the inverse modeling, different starting values of the material properties were explored, and the sensitivity to the stopping criterion parameters was also performed.

Parametric Study

Following simulations were performed:

1. The base case simulates the drug uptake through the epidermis as released from a finite drug reservoir based on experimentally determined material properties (Table 1) (Rim et al., 2005). Simulations were done for c_pt,ini^α = 80 kg m⁻³.

2. Indirect multiparameter inverse modeling to determine the optimal set of material properties (D_pt^α, D_ep^α, K_ep^α) that gives the best agreement with the experimental data during a 72-h drug uptake period (Sensitivity of Drug Uptake Kinetics to Material Properties) for c_pt,ini^α = 80 kg m⁻³. Indirect modeling implies that a large set of combinations of these properties was explored; in this case, 30,618 combinations/simulations. The main aim here was to get insight into how the relationship between these parameters affects solution accuracy.

3. Direct multiparameter inverse modeling to determine the optimal set of material properties (D_pt^α, D_ep^α, K_ep^α) that gives the best agreement with the experimental data (Finding the Optimal Set of Material Properties by Multiparameter Inverse Modeling) for c_pt,ini^α = 80 kg m⁻³. Here automated least-squares optimization was used to progress fast to the most optimal solution. The main aim here is to identify how much faster a solution can be found compared to indirect inverse modeling. We also performed direct inverse modeling for c_pt,ini^α = 60 kg m⁻³.

4. Cross verification of the obtained parameter set by direct inverse modeling (for c_pt,ini^α = 80 kg m⁻³) by comparison with independent experimental data of patches with a different drug concentration (c_pt,ini^α = 60 kg m⁻³). This independent verification proves the generalized use of the model, so if the calibrated model can be used for a wider operational range of patches. We also performed the same cross-verification by evaluating a patch with an initial drug concentration of c_pt,ini^α = 80 kg m⁻³ using the optimal model parameters obtained from direct inverse modeling for c_pt,ini^α = 60 kg m⁻³.

Metrics to Evaluate TDD

The simulated drug delivery process was analyzed quantitatively by calculating several metrics. From the experiments, the flux taken up by the Franz diffusion cell (g_Fr,up [kg m⁻² s⁻¹]) was measured as a function of time for discrete points in time. This experiment was done by measuring the total amount of drugs taken up in the respective timeframe into the Franz diffusion cell through surface 1 [kg s⁻¹]. Afterward, this flow rate was scaled with the patch’s surface area (surface 2, R_pt = 9 mm). This flux quantity was used as the primary metric to compare with the 1D simulations, where it was determined in the same way. The Root Mean Square Deviation (RMSD) between measured fluxes and the predicted fluxes by the simulations over the N measured points at each time j was calculated as:

$\mathrm{RMSD}=\sqrt{\frac{{\displaystyle \sum _j{\left(g_{Fr,up}^{\mathrm{SIM}}\left(t_j\right)-g_{Fr,up}^{\mathrm{EXP}}\left(t_j\right)\right)}^2}}{N}},(3)$

This RMSD [kg m⁻² s⁻¹] is the quadratic mean over these differences in uptake fluxes between experiments and simulations. Thereby it is a direct measure of the accuracy of the simulations. Using the RMSD, we can quickly elucidate how far off these parameters are from the experimental data and the original dataset (base case). The advantage of the RMSD is that it aggregates each simulation’s accuracy into one single value, being the prediction error of the simulation. An additional benefit of the RMSD for comparison is that it has the same units as the flux [kg m⁻² s⁻¹], so it can be directly contrasted to the measured fluxes. The RMSD was calculated considering the entire 72 h timeframe (N = 9 data points) and a shorter 30 h timeframe (RMSD′, N = 4).

In addition to the RMSD, the average uptake flux across the skin into the diffusion cell over the 72 h is quantified ${\overline{g}}_{Fr,up}^{\mathrm{SIM}}$ [kg m⁻² s⁻¹]. The reason is that the normalized RMSD (NRMSD) can be calculated as the ratio of the RMSD to this average flux, which can be expressed as a percentage as well:

$\mathrm{NRMSD}=\frac{\mathrm{RMSD}}{{\overline{g}}_{Fr,up}^{\mathrm{SIM}}},(4)$

Results and Discussion

Sensitivity of Drug Uptake Kinetics to Material Properties

We aim to quantify how the model parameters—the material properties D_pt^α, D_ep^α, and K_ep^α—affect the simulations’accuracy, so the RMSD. We mainly aim to elucidate if there are particular combinations of model parameters that enhance model accuracy. These results shed light on which transport kinetics are decisive for achieving an accurate mechanistic model. Different subsets of the material properties were analyzed within the sizable parametric space explored for D_pt^α, D_ep^α, and K_ep^α. First, we analyze the impact of epidermal diffusion vs. partitioning in the epidermis. Second, the effect of diffusion in the epidermis vs. diffusion in the patch is targeted. Finally, the most accurate combination of parameters in this design space that leads to the lowest RMSD is identified.

The Impact of Epidermal Diffusion and Partitioning

Parameter variations in the diffusion (D_ep^α) and partition coefficient of the epidermis (K_ep^α) are explored. The diffusion coefficient in the patch (D_pt^α) is kept the same as the base case. The resulting RMS

D values for each combination of D_ep^α and K_ep^α are shown in Figure 2. The drug uptake profiles for the base case (c_pt,ini^α = 80 kg m⁻³) and the combination of D_ep^α and K_ep^α out of this parametric space that lead to the lowest RMSD are depicted in Figure 3. We evaluated only a discrete number of combinations of D_ep^α and K_ep^α. As such, a more optimal combination of D_ep^α and K_ep^α with a lower RMSD could be found when refining even more. This step was done in Finding the Optimal Set of Material Properties by Multiparameter Inverse Modeling by direct inverse modeling. However, the aim of the current section was to analyze the processes instead of finding the perfect set of parameters.

FIGURE 2

FIGURE 2. RMSD over the uptake period as a function of D_ep^α and K_ep^α (for D_pt^α of the base case = 1.00 × 10^–13 m² s⁻¹) for an initial patch concentration of 80 kg m⁻³. The parameters for the base case and most optimal parameter set are also shown.

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 an initial patch concentration of 80 kg m⁻³ as a function of time. Results are shown of the base case, the best combination of epidermal parameters (Sim.—epidermis, D_ep^α and K_ep^α, for D_pt^α of the base case = 1.00 × 10^–13 m² s⁻¹), and the best combination of D_ep^α, K_ep^α, and D_pt^α (Sim.—epidermis and patch) from indirect multiparameter inverse modeling.

The overall accuracy of the simulations is quantified with the RMSD (Figure 2). The RMSD of the base case (0.40) already lies quite close to the minimal RMSD value that is achieved for a more optimal combination of epidermal partition and diffusion coefficients (0.33). However, a significant improvement in RMSD was still made. The normalized RMSD was 22%, where the base case yielded 26%. In this case, the corresponding D_ep^α and K_ep^α were 4.8 × 10^–14 m² s⁻¹ and 0.067, respectively, in contrast to 3.00 × 10^–14 m² s⁻¹ and 0.14 for the base case. Interesting observations can be made when analyzing the RMSD landscape (Figure 2). A “valley/canyon” with very low values of the RMSD is present. This valley stretches over a wide range of D_ep^α and K_ep^α. This contrasts with a typical optimization problem where an exact global minimum is present in the multiparameter space. As such, there is a broad range of possible combinations of diffusion and partition coefficients that lead to a low RMSD, so a relatively accurate simulated drug uptake process. This is problematic as an unphysical value of D_ep^α, for example, similar to that of the patch, can give an acceptable agreement with the experimental data when tuning/calibrating the model with the appropriate K_ep^α.

The simulation accuracy throughout the drug uptake process is assessed by comparing the predicted drug uptake flux over time (Figure 3). The simulations capture the initial period of the uptake accurately, namely the substantial increase in the flux. This increase occurs when drugs are transported into the epidermis. This leads to an increased drug concentration in the epidermis, so a part of the drugs is stored there. After this initial phase, drugs do not accumulate anymore in the epidermis, and a quasi-steady-state condition sets in where the amount of drugs entering the epidermis also diffuses out (Figure 4). Because the skin’s capacity to store drugs is reached, the drug primarily diffuses through the epidermis, and the stored amount remains relatively constant. Since the drug concentration in the patch (a finite reservoir) is decreasing over time, steady-state conditions with constant flux are not reached. Instead, the uptake flux slowly decreases over time, simply because the concentration gradient over the epidermis decreases. The experiments predict a much steeper decline in the flux than in the simulations, in which the predicted decline has a relatively constant slope, so is linear. During this period, the simulations do not lie within the experimental error bars. The simulations seem to miss a critical physical process, irrespective of the combination of D_ep^α and K_ep^α that is used. An answer is sought and found in the next section.

FIGURE 4

FIGURE 4. Drug flux released by the patch into the epidermis and taken up by blood as a function of time for the base case simulation and the optimal parameter combination for inverse modeling. The diffusion in the reservoir vs. in the epidermis.

The Diffusion in the Reservoir vs. in the Epidermis

We aim to unveil whether the diffusion kinetics in the patch could help explain and improve the residual differences found between experiments and simulations in The Impact of Epidermal Diffusion and Partitioning. To this end, parameter variations in the diffusion coefficient of the patch D_pt^α are included, in addition to those in D_ep^α and K_ep^α. A subset of the resulting RMSD values is given in Figure 5, as depicted for different K_ep^α. These contour plots show how changing the ratio of the resistance to drug diffusion in the patch (D_pt^α), relative to that in the epidermis (D_ep^α), affects the drug uptake. The drug uptake profiles for the base case (c_pt,ini^α = 80 kg m⁻³) and for the combination of D_ep^α, D_pt^α, and K_ep^α out of this parametric space that lead to the lowest RMSD, are depicted in Figure 3.

FIGURE 5

FIGURE 5. RMSD over the uptake period as a function of D_ep^α and D_pt^α for multiple values of K_ep^α for an initial patch concentration of 80 kg m⁻³. The parameters for the base case and the most optimal parameter set are also shown.

The overall accuracy of the simulations is evaluated using the RMSD (Figure 5). A more optimal combination of D_ep^α, D_pt^α, and K_ep^α made a significant improvement in RMSD (0.141) compared to only varying the epidermal model parameters (RMSD = 0.33, The Impact of Epidermal Diffusion and Partitioning) or the base case (RMSD = 0.4). The normalized RMSD was 9.3%. The base case yielded 26%, and the combination of D_ep^α and K_ep^α produced 22%. The corresponding D_ep^α, D_pt^α, and K_ep^α were in this case 3.0 × 10^–14, 6.91 × 10^–16, and 2.15, respectively, in contrast to 3.00 × 10^–14, 1 × 10^–13, and 0.14 for the base case. Here, “valleys/canyons” with low values of the RMSD are present in all graphs. A broad range of possible combinations of diffusion coefficients (D_ep^α and D_pt^α) lead to a simulated drug uptake process with a low RMSD. This is true for each of the selected partition coefficients since the shape of the surface/contour plots are similar. This finding also means that partitioning between patch and epidermis does not significantly affect the relative importance of the diffusive processes in these two domains. This is confirmed in Figure 6, where the minimal values of the RMSD are shown over the full parametric space for each value of D_ep^α, D_pt^α, and K_ep^α. Both D_ep^α and D_pt^α have an explicit minimum. However, for K_ep^α, there is an entire range where a RMSD below 0.15 is reached. Both values above and below one lead to a low RMSD. When the drug capacity of the epidermis is higher than one, namely the capacity of the patch (K_pt^α = 1), the partition coefficient (K_pt/ep^α = K_pt^α/K_ep^α = c_pt^α/c_ep^α) becomes smaller than one. This indicates that the patch’s drug concentration under equilibrium conditions is lower than the drug concentration in the epidermis. The drug will be preferably in the epidermis than in the patch, even in the absence of a concentration gradient.

FIGURE 6

FIGURE 6. Minimal values of the RMSD over the full parametric space for each value of (A) K_ep^α, (B) D_ep^α, (C) D_pt^α. The values of the base case are also shown.

The simulation accuracy throughout the drug uptake process is assessed by comparing the predicted drug uptake flux over time (Figure 3). The best combination of D_ep^α, D_pt^α, and K_ep^α leads to an accurate drug uptake prediction throughout the entire 72 h timeframe, namely almost within the experimental error bars. Both the steep increase in the flux during the initial uptake period and the subsequent nonlinear decline are captured. Based on these results, we can elucidate the physical reasons for the excellent agreement for the specific set of D_ep^α, D_pt^α, and K_ep^α. The lowest RMSD and the best agreement of the uptake flux profiles were obtained for a lower D_pt^α, compared to the base case. D_pt^α was determined initially experimentally. Its magnitude resulted in much easier diffusive drug transport in the patch than in the epidermis (1.00 × 10^–13 m² s⁻¹ vs. 3.00 × 10^–14 m² s⁻¹). The primary resistance to diffusive transfer from the patch through the epidermis into the diffusion cell was located in the epidermis for the base case. This diffusive resistances R_i [s m⁻¹] is defined as:

$R_i=\frac{d_i}{D_i^{\alpha }},(5)$

where d_i is the thickness of the patch or the epidermis [m]. The corresponding resistances for patch and epidermis were 5.08 × 10⁸ and 1.69 × 109 s m⁻¹, respectively, for the base case. A better agreement with experiments is obtained when the diffusion in the reservoir holding the drug is restricted more so that it also has a notable resistance to drug transport. As a result of the reduced transport in the patch, even a non-uniform concentration in the patch will arise, in contrast to the base case, as depicted in Figure 7. Apart from the decreasing concentration in the patch, which leads to a reduced gradient over the skin, the concentration gradient inside the patch also plays a role. As a result, both the steep increase in the drug uptake flux and the nonlinear decline are predicted correctly (Figure 4).

FIGURE 7

FIGURE 7. Color contours of drug concentration in the drug reservoir and epidermis for different points in time for (A) the base case and, (B) simulations with optimal material properties as determined from inverse modeling.

This lower diffusion coefficient of the patch is suggested by the physics-based simulations, as this is required to obtain the measured drug uptake process accurately, namely the steep rise in drug flux and its subsequent sharp nonlinear decline. With the experimentally determined material properties of the patch, the realistic drug uptake behavior could not be accurately captured by the model (Figure 3). The question remains why there is a mismatch with the measured diffusion coefficient of the patch, with which the physical drug uptake behavior could not be captured. This diffusion coefficient was however measured by a drug release experiment of the patch into an aqueous receptor medium. This case exposes the patch to different conditions as when placed on the skin and might induce a higher diffusion rate or swelling of the patch. Therefore, it seems that the measured diffusion coefficient might be too high. Our simulations confirm that such a high diffusion coefficient does not allow to reproduce the experimental drug-uptake behavior. Therefore, a cross-verification of experimentally determined diffusion coefficients of transdermal patches by using techniques other then placement in an aqueous receptor medium could be advised.

Finding the Optimal Set of Material Properties by Multiparameter Inverse Modeling

We determine the optimal set of material properties (D_pt^α, D_ep^α, K_ep^α) by indirect and direct multiparameter inverse modeling. For indirect modeling, we searched the combination in a sizable parametric space of 30,618 combinations that gives the best agreement with the experimental data. The parameter set with the best agreement was identified in The Diffusion in the Reservoir vs. in the Epidermis. For direct inverse modeling, we converged in a more straightforward way to the optimal solution, without exploring the full parametric space, so with a much lower amount of simulations. The optimal combinations of material properties can be found in Table 2, together with the RMSD for each case. The drug uptake flux profile for direct and indirect inverse modeling and the RMSD at each time step are shown in Figure 8.

TABLE 2

TABLE 2. Material transport properties for the base case simulations from experiments (Rim et al., 2005) and for the optimal solution for indirect and direct inverse modeling and for cross verification simulations, using the optimal parameters to simulate drug uptake with a different initial concentration.

FIGURE 8

FIGURE 8. (A) Drug fluxes of fentanyl at surface 1 (g_bl,up), so leaving the epidermis, from experiments (Exp.) and simulations (Sim.) for an initial patch concentration of 80 kg m⁻³ as a function of time. Results are shown of the base case and the optimal combination of D_ep^α, K_ep^α, and D_pt^α by direct and indirect inverse modeling. (B) Corresponding normalized RMSD at each experimental time step.

Indirect inverse modeling gives a much better solution than with the experimentally determined values. We improved the normalized RMSD from 26 to 9.3% and reduced the maximal deviation (per point) from 41 to 19%. Here this maximal deviation was calculated as the local NRMSD for each point (Eqs 3, 4), N = 1). However, the sizable parametric space is systematically explored, which is time-consuming, and only at discrete combinations of parameters D_pt^α, D_ep^α, K_ep^α. As only a discrete amount of combinations was analyzed, indirect inverse modeling gives the best (discrete) solution in this space but not the most optimal solution in the continuous space. This can be improved by running another parametric study around the optimal set of parameters at a higher resolution, so smaller steps between the parameter values. Nevertheless, the number of simulations to obtain this solution this way is vast.

Direct inverse modeling gives an even better RMSD than indirect modeling and much faster as well. We improved the normalized RMSD from 26 to 8.8% and reduced the maximal normalized RMSD deviation from 41 to 16%. Only 306 simulations were needed to attain the optimal solution, which is only 1% of what we had with indirect inverse modeling (30,618). Note, however, that additional simulations were performed to check the solution’s sensitivity to the starting point, which adds other computations.

Note that we also performed direct inverse modeling for c_pt,ini^α = 60 kg m⁻³. The optimal combinations of material properties are also depicted in Table 2, together with the RMSD. We improved the normalized RMSD from 24 to 2.4% and reduced the maximal normalized RMSD deviation from 26 to 3.5%.

Cross Verification of Optimized Model

We verify if the optimal set of transport parameters that were derived for the epidermis and patch in Finding the Optimal Set of Material Properties by Multiparameter Inverse Modeling for a specific initial drug concentration (c_pt,ini^α = 80 kg m⁻³), can be used for a more comprehensive operational range of patches. To this end, we simulate drug uptake with another initial concentration (c_pt,ini^α = 60 kg m⁻³), for which also experimental data is available, but using the optimal parameter dataset from direct inverse modeling (Table 2), determined for (c_pt,ini^α = 80 kg m⁻³). Vice versa, we also performed the same cross verification for an initial drug concentration of c_pt,ini^α = 80 kg m⁻³ using the optimal model parameters obtained from direct inverse modeling for c_pt,ini^α = 60 kg m⁻³. The drug uptake flux profiles for this simulation and the experimental data are shown in Figure 9, and the RMSD is shown in Table 2. The results show that there is a good agreement for both concentrations, which is in both cases better than the simulations where experimentally determined material properties were used. For c_pt,ini^α = 60 kg m⁻³ (Figure 9A), the concentration profile lies entirely within the experimental uncertainty. This means that a model, calibrated via inverse modeling on a specific concentration, can successfully be used to simulate patches with other doses as well. As such, the calibration, based on inverse modeling, only needs to be done once, and the calibrated model is valid for a more comprehensive operating range. Note that in typical physics-based modeling, normally the transport properties are determined experimentally for a single set of conditions, or fitted to match a certain set of experiments (Rim et al., 2005). The variations of the material properties over a wide range of operating conditions are rarely explored when using physics-based modeling. Diffusion or partition coefficients that are a function of drug concentration are rarely determined (Lundborg et al., 2018), and constant values are used. Here we showed that this variation, over a limited range of operating conditions, is already rather limited.

FIGURE 9

FIGURE 9. Drug fluxes of fentanyl at surface 1 (g_bl,up), so leaving the epidermis, from experiments (Exp.) for an initial patch concentration of 60 kg m⁻³(A) and 80 kg m⁻³(B) as a function of time and simulations (Sim.). (A) Simulation results are shown of the base case with an initial patch concentration of 60 kg m⁻³ and the optimal combination of D_ep^α, K_ep^α, and D_pt^α, but where the parameters are derived for an initial patch concentration of 80 kg m⁻³. (B) Simulation results are shown of the base case with an initial patch concentration of 80 kg m⁻³ and the optimal combination of D_ep^α, K_ep^α, and D_pt^α, but where the parameters are derived for an initial patch concentration of 60 kg m⁻³.

Discussion and Outlook

A good agreement with experimental data was obtained in the present study for the optimal set of material property data, as determined by inverse modeling. This indicates that the primary physical processes at play are included in, and accurately captured by, the mechanistic model, namely diffusion and partitioning within the epidermis and transdermal patch. However, it is possible that secondary processes affect drug uptake as well, which are now not included yet. The processes that could be additionally modeled are 1) swelling or shrinkage of the epidermis caused by changing hydration due to the presence of the patch, since the skin can swell significantly; 2) chemical metabolization (or binding) of the fentanyl molecule in the epidermis; 3) physical binding of the drug molecules so adsorption of the drug molecule into the epidermis; 4) differential diffusion in the epidermis, which can be accounted for by splitting it up into a stratum corneum and a viable epidermis, since the majority of the resistance to diffusion and the primary drug storage occurs in the SC; 5) diffusion and partition coefficients that are a function of drug concentration, instead of constant values. The importance of these secondary effects should be assessed on a case by case basis. For fentanyl, the high bioavailability through the skin (e.g., 92% in (ASHP, 2017)) implies, for example, that physical adsorption or chemical metabolization would not play a decisive role.

In addition, temperature can also affect the permeability of the skin for the drugs, and induce an increase of the drug penetration (US Food and Drug Administration, 2005; Iikura et al., 2019). Using a pharmacokinetic model, serum fentanyl concentrations were shown to increase by approximately one-third for patients with a body temperature of 40°C due to temperature-dependent increases in fentanyl released from the system and the increased skin permeability. Also, it has been shown that higher temperatures lead to a higher permeation coefficient of the skin for hydrophilic drugs. Such studies in fentanyl transdermal therapy showed that an increase in temperature might lead to an overdose for the patient. For transdermal therapy where an elevated temperature is at play, the physics-based model needs to account for this to ensure reliability. This can be done by upgrading the model, based on the experimental data, to account for the effect of temperature on the permeation of drugs through the skin. As such, the mechanistic model is able to predict the amount of absorbed drugs at elevated temperatures and the associated risk. The inverse modeling approach can thus also be used to quantify the optimal skin penetration parameters at different temperatures.

The inverse modeling presented could significantly reduce experimental work and resources (typically in vitro) that are associated with mechanistic modeling since it is no longer required to determine each of the model parameters separately. This inverse modeling only needed a single experiment on the resulting process—drug uptake via a Franz diffusion cell experiment. An experiment at another patch drug concentration was then used for verification. Thereby, the experimental determination of the diffusion coefficients of the epidermis and patch and the partition coefficient was avoided (D_pt^α, D_ep^α, and K_ep^α), which would imply three additional experiments on different equipment. The idea behind such inverse modeling to determine the material properties in-silico is that these properties are derived on a “training dataset” but have a sufficient accuracy to be generally applied afterward for simulations of similar cases. Verification is then used for at least one independent dataset. This can save time and resources, hence speeding up the simulation process.

The potential value of such in-silico determination of material properties becomes even larger if one wants to use mechanistic modeling for therapy for individual patients. Such mechanistic models are a fundamental building block of so-called digital twins (Defraeye et al., 2020). These twins could be used for steering next-generation transdermal drug delivery systems. In clinics, significant additional resources would be required to do experiments on skin samples to calibrate the model for a specific patient. Furthermore, such lab experiments measure material properties on separate parts of the modeled system (e.g., epidermis) and different conditions (e.g., aqueous solution), which differ from reality, living human tissue. Such in-silico skin models can be used to fit the optimal model parameters, just by measuring the response of the patient to a specific biomarker that is representative of the drug of interest. Once calibrated for a particular patient, the in-silico skin model can be used to assess/design therapy within a specific operating range. This concept however still requires an implementation plan for the clinics as well as a rigorous clinical verification.

Indirect inverse modeling provided a clear insight into the impact and sensitivity of different model parameters and elucidated the optimal parameter set in-silico. Running 30,618 simulations was possible in the present study as a 1D model was used, but it will be computationally very demanding for more complex models. Direct inverse modeling is the fastest way to reveal the optimal set of parameters were and required 100 times fewer simulations to be performed. On the other hand, direct inverse modeling does not enable us to analyze and interpret the shape of the parametric space in detail, so the role of the different physics.

Finally, the physics-based modeling is set up for small drug molecules that are delivered transdermally via first-generation systems (Prausnitz and Langer, 2008). Apart from fentanyl, the modeling approach can be used already for other lipophilic drugs with a small molecular weight. Examples are buprenorphine, rotigotine, rivastigmine, or ibuprofen. For the delivery of larger molecules like insulin, more advanced next-generation drug delivery systems are required (Bartosova and Bajgar, 2012; Lee et al., 2018). These systems disrupt the stratum corneum using microneedles or thermal ablation increasing the drug uptake via chemical permeation enhancers or iontophoresis (Bartosova and Bajgar, 2012).

Conclusion

We evaluated multiparameter inverse modeling as a fast and reliable alternative to determine the drug diffusion and partition coefficients for transdermal delivery of fentanyl using mechanistic modeling. We found that capturing the drug uptake process accurately by simulations was only possible if the diffusion in the reservoir holding the drug was restricted to some extent, via D_pt^α, in contrast to what the experimentally measured values suggested. This implied that a non-uniform concentration in the patch arose with distinct gradients. Only then, the initial steep increase in the drug uptake flux and the later nonlinear decline were predicted correctly.

Indirect inverse modeling was used to determine the skin and patch’s optimal material properties, namely the diffusion, and partitioning coefficients. A much better agreement with the experimental drug uptake was obtained, compared to when using the experimentally determined material properties. We improved the normalized RMSD from 26 to 9.3% and reduced its maximal deviation (per point) from 41 to 19%. A large range of appropriate combinations of diffusion and partition coefficients was identified that led to an accurately simulated drug uptake process. Hence, in this optimization problem, there was no clear minimum but rather a valley of optimal values. Direct inverse modeling gave an even better RMSD than indirect inverse modeling. We improved the normalized RMSD from 26 to 8.8% and reduced its maximal deviation from 41 to 16%. In addition, only 306 simulations were needed to attain the optimal solution compared to 30,618 with indirect inverse modeling, so a speedup of factor 100.

By cross verification, we showed that our model, when calibrated via inverse modeling for a specific concentration, can successfully be used to simulate patches with other drug doses as well. As such, the calibration, based on inverse modeling, only needs to be done once to have a valid model for a more comprehensive operating range. We thereby successfully demonstrated a fast way to build up an accurate mechanistic model for drug delivery. This opens the door to clinically ready patient-specific therapy tailored to patients.

Data Availability Statement

The raw data supporting the conclusion 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 FB; TD supervised FB; TD wrote the original draft of the paper and did the visualization, with key input from FB; RR performed critical review and editing.

Funding

This work was supported by the Novartis Research Foundation (grant “Virtual twinning for intelligent, personalized transdermal drug delivery”). We acknowledge the support of Chandrima Shrivastava in the first exploratory work on data processing of the simulation results. 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.

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.

Nomenclature

Symbols

A_pt Active area of the patch [m²]

c_i^α Drug concentration of substance α in material i [kg m⁻³]

c_pt,ini^α Initial concentration in the patch [kg m⁻³]

d_ep Thickness of epidermis [m]

d_pt Thickness of patch [m]

D_i^α Diffusion coefficient/diffusivity of substance α in material i [m² s⁻¹]

g_bl,up(t) Uptake flux across the skin into the blood at a specific point in time [kg m⁻² s⁻¹]

g_Fr,up(t) Uptake flux across the skin into the Franz diffusion cell at a specific point in time [kg m⁻² s⁻¹]

K_A/B^α Partition coefficient between material A and B for substance α

K_o/w^α Partition coefficient between octanol and water for substance α

K_i^α Drug capacity of substance α in material i [−]

N Number of measured points

R_i Diffusive resistance of a material [s m⁻¹]

R_ep Radius of the epidermis in the model [m]

R_ep Radius of the transdermal patch [m]

S_s^α Volumetric source term for substance α [kg m⁻³ s⁻¹]

t Time [s]

Greek symbols

α substance indicator

ψ^α drug potential of substance α [kg m⁻³]

Subscripts

bl Blood

i Material indicator

ini Initial

ep Epidermis

up Uptake

pt Patch

Sim. In simulations

Exp. In experiments

Abbreviations

TDDS Transdermal drug delivery systems

RMSD Root mean square deviation

NRMSD Normalized root mean square deviation

References

Akkaram, S., Beeson, D., Agarwal, H., and Wiggs, G. (2007). Inverse modeling technology for parameter estimation. Struct. Multidisc Optim. 34 (2), 151–164. doi:10.1007/s00158-006-0067-1

CrossRef Full Text | Google Scholar

Amarah, A. A., Petlin, D. G., Grice, J. E., Hadgraft, J., Roberts, M. S., and 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

ASHP(2017). American society of health-system pharmacists: drug information. Bethesda, MD: ASHP.

Barbero, A. M., and Frasch, H. F. (2017). Effect of stratum corneum heterogeneity, anisotropy, asymmetry and follicular pathway on transdermal penetration. J. Control. Release 260, 234–246. doi:10.1016/j.jconrel.2017.05.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Barratt, M. D. (1995). Quantitative structure-activity relationships for skin permeability. Toxicol. In. Vitro 9 (1), 27–37. doi:10.1016/0887-2333(94)00190-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartosova, L., and 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

Caccavo, D., Barba, A. A., d'Amore, M., De Piano, R., Lamberti, G., Rossi, A., et al. (2017). Modeling the modified drug release from curved shape drug delivery systems—dome Matrix®. Eur. J. Pharm. Biopharm. 121, 24–31. doi:10.1016/j.ejpb.2017.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Casey, M., and Wintergerste, T. (2000). Special interest group on “quality and trust in industrial CFD” best practice guidelines. 1st Edn. Cambridge, United Kingdom; ERCOFTAC.

Chakravarty, K., and Dalal, D. C. (2019). A nonlinear mathematical model of drug delivery from polymeric matrix. Bull. Math. Biol. 81 (1), 105–130. doi:10.1007/s11538-018-0519-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chavoshi, S., Rabiee, M., Rafizadeh, M., Rabiee, N., Shamsabadi, A. S., Bagherzadeh, M., et al. (2019). Mathematical modeling of drug release from biodegradable polymeric microneedles. BDM 2 (2), 96–107. doi:10.1007/s42242-019-00041-y

CrossRef Full Text | Google Scholar

Chen, L., Han, L., Saib, O., and Lian, G. (2015). In silico prediction of percutaneous absorption and disposition kinetics of chemicals. Pharm. Res. 32, 1779–1793. doi:10.1007/s11095-014-1575-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Defraeye, T., Bahrami, F., Ding, L., Malini, R. I., Terrier, A., and Rossi, R. M. (2020). Predicting transdermal fentanyl delivery using mechanistic simulations for tailored therapy. Front. Pharmacol. 11, 585393. doi:10.3389/fphar.2020.585393

PubMed Abstract | CrossRef Full Text | Google Scholar

FDA(2016). Reporting of computational modeling studies in medical device submissions—guidance for industry and food and drug administration staff. Silver Spring, MD: Food and Drug Administration.

Ferreira, J. A., de Oliveira, P., and Pena, G. (2017). Transdermal iontophoresis—a quantitative and qualitative study. Comput. Math. Appl. 74 (10), 2231–2242. doi:10.1016/j.camwa.2017.07.001

CrossRef Full Text | Google Scholar

Filipovic, N., Saveljic, I., Rac, V., Graells, B. O., and Bijelic, G. (2017). Computational and experimental model of transdermal iontophorethic drug delivery system. Int. J. Pharm. 533 (2), 383–388. doi:10.1016/j.ijpharm.2017.05.066

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Gajula, K., Gupta, R., Sridhar, D. B., and Rai, B. (2017). In-silico skin model: a multiscale simulation study of drug transport. J. Chem. Inf. Model. 57 (8), 2027–2034. doi:10.1021/acs.jcim.7b00224

PubMed Abstract | CrossRef Full Text | Google Scholar

Halliday, A. (2020). COMSOL blog: how to use the parameter estimation study step for inverse modeling. Available at: https://www.comsol.com/blogs/how-to-use-the-parameter-estimation-study-step-for-inverse-modeling/ (Accessed July 09, 2020).

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 (2), 352–367. doi:10.1016/j.ejpb.2007.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Iikura, H., Uchida, K., Ogawa-Fuse, C., Bito, K., Naitou, S., Hosokawa, S., et al. (2019). Effects of temperature and humidity on the skin permeation of hydrophilic and hydrophobic drugs. AAPS PharmSciTech 20 (7), 1–9. doi:10.1208/s12249-019-1481-1

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 (44), D1202–D1213.

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, S., Bae, H., Jo, J., and Yoon, S. (2019). Comprehensive ensemble in QSAR prediction for drug discovery. BMC Bioinformatics 20 (1), 521. doi:10.1186/s12859-019-3135-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, R. H., Nielsen, F., Sørensen, J. A., and Nielsen, J. B. (2003). Dermal penetration of fentanyl: inter- and intraindividual variations. Pharmacol. Toxicol. 93 (5), 244–248. doi:10.1046/j.1600-0773.2003.pto930508.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H., Song, C., Baik, S., Kim, D., Hyeon, T., and 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

Lundborg, M., Wennberg, C. L., Narangifard, A., Lindahl, E., and 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

Madhihah, S., Malik, A., Abdullah, I., and Mahali, S. M. (2018). Analytic solution for hollow microneedles assisted transdermal drug delivery model. Int. J. Appl. Eng. Res. 13 (1), 737–742.

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

Muijsers, R. B., and Wagstaff, A. J. (2001). Transdermal fentanyl: an updated review of its pharmacological properties and therapeutic efficacy in chronic cancer pain control. Drugs 61 (15), 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., and 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 (2), 332–338. doi:10.1016/j.ejpb.2008.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Polat, B. E., Deen, W. M., Langer, R., and Blankschtein, D. (2012). A physical mechanism to explain the delivery of chemical penetration enhancers into skin during transdermal sonophoresis—insight into the observed synergism. J. Control Release 158 (2), 250–260. doi:10.1016/j.jconrel.2011.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Pontrelli, G., and De Monte, F. (2014). 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., and Langer, R. (2008). Transdermal drug delivery. Nat. Biotechnol. 26 (11), 1261–1268. doi:10.1038/nbt.1504

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rim, J. E., Pinsky, P. M., and Van Osdol, W. W. (2009). Multiscale modeling framework of transdermal drug delivery. Ann. Biomed. Eng. 37 (6), 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 (3), 405–413. doi:10.1115/1.2910291

CrossRef Full Text | Google Scholar

Ronnander, P., Simon, L., Spilgies, H., and Koch, A. (2018). Modelling the in-vitro dissolution and release of sumatriptan succinate from polyvinylpyrrolidone-based microneedles. Eur. J. Pharm. Sci. 125, 54–63. doi:10.1016/j.ejps.2018.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Selzer, D., Hahn, T., Naegel, A., Heisig, M., Kostka, K. H., Lehr, C. M., et al. (2013). Finite dose skin mass balance including the lateral part: comparison between experiment, pharmacokinetic modeling and diffusion models. J. Control. Release 165 (2), 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., and Schaefer, U. F. (2015). 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

Shirazi, R. N., Islam, S., Weafer, F. M., Whyte, W., Varela, C. E., Villanyi, A., et al. (1900). Multiscale experimental and computational modeling approaches to characterize therapy delivery to the heart from an implantable epicardial biomaterial reservoir. Adv. Healthc. Mater. 16, 1900228. doi:10.1002/adhm.201900228

CrossRef Full Text | Google Scholar

Toropova, A. P., and Toropov, A. A. (2017). The index of ideality of correlation: a criterion of predictability of QSAR models for skin permeability? Sci. Total Environ. 586, 466–472. doi:10.1016/j.scitotenv.2017.01.198

PubMed Abstract | CrossRef Full Text | Google Scholar

U. S. Food and Drug Administration(2005). Duragesic label. Silver Spring, MD: Food and Drug Administration.

U. S. Food and Drug Administration(2020). Adaptive designs for clinical trials of drugs and biologics: guidance for Industry. Silver Spring, MD: Food and Drug Administration.

Wiedersberg, S., and 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

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

CrossRef Full Text | Google Scholar

Zhan, X., Mao, Z., Chen, J., and Zhang, Y. (2015). Acrylate copolymer: a rate-controlling membrane in the transdermal drug delivery system. E-Polymers 15 (1), 55–63. doi:10.1515/epoly-2014-0123

CrossRef Full Text | Google Scholar