In-Silico Modelling of Transdermal Delivery of Macromolecule Drugs Assisted by a Skin Stretching Hypobaric Device

Objectives To develop a simulation model to explore the interplay between mechanical stretch and diffusion of large molecules into the skin under locally applied hypobaric pressure, a novel penetration enhancement method. Methods Finite element method was used to model the skin mechanical deformation and molecular diffusion processes, with validation against in-vitro transdermal permeation experiments. Simulations and experimental data were used together to investigate the transdermal permeation of large molecules under local hypobaric pressure. Results Mechanical simulations resulted in skin stretching and thinning (20%–26% hair follicle diameter increase, and 21%–27% skin thickness reduction). Concentration of dextrans in the stratum corneum was below detection limit with and without hypobaric pressure. Concentrations in viable epidermis and dermis were not affected by hypobaric pressure (approximately 2 μg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bullet$$\end{document}∙ cm−2). Permeation into the receptor fluid was substantially enhanced from below the detection limit at atmospheric pressure to up to 6 μg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bullet$$\end{document}∙ cm−2 under hypobaric pressure. The in-silico simulations compared satisfactorily with the experimental results at atmospheric conditions. Under hypobaric pressure, satisfactory comparison was attained when the diffusion coefficients of dextrans in the skin layers were increased from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}∼ 10 μm2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bullet$$\end{document}∙ s−1 to between 200–500 μm2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bullet$$\end{document}∙ s−1. Conclusions Application of hypobaric pressure induces skin mechanical stretching and enlarges the hair follicle. This enlargement alone cannot satisfactorily explain the increased transdermal permeation into the receptor fluid under hypobaric pressure. The results from the in-silico simulations suggest that the application of hypobaric pressure increases diffusion in the skin, which leads to improved overall transdermal permeation.

Mass fraction of protein ( pro ) and lipids ( lip ) respectively in the stratum corneum Dimensionless

Introduction
The dermal/transdermal drug delivery global market will hit $87.32 billion by 2030 given its advantages with respect to other drug delivery methods such as the oral and the parenteral route [1]. Among these advantages, dermal/transdermal drug delivery avoids gastrointestinal degradation and firstpass metabolism, encouraging patient compliance, and facilitating sustained and controlled local release [2][3][4][5]. Common strategies involve passive diffusion by means of patches and semisolid formulations such as creams, gels, and ointments [6]. These methods are available commercially for drugs with small molecular weight such as nicotine, menthol, or estradiol, among others. The number of drugs that can be administered is, however, limited by the barrier properties of the stratum corneum, which gives rise to low skin permeability for modern drugs used in advanced therapies with molecular weights greater than ca. 500 Da [7]. For instance, hormones have molecular weights of approximately 10 kDa, while vaccines are in the range between 145 and 160 kDa, and liposomes and nanospheres are even bigger. Efficient transdermal delivery of advanced therapies in convenient doses demands, thus, the development of novel technologies for dermal/transdermal drug delivery.
To help by-pass the stratum corneum barrier, minimally invasive techniques have been developed, which present advantages over common passive diffusion methods [8,9]. These techniques provide improved safety, better compliance, painless administration, and faster delivery. Consequently, they are supported by prominent public health organizations including the WHO [10]. Microneedles, for instance, have been extensively studied, although their commercial development remains limited at present [11]. Other methods with a promising future are those using pressure gradients to overcome the barrier properties of the outer skin layer, such as jet injection, which uses a blast of liquid at high pressure to overcome the stratum corneum and deliver an active principle. Recent advances have been demonstrated, too, using a blast of air instead of liquid, showing excellent results [12].
Recently, a novel method using hypobaric instead of hyperbaric pressure has been reported for large molecules. The application of hypobaric pressure on the skin surface allows the formation of a skin dome, which results in skin stretching and that facilitates transdermal administration of large molecules. Successful demonstration of direct delivery of large molecules into the skin both ex vivo and in vivo using the hypobaric device has been reported recently [13]. Also imaging techniques have been used to measure the skin deformation under hypobaric pressure [14]. However, the mechanism of how skin geometry changes enhanced permeation has not been elucidated. In this work, by using a combination of in-vitro experiments and in-silico modelling and simulations, we quantify the effect of mechanical deformation and delve into the mechanisms behind the observed link between mechanical deformation and enhanced permeation under hypobaric pressure. The study contributes to advance the knowledge on the interplay between the mechanical and diffusion characteristics of the skin under local hypobaric pressure.

