Modelling Transdermal Permeation of Volatiles from Complex Product Formulations

Abstract

Background: The evaporation of volatile ingredients from topical formulations strongly influences transdermal permeation and overall bioavailability, yet coupled evaporation-permeation dynamics are mostly simplified or neglected in existing models. Methods: We developed a mechanistic framework that couples Fickian gas-phase evaporation and transdermal permeation, both driven by the activity coefficients of volatiles. The model equations are implemented in a hybrid MATLAB-Python architecture with the volatile activity computed on-the-fly using UNIFAC and the gas-phase diffusivity calculated by the kinetic equation of Fuller-Schettler-Giddings (FSG). Initial validation used published IVPT data for 4-Tolunitrile and Nitrobenzene. Results: For 4-Tolunitrile, the FSG-based model estimated an initial evaporation coefficient of K_evap,i = 7.9348 × 10^-10 mol·cm^-2·s^-1, and parameter optimization converged to 8.3929 × 10^-11 mol·cm^-2·s^-1 (≈1/10 of the FSG estimate). The optimized model predicted an accumulation amount of 19.15% versus an experimental value of 16.97% in the receptor fluid (RF) at 24 h. For Nitrobenzene, the FSG initial estimation value of K_evap,i = 6.6480 × 10^-10 mol·cm^-2·s^-1 was optimized to 8.1174 × 10^-11 mol·cm^-2·s^-1 (≈1/8 of the FSG value), and the predicted amount of 24 h RF is 27.61% (experimental 23.19%). Both optimized K_evap,i values are roughly one order of magnitude lower than the initial FSG estimates, but >20× larger than Stokes-Einstein (SE)-derived values. Sensitivity scans show that further tuning of internal skin parameters (e.g., diffusion coefficient (D_SC,i) and partition coefficient (P_SCw,i)) produced only marginal improvements in RF prediction once K_evap,i was optimized. Conclusions: The coupled evaporation-permeation framework reproduces key IVPT kinetics for volatile solutes when the effective evaporation coefficient is calibrated. The kinetic-theory estimates (FSG-based) are a reasonable starting point, but typically overestimate the evaporation rate constant under finite-dose unoccluded IVPT conditions. By implementing the on-the-fly computation of volatile activity using UNIFAC, the approach is extensible to modelling transdermal permeation of volatiles from multicomponent/non-ideal formulations.

Keywords: diffusion; evaporation; permeation; stratum corneum; transdermal delivery.

Figure 1

Diagram of the simulation model, which follows the setting of the reference paper.

Figure 2

Solute distributions under initial calculated value of 4-Tolunitrile: K_evap,i = 7.9348 × 10⁻¹⁰ mol·cm⁻²·s⁻¹, D_SC,i = 6.9499 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 5.9352. (a) Spatial solute distribution across all domains, including the evaporation region, vehicle, skin layers, and receptor fluid. (b) Solute distribution in the receptor fluid.

Figure 3

Single-parameter K_evap,i sensitivity analysis of 4-Tolunitrile with D_SC,i = 6.9499 × 10⁻¹⁰ cm²·s⁻¹ and P_SCw,i = 5.9352. The black dot denotes the point yielding the minimum MSE among all sampled values, and the corresponding parameter value is annotated directly adjacent to the marker.

Figure 4

Solute distributions under preliminary K_evap,i optimization of 4-Tolunitrile: K_evap,i = 8.3929 × 10⁻¹¹ mol·cm⁻²·s⁻¹. (D_SC,i = 6.9499 × 10⁻¹⁰ cm²·s⁻¹ and P_SCw,i = 5.9352). (a) Spatial solute distribution across all domains, including the evaporation region, vehicle, skin layers, and receptor fluid. (b) Solute distribution in the receptor fluid.

Figure 5

Single-parameter sensitivity analysis of 4-Tolunitrile with K_evap,i fixed at its optimized value (8.3929 × 10⁻¹¹ mol·cm⁻²·s⁻¹). Panels show MSE responses when scanning (a) D_SC,i (P_SCw,i fixed at 5.9352) and (b) P_SCw,i (D_SC,i fixed at 6.9499 × 10⁻¹⁰ cm²·s⁻¹).

Figure 6

Solute distributions under the initial calculated value of Nitrobenzene: K_evap,i = 6.6480 × 10⁻¹⁰ mol·cm⁻²·s⁻¹, D_SC,i = 5.2677 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 4.8645. (a) Spatial solute distribution across all domains, including the evaporation region, vehicle, skin layers, and receptor fluid. (b) Solute distribution in the receptor fluid.

Figure 7

Single-parameter K_evap,i sensitivity analysis of Nitrobenzene with D_SC,i = 5.2677 × 10⁻¹⁰ cm²·s⁻¹ and P_SCw,i = 4.8645. The black dot denotes the point yielding the minimum MSE among all sampled values, and the corresponding parameter value is annotated directly adjacent to the marker.

Figure 8

Single-parameter sensitivity analysis of Nitrobenzene with K_evap,i fixed at its optimized value (8.1174 × 10⁻¹¹ mol·cm⁻²·s⁻¹). Panels show MSE responses when scanning: (a) D_SC,i (P_SCw,i fixed at 4.8645) and (b) P_SCw,i (D_SC,i fixed at 5.2677 × 10⁻¹⁰ cm²·s⁻¹).

Figure 9

Comparison of RF concentration profiles generated under two parameter sets for 4-Tolunitrile and Nitrobenzene. (a) 4-Tolunitrile: (i) The initial prediction (K_evap,i = 7.9348 × 10⁻¹⁰ mol·cm⁻²·s⁻¹, D_SC,i = 6.9499 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 5.9352); (ii) the K_evap,i-calibrated case (K_evap,i = 8.3929 × 10⁻¹¹ mol·cm⁻²·s⁻¹, D_SC,i = 6.9499 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 5.9352). (b) Nitrobenzene: (i) The initial prediction (K_evap,i = 6.6480 × 10⁻¹⁰ mol·cm⁻²·s⁻¹, D_SC,i = 5.2677 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 4.8645); (ii) the K_evap,i-calibrated case (K_evap,i = 8.1174 × 10⁻¹¹ mol·cm⁻²·s⁻¹, D_SC,i = 5.2677 × 10⁻¹⁰ cm²·s⁻¹, and P_SCw,i = 4.8645).

Figure 10

Solute distributions under preliminary K_evap,i optimization of Nitrobenzene: K_evap,i = 8.1174 × 10⁻¹¹ mol·cm⁻²·s⁻¹. (D_SC,i = 5.2677 × 10⁻¹⁰ cm²·s⁻¹ and P_SCw,i = 4.8645). (a) Spatial solute distribution across all domains, including the evaporation region, vehicle, skin layers, and receptor fluid. (b) Solute distribution in the receptor fluid.