Experimental
In vitro skin permeation studies were conducted using full thickness porcine ear skin (thickness 1.36 ± 0.16 mm, 10 measurements, obtained from a local abattoir). Figure 1 shows a photograph of the diffusion cell set-up and a schematic of the experimental procedure. The skin was prepared by the removal of any hair and by the careful removal of the subcutaneous fat layer. The skin was then cut into 3.2 cm diameter circles and mounted onto the receptor compartments of upright Franz diffusion cells (2.21 cm 2 inner diameter, 9.5 ml receptor fluid, University of Southampton, UK) with the stratum corneum facing upwards. For the duration of the hypobaric pressure application, the receptor phase consisted of a phosphatebuffered saline (PBS) wetted sponge placed underneath the skin to avoid back flow of any fluid to the donor compartment. The hypobaric chamber was attached to the skin and then the mounted cells were placed on a submersible stirring plate in a pre-heated water bath (Grant Instruments, Cambridge, UK) set at 37°C, to obtain a temperature of 32°C at the skin surface. The studies were initiated by the application of an infinite dose (0.5 ml) of the purified FITC-dextran solutions (500 mg/ml) to the skin surface through the chamber application port. The studies were conducted under control conditions (atmospheric pressure) and hypobaric pressure conditions (-4.5 psig and -8.5 psig) applied for 20 min, 40 min and 1 h. Purified FITC-dextran with measured molecular weights of 31.2 kDa and 95.6 kDa were used [13].
At the end of the permeation studies, the Franz diffusion cells were disassembled, the skin was removed and wiped with water, two tape strips were removed to ensure the removal of any residual formulation. The skin layers were separated as described previously [14]. The stratum corneum was obtained by tape stripping (ca. 20 strips until the skin was translucent) using adhesive tape (Scotch 845 book tape, 3 M, Bracknell, UK). The sponge, which acted as the receptor, and the tape strips were incubated in PBS (pH 7.4) on a shaking water bath overnight at 32°C. The epidermis and dermis were physically separated using a scalpel. Each layer was cut into small pieces and incubated in PBS (pH 7.4) overnight on a shaking water bath at 32°C. These mixtures were then homogenised (Ultra Turrax, Fisher Scientific) and then filtered using 0.45 μm syringe filters.
Quantitative determination of FITC-dextran was achieved by fluorescent spectroscopy, using a Spark plate reader (Tecan Ltd) at an excitation wavelength of 485 ± 20 nm and emission wavelength of 535 ± 20 nm.

In-Silico Modelling of Transdermal Permeation
A baseline transdermal permeation model was developed by adapting an earlier model [15] to the specifics of porcine skin and large molecules using the Finite Element Method (FEM) commercial solver COMSOL Multiphysics 5.5. The baseline model was compared to the experimental in-vitro data at atmospheric pressure for validation. A 2-D geometry was used to represent the entirety of the application area. The 2-D geometry features an infundibulum zone which accounts for all the infundibula in the application area in terms of superficial area and depth into the skin layers. The output from the simulation is the time evolution of the mass in each layer of the skin per unit area of application. This is calculated by integrating the concentration within a specific layer of the skin and then divided by the application area ( A app ). A schematic of the 2-D geometry can be seen in Fig. 2, which includes an explanation on the choice for the horizontal dimensions, and a colour-code for the boundary conditions. The donor with infinite dose was not included in the simulation. Instead, the interface between the donor and the stratum corneum was modelled as a constant concentration boundary (i.e., concentration in the donor corrected by the stratum corneum to donor partition coefficient). The interface between the dermis and the receptor fluid was set as a variable concentration boundary (i.e., specified as the concentration in the receptor corrected by the corresponding partition coefficient). Partition conditions were set at all boundaries separating the skin layers and no flux condition was imposed at the rest of the boundaries. The initial condition was zero concentration in all skin layers and the receptor fluid compartment. The thickness of the layers and depth of the infundibulum correspond to average values reported for porcine ear skin [16]. The model includes the basement membrane at the epidermal-dermal junction as an extra layer in the 2-D geometry, following studies that showed its significant role in hindering the permeation of large molecules [17].
The model for transdermal permeation of large molecules solves the diffusion equation for each layer of the skin [18,19] where c i is the concentration in a specific (the i-th) layer of the skin (e.g., stratum corneum, viable epidermis, et cetera), t denotes time and � ⃗ J corresponds to Fick's diffusive flux: (1) with D i being the diffusion coefficient in the i -th layer of the skin. Diffusion across different layers of the skin is dealt with by introducing the partition coefficients between the layers, as reported previously [15] by using well-established methodology in the literature. The receptor fluid is modelled as a shapeless compartment governed by the equation Equation (2) describes the mass balance on the receptor compartment, where the accumulation term equals the mass input due to the average flux J in the dermis/receptor interface. The concentration in the receptor fluid is c RF , and A app_def is the contact area between the dermis and the receptor fluid. At atmospheric pressure A app_def = A app , while under hypobaric pressure A app_def > A app as the formation of a dome increases this contact area. V RF is the volume of the receptor fluid. The input from the skin to the receptor fluid compartment described by Eq. 2 is the average flux on the dermis/receptor fluid interface J . The concentration on the dermis/receptor fluid interface is imposed as a function of the concentration in the receptor fluid as c RF K De∕w , where K De∕w represents the partition coefficient between the dermis and the PBS solution used in the receptor fluid.
The follicular infundibulum was assumed to be initially filled with the donor PBS solution to mimic the conditions of an in-vitro test. The hair shaft is assumed to be impermeable and was not included in the model geometry.

Partition Coefficients
Partition coefficients were estimated using the established QSPR (quantitative structure-property relationship) models from the literature. For the case of the stratum corneum, the following expression was used [20], which averages the properties of the proteins, lipids, and water: where represents the volume fraction and represents the density. The overall density of the hydrated stratum corneum is calculated using the equation where x denotes mass fraction and denotes the mass of water absorbed per unit mass of dry stratum corneum. Fig. 2 Schematic of the 2-D computational domain for the transdermal permeation simulations at atmospheric pressure. Not to scale. The 2.21 cm 2 application area contains 40 follicular infundibula on average and each follicle has an average infundibulum diameter of 201 μm and hair fibre average diameter of 77.5 μm [16]. The total width selected was 120 μm (119.9 μm excluding the impermeable hair fibre), which accommodates a sufficient number of corneocytes [15]. Keeping the ratio of areas between application area, infundibulum, and hair to be the same as the real anatomy, the width of the infundibulum zone is 0.58 μm and the width of the hair is 0.1 μm.
The viable epidermis and the dermis were assumed to have the behaviour of a hydrogel [21][22][23]. Measurements for dextrans in PBS and hydrogel reported in the literature have been implemented for the partition of the viable epidermis to water K VE∕w and the dermis to water K De∕w partition coefficients [24,25]. This assumption was supported by the established QSPRs in the literature [26], assuming no protein binding ( f u =1) and no ionisation ( f non =1): The partition coefficient between the basement membrane (a lipid bilayer) and water K BM∕w was obtained using the following expression for lipids [20]:

Diffusion Coefficients
The diffusion coefficient of dextrans in the infundibulum was set to be that for water and calculated by the Stokes-Einstein equation at 32°C (dynamic viscosity of water 0.7644 mPa • s). For the stratum corneum and the basement membrane, values which ensured a negligible mass of dextrans in these layers were selected, a phenomenon observed in the experiments. The diffusion parameters for the viable epidermis D VE and dermis D De were selected according to values reported in the literature for hydrogels. Diffusion coefficients with values around 10 m 2 ∕s for dextrans in agarose gel discs have been reported [25]. Other studies show measured diffusion coefficients of 70, 500, and 2000 kDa dextrans in ex-vivo human skin using intradermal microneedles, finding values between 2-10 m 2 ∕s for the viable epidermis, with a significant reduction in the dermal/epidermal junction, and increased again in the reticular dermis to its highest levels (10-20 m 2 ∕s ) [27].

FEM Modelling of the Skin Mechanics
The application of hypobaric pressure on the skin surface results in the formation of a dome, the stress-strain characteristics of which have been previously reported using a validated static simulation model [28]. A variation of the reported hyperelastic skin mechanical deformation simulation set-up including the follicular gaps was developed in COMSOL Multiphysics 5.5. The model solves the quasistatic equilibrium equation where S denotes surface, �⃗ F is the force vector field, and F is the sum of the identity matrix I and the displacement (5) tensor as F = I + ∇� ⃗ u . Initial rest conditions were assumed. A load boundary was used on the hypobaric pressure application area. A restrained movement boundary condition was used on top (except the load boundary) and on each side of the geometry. Free boundaries were used on the bottom boundary and inside the follicular gaps. The one-term Ogden hyperelastic constitutive equation was used where denotes stress, lambda denotes stretch ratio, and = 0.03 MPa and = 13.57 are constants which are fitted from experimental skin tensile tests on porcine skin samples [28].
The model was then used to predict the opening of the hair follicle depending on its position on the dome. The average skin thickness and application area were adapted to the values used in this study (1.36 mm and 2.21 cm 2 respectively). The results from the skin mechanics simulation were used to modify the geometry of the baseline transdermal permeation model when stretching occurs due to the application of hypobaric pressure.

Experimental Results
The experimental results are gathered in Fig. 3, which shows the time evolution of the mass of dextran in each layer of the skin per unit area of application under different hypobaric pressures. Mass of dextran in the stratum corneum below detection limit was expected at atmospheric conditions. However, the results suggest that mechanical deformation did not result in detectable disruption of the stratum corneum barrier, hence the negligible concentration of dextrans in the stratum corneum observed under hypobaric pressure. The application of hypobaric pressure also resulted in little impact in the mass of dextrans in the viable epidermis and dermis but resulted in considerably increased transdermal permeation into the receptor fluid. This suggests that hypobaric pressure has increased the flux in and out of the skin both significantly, resulting in remarkable increase of penetration through the skin but negligible change in skin deposition. Such results also support the increase of diffusion parameters in our simulation presented later. The results presented in Fig. 3 provide evidence of the adequacy of skin stretching using hypobaric pressure as an alternative minimally invasive novel method to facilitate transdermal delivery of macromolecules [29].

Modelling Transdermal Permeation at Atmospheric Pressure: Baseline Case
The simulation parameters for the baseline simulation for the transdermal permeation of dextrans at atmospheric pressure are presented in Table I. Figure 4 shows the comparison between the simulated and the experimental results. The simulation results replicated the trend in the experiments by which greater concentrations of dextrans are observed in the viable epidermis than in the dermis. This points to the role that the basement membrane has as a barrier for diffusion into the dermis [30]. The use of independently derived values for the diffusion and partition coefficients (not calibrated from our experiments) resulted in a model with an adequate prediction capability. Consequently, reasonable agreement between the experimental and simulated results was observed for the baseline case, though the insilico model somewhat underpredicts the kinetics in epidermis and dermis. A grid independence study was carried out, using three levels of grid density with 11,938 elements (Extra Coarse option), 23,930 elements (Normal option), and 30,327 (Extremely Fine option). The results show differences < 1% between the three grids used for the mass in the viable epidermis and dermis after 1 h permeation time for the 31.2 kDa dextran (baseline case). The simulations were conducted with the Normal option grid.

Modelling Transdermal Permeation Under Local Hypobaric Pressure
Simulation results of the mechanical deformation of the skin under hypobaric pressure can be seen in Fig. 5, where a visual of the dome formation colour coded by radial strain (i.e., stretching in the direction of the black arrows) has been included. A grid independence study was performed for an applied hypobaric case of -4.5 psig. Three grid levels were used: 2,027 elements (Normal option), 40,769 elements (Extra Fine option), and 189,956 (Extremely Fine option). Differences < 1% between the three grids were observed for the displacement of the dome apex. The simulations were conducted with the Extra Fine option grid. Hypobaric   [27]. b Using Eq. (6) and a Value for the Octanol to Water Partition Coefficient of K o∕w =0.17 [31]. c Using Eq. (3) and the Following Parameters: pro = 0.1476, lip = 0.0671, w = 0.7853, pro = 1.37 g/cm 3 , lip = 0.9 g/cm 3 , w = 1 g/cm 3 , w pro = 0.77, w lip = 0. 23  pressure gave rise to a positive radial strain distribution across the skin, which resulted in increased diameter of the infundibula. At the dome apex (where stretching is maximum), the increase in diameter was 20% at -4.5 psig and 26% at -8.5 psig. Corresponding changes in area of the skin are gathered in Table II. Thickness reduction of the skin caused by Poisson effect was also maximum at the dome apex (21% overall thickness reduction at -4.5 psig and 27% at -8.5 psig). Accordingly, skin deformation was introduced in the geometry of the transdermal permeation simulations, including the expansion of the application area, and thinning of the skin layers, resulting in the dimensions shown in Fig. 6. Transdermal permeation simulations using the deformed geometry were run while keeping the parameters of the baseline case to check the effect of skin deformation only on the permeation process and the results are gathered in Fig. 7.
The simulated results for the viable epidermis, dermis, and stratum corneum present good agreement with the experiments in Fig. 7, where the geometric changes induced by hypobaric pressure are incorporated into the model while keeping the diffusion and partition parameters to the same value as in the baseline case. However, the simulated data series for the receptor fluid present values close to zero, unlike the experimental results, thus suggesting that the geometry change only cannot explain satisfactorily the increased transdermal permeation observed in the experiments under hypobaric conditions. Some studies in the literature suggest that skin deformation by means of flexing and massage [32] and motion [33] may result in enhanced dermal permeation. The effect of mechanical stress on the diffusion properties is documented in the literature in the field of materials science [34]. Accordingly, we postulate that besides skin mechanical stretching, the application of hypobaric pressure may cause a change in the properties of the skin diffusive barrier leading to increased transdermal permeation.
Based on this postulation, an additional set of transdermal permeation simulations were run on the deformed geometry by increasing the diffusion coefficients. First, we calibrated the diffusion coefficients to the experimental data of 31.2 kDa / Fig. 4 Comparison of simulation results of baseline model against experimental data at atmospheric pressure.   1.3). Finally, the extrapolation from 31.2 kDa to 95.6 kDa was achieved as follows. As the ratio between the diffusivities for the two molecular weights at atmospheric pressure is 1.5 (Table I), the same ratio is assumed regardless of the hypobaric pressure applied. In other words, the diffusivities for 95.6 kDa modules are decreased from the values for 31.2 kDa by a factor of 1.5 in all cases. This approach minimises the need for calibration to only one scenario and allows testing the extrapolation capability of the model. The diffusion parameters obtained from this calibration-extrapolation study are gathered in Table III.
The results of the calibration-extrapolation can be seen in Fig. 8. The simulated results show values close to zero for the mass in the stratum corneum, in agreement with the experiments, while reasonable agreement is seen  for the viable epidermis and dermis. The calibrated and extrapolated diffusion coefficients resulted in increased overall permeability, hence enhanced transdermal permeation, with the mass in the receptor fluid attaining similar levels to that observed in the experiments.

Conclusions
A computer simulation model of the permeation of large molecules into the skin under atmospheric and hypobaric conditions has been developed. In-vitro experimental data were also obtained, which allowed validation of the simulation model at atmospheric pressure. Comparison between experiments and simulations under hypobaric conditions resulted in valuable understanding on the interplay between mechanical and diffusive properties when permeation of large molecules is assisted by a novel method based on hypobaric pressure.
Implementing skin deformation and subsequent follicular aperture in the simulations did not fully explain the increased transdermal dextran permeation observed under hypobaric pressure in the experiments. Reasonable agreement between the experimental data and the simulated mass of dextran permeated into the receptor fluid at hypobaric pressure is observed however when both the skin geometry and the diffusion coefficients of the skin layers in the simulation are modified. Specifically, when the diffusion coefficients are increased considerably, suggesting that substantial change of the skin diffusive barrier is likely to occur under hypobaric pressure along with the obvious skin deformation.
The computer simulation model developed in this work and the conclusions thereof are being used for the design of an innovative needle-free hypobaric device for topical administration of large molecules to the skin. Experimental investigation of the mechanisms of change in the skin diffusive properties will be conducted in the future.

Fig. 8
Comparison between experimental and simulated transdermal permeation data. Deformation is implemented in these simulations and diffusion coefficients are calibrated according to guidelines in Table III to obtain reasonable  agreement with experimental  data