Impedance-Controlled Molecular Transport Across Multilayer Skin Membranes

1. Introduction

Transdermal drug delivery (TDD) is a noninvasive strategy for local and systemic drug administration, providing controlled dosing, avoidance of first-pass metabolism, and improved patient compliance [1,2,3,4,5]. The skin functions as a multilayer membrane, where molecular diffusion is governed by the coupled properties of its layers rather than by the epidermal barrier alone [6,7,8,9,10] (Figure 1). While the stratum corneum (SC) and viable epidermis (VE) are primary barriers, drug retention, delayed clearance, and redistribution within the dermis and hypodermis significantly influence overall transport kinetics [11,12,13,14,15]

Most analytical models represent deeper skin layers using simplified boundary conditions [16,17]. Commonly, a Dirichlet boundary condition is imposed at the viable epidermis–dermis interface, treating the dermis as an ideal sink where drug molecules are instantly removed by the capillary network [1,12,18,19]. While mathematically covenient, these models neglect finite transport, storage [13,18,20], and clearance within the backing layers, and cannot distinguish short-time (diffusively thick) from long-time (diffusively thin) transport regimes.

Alternative approaches introduce interfacial resistances or Robin-type boundary conditions to represent imperfect coupling between the epidermis and dermis [9,10,12,13]. While partially relaxing the ideal sink assumption, these models rely on phenomenological parameters and do not capture the dynamic influence of backing-layer thickness, conductivity, or storage properties.

In this work, we model transdermal transport as a dynamically coupled multilayer diffusion problem [21], rather than as a single barrier with prescribed boundary conditions. The skin is represented by a thin epidermal membrane coupled to a finite diffusive backing layer (dermis and hypodermis), each defined by its intrinsic transport conductivity, storage capacity, and thickness.

By introducing diffusion impedance and transfer-matrix analysis, the interaction between skin layers is obtained directly from the continuity of concentration and diffusive flux at internal interfaces. This eliminates artificial interfacial boundary conditions while retaining only physically real external boundary interactions.

In this framework, the effect of deeper layers is naturally captured by diffusion impedance, relating interfacial concentration and flux in the Laplace domain. This provides a clear measure of how backing-layer properties and distal boundary conditions control effective transport across the epidermis. Short-time dynamics are dominated by backing-layer geometry and boundary effects, while long-time behavior depends mainly on conductive and storage properties, illustrating the regime-dependent impact of macroscopic skin heterogeneity on cumulative drug uptake.

For tractable analytical results, the multilayer skin system is represented by a simplified two-layer model: the stratum corneum and viable epidermis form a thin effective membrane, while the dermis and hypodermis are combined into a finite diffusive backing layer. The drug patch serves as an external source, and the distal boundary represents the interface with the systemic receiver, which can be modeled using absorbing or reflecting conditions depending on the physiological scenario.

While the skin shows pronounced microstructural heterogeneity [14,17,18], especially within the stratum corneum, this study deliberately focuses on the macroscopic layered organization. This isolates the role of deeper layers in transdermal transport through measurable quantities, such as cumulative drug uptake, while avoiding the additional complexity of microheterogeneity [22,23,24,25,26]. Extensions to anomalous and memory-dependent transport processes can be naturally incorporated within the same impedance-based framework [27].

The remainder of this paper is organized as follows. Section 2 presents the theoretical formulation and introduces diffusion impedance and transfer-matrix analysis. Section 3 derives analytical expressions for cumulative drug uptake and discusses the relevant transport regimes. The main conclusions are summarized in Section 4 and Section 5.

2. Theory: The Concept of Diffusion Impedance

2.1. Model Description

We base our analysis on Fickian diffusion and focus on the macroscopic layered structure of the skin, neglecting the heterogeneous and fractal microstructure of individual layers. This deliberate choice allows the role of macroscopic (layered) heterogeneity to be isolated and analyzed independently of microstructural effects, which may give rise to anomalous or memory-dependent transport regimes [17]. To quantify the influence of finite dermal conductivity and storage capacity on transdermal transport, we consider a two-layer analytical model. The stratum corneum (SC) and viable epidermis (VE) are treated as a thin effective barrier, while the dermis and hypodermis are represented as a finite backing layer.

The transdermal patch acts as an external source at the membrane surface, while the distal boundary of the backing layer corresponds to the interface with the systemic receiver, which can be modeled as either an ideal sink or a perfectly reflecting boundary (see Figure 2). Continuity of both concentration and diffusive flux is imposed at the interface between the barrier and the backing layer, ensuring physically consistent coupling. Within this formulation, the multilayer system can be viewed as a sequence of coupled transport elements, whose interaction is most naturally described using transfer matrices and impedance concepts [21,27].

The initial condition assumes no drug molecules present in the system prior to patch application, consistent with standard TDD experiments.

The mathematical formulation of the problem is given by the following system of partial differential equations.

For layer 1 (viable epidermis, VE):

$${\displaystyle \frac{{\partial }^2{n}_1(x,t)}{\partial x^2}}-{\displaystyle \frac{{\gamma }_1}{D_1}}{\displaystyle \frac{\partial n_1(x,t)}{\partial t}}=0,$$

(1)

$$j_1(x,t)=-D_1{\displaystyle \frac{\partial n(x,t)}{\partial x}},$$

(2)

For layer 2 (dermis + hypodermis):

$${\displaystyle \frac{{\partial }^2{n}_2(x,t)}{\partial x^2}}-{\displaystyle \frac{{\gamma }_2}{D_2}}{\displaystyle \frac{\partial n_2(x,t)}{\partial t}}=0,$$

(3)

$$j_2(x,t)=-D_2{\displaystyle \frac{\partial n_2(x,t)}{\partial x}},$$

(4)

where γ1 and γ2 denote the mass storage capacities of the viable epidermis and dermis, respectively, and D1 and D2 are the corresponding diffusion coefficients.

The donor phase is described by a constant concentration, which is applied at the moment the transdermal patch is placed on the skin [1]:

$$n_1(x=0,t)=n_{0}h(t),$$

(5)

At the distal boundary of the backing layer, either an ideal receiver (ideal sink) is assumed:

$$n_2(x=l_1+l_2,t)=0,$$

(6)

or a completely impermeable boundary (adiabatic condition):

$$j_2(x=l_1+l_2,t)=0$$

(7)

At the heterointerface between the barrier and the backing layer, continuity of concentration and diffusive flux is imposed:

$$n_1(x=l_1,t)=n_2(x=l_1,t)$$

(8)

$$j_1(x=l_1,t)=j_2(x=l_1,t),$$

(9)

which represents the most physically justified choice for coupling between layers.

It is further assumed that no drug molecules are present in the system prior to patch application, leading to zero initial conditions:

$$n_1(x,t=0)=n_2(x,t=0)=0,$$

(10)

2.2. Laplace Formulation and Diffusion Impedance

Applying the Laplace transform to the governing equations yields the following system describing the two-layer structure:

$${\displaystyle \frac{d^2{\overline{n}}_k(x,s)}{d{x}^2}}-{\overline{\sigma }}_k^2\overline{n}(x,s)=0,$$

(11)

$${\overline{j}}_k(x,s)=-{\displaystyle \frac{1}{{\overline{\sigma }}_k{\overline{Z}}_{ck}}}{\displaystyle \frac{d{\overline{n}}_k(x,s)}{dx}},$$

(12)

where

$${\overline{\sigma }}_k=\sqrt{{\displaystyle \frac{{\gamma }_k}{D_k}}}\sqrt{s},$$

(13)

$${\overline{Z}}_{ck}={\displaystyle \frac{1}{\sqrt{{\gamma }_k{D}_k}}}{\displaystyle \frac{1}{\sqrt{s}}},$$

(14)

with k=1,2.

In this representation, each layer is characterized by an effective propagation coefficient and a characteristic diffusion impedance, allowing the multilayer structure to be treated analogously to a cascaded transport network [27].

The boundary condition at the donor interface corresponds to a constant drug concentration applied by the patch:

$${\overline{n}}_1(x=0,s)=n_0{\displaystyle \frac{1}{s}},$$

(15)

At the distal side of the backing layer, either an ideal sink (Dirichlet) or a reflecting (Neumann) boundary condition can be used,

$${\overline{n}}_2(x=l_1+l_2,s)=0,$$

(16)

$${\overline{j}}_2(x=l_1+l_2,s)=0,$$

(17)

while continuity conditions are enforced at the membrane–backing interface:

$${\overline{n}}_1(x={l_1}^{-},s)={\overline{n}}_2(x={l_1}^{+},s),$$

(18)

$${\overline{j}}_1(x={l_1}^{-},s)={\overline{j}}_2(x={l_1}^{+},s),$$

(19)

Solving this system provides analytical expressions for concentration and flux at the interface:

For x<l₁ we obtained

$$\left[\begin{array}{c}{\overline{n}}_1(x,s)\\ {\overline{j}}_1(x,s)\end{array}\right]=\left[\begin{array}{cc}ch({\overline{\sigma }}_{1}x)& -{\overline{Z}}_{c1}sh({\overline{\sigma }}_{1}x)\\ -{\displaystyle \frac{1}{{\overline{Z}}_{c1}}}sh({\overline{\sigma }}_{1}x)& ch({\overline{\sigma }}_{1}x)\end{array}\right]\left[\begin{array}{c}{n}_0\\ \overline{j}(0,s)\end{array}\right]$$

(20)

For x>l₁

$$\left[\begin{array}{c}{\overline{n}}_2(x,s)\\ {\overline{j}}_2(x,s)\end{array}\right]=\left[\begin{array}{cc}ch({\overline{\sigma }}_2(x-l_1))& -{\overline{Z}}_{c2}sh({\overline{\sigma }}_2(x-l_1))\\ -{\displaystyle \frac{1}{{\overline{Z}}_{c2}}}sh({\overline{\sigma }}_2(x-l_1))& ch({\overline{\sigma }}_2(x-l_1))\end{array}\right]\left[\begin{array}{c}{\overline{n}}_1(l_1,s)\\ {\overline{j}}_1(l_1,s)\end{array}\right]$$

(21)

To explicitly account for the effect of the backing, we define the diffusion impedance of the backing layer as

$${\overline{Z}}_2={\displaystyle \frac{{\overline{n}}_1\left(l_1,s\right)}{{\overline{j}}_1(l_1,s)}},$$

(22)

where $\overline{n}\left(l_1,s\right)$ and $\overline{j}(l_1,s)$ are the Laplace-domain concentration and flux at the interface between the membrane and the backing, respectively.

This definition formally encapsulates the entire response of the backing layer—including its material properties, thickness, and distal boundary condition—into a single impedance function.

Substituting Eq. (22) into Eq. (21) and solving the resulting transfer-matrix equation yields with boundary conditions Eqs (15)-(17), the total flux and concentration at the barrier–backing interface:

$$\overline{j}(l_1,s)=n_0{\displaystyle \frac{1}{s}}{\displaystyle \frac{1}{{\overline{Z}}_{c1}sh({\overline{\sigma }}_1{l}_1)+{\overline{Z}}_{2}ch({\overline{\sigma }}_1{l}_1)}},$$

(23)

$$\overline{n}(l_1,s)=n_0{\displaystyle \frac{1}{s}}{\displaystyle \frac{{\overline{Z}}_2}{{\overline{Z}}_{c1}sh({\overline{\sigma }}_1{l}_1)+{\overline{Z}}_{2}ch({\overline{\sigma }}_1{l}_1)}},$$

(24)

The backing impedance can be obtained by solving the differential equation Eq (21) with boundary condition at the distal side of deeper layer of skin, given by (Eq. (16) for reflection or by Eq. (17) for antireflection condition.

By applying Eq. (21) and the Dirichlet boundary condition (reflective) at the distal side of the backing layer (Eq. (16)), we obtain the following expression for the diffusion impedance of the backing

$${\overline{Z}}_2={\displaystyle \frac{{\overline{n}}_1(l_1,s)}{{\overline{j}}_1(l_1,s)}}={\overline{Z}}_{c2}{\displaystyle \frac{sh({\overline{\sigma }}_2{l}_2)}{ch({\overline{\sigma }}_2{l}_2)}}$$

(25)

Alternatively, by applying Eq. (21) together with the Neumann boundary condition (antireflective)(Eq. (17)), the diffusion impedance becomes

$${\overline{Z}}_2={\displaystyle \frac{{\overline{n}}_1(l_1,s)}{{\overline{j}}_1(l_1,s)}}={\overline{Z}}_{c2}{\displaystyle \frac{ch({\overline{\sigma }}_2{l}_2)}{sh({\overline{\sigma }}_2{l}_2)}},$$

(26)

It is important to note that both the characteristic diffusion impedance and the effective propagation coefficient depend on the transport properties of the backing layer, such as conductivity and mass storage capacity (Eqs. (13)–(14)). In cases where the relationship between thermodynamic scalar fields and their fluxes deviates from classical Fickian behavior—such as in subdiffusive, fractional, or wave-like transport processes [22,23,24,25,26]—these quantities acquire modified frequency dependence determined by the fractional order of differentiation or relaxation times associated with inertial memory effects.

This observation highlights the generality of the proposed framework: by appropriately redefining the coefficient of propagation (Eq 13) and characteristic impedance (Eq 14), the same analytical structure remains valid for fractional, subdiffusive, or wave-like transport models, without altering the macroscopic transfer-matrix formulation [18,27].

Although Robin boundary conditions formally introduce an interfacial impedance 1/h, where h is commonly related to partition and mass transfer coefficients [9,10], this impedance is assumed to be constant. As a result, the Robin formulation cannot distinguish whether this effective resistance originates from a finite backing layer with a distal boundary or from a semi-infinite diffusive medium. In both cases, the entire backing region is reduced to the same effective resistance h/D₂, thereby obscuring the influence of layer geometry and distal boundary conditions on the transport process.

3. Cumulative Drug Amount and Effective Transmission

A key quantity in transdermal drug delivery (TDD) applications is the cumulative amount of drug passing through the epidermal barrier, defined as the time integral of the diffusive flux at the VE–dermis interface [1,12,17]:

$$Q(t)=A{\displaystyle \underset{0}{\overset{t}{\int }}j(x=l_1,t)},$$

(27)

where A denotes the cross-sectional area of the patch application site. The quantity Q(t) has the dimension of substance amount and represents the total number of drug molecules that have crossed the epidermal barrier up to time t.

To evaluate this quantity for the two-layer system shown in Figure 2, we first apply the Laplace transform to Eq. (27). The spectral representation of the cumulative amount of drug transmitted across the barrier is given by:

$$\overline{Q}(s)={\displaystyle \frac{A}{s}}{\overline{j}}_1(x=l_1,s),$$

(28)

Substituting the interface flux obtained in Eq. (23) into Eq. (28) yields:

$$\overline{Q}(s)=A{n}_0{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{1}{{\overline{Z}}_{c1}sh({\overline{\sigma }}_1{l}_1)+{\overline{Z}}_{2}ch({\overline{\sigma }}_1{l}_1)}}$$

(29)

As evident from Eq. (29), the cumulative amount of drug entering the vascularized layers of the skin depends not only on the properties of the epidermal barrier (thickness, conductivity, and storage capacity), but also on the transport properties of the backing layer, fully encoded exclusively through its diffusion impedance.

The inverse Laplace transform provides the time-dependent cumulative drug amount Q(t), which naturally captures the influence of layer thicknesses, transport parameters, and distal boundary conditions, without introducing artificial interfacial resistances.

3.1. Ideal Receiver as a Limiting Case

If the backing layer is assumed to act as an ideal receiver (D2→∞), the corresponding diffusion impedance vanishes, ${\overline{Z}}_{c2}\to 0$ (Eq. 14), and Eq. (29) reduces to

$$\overline{Q}(s)=n_0{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{1}{{\overline{Z}}_{c1}sh({\overline{\sigma }}_1{l}_1)}},$$

(30)

This is equivalent to imposing a Dirichlet boundary condition at the VE–dermis interface [1]. In this regime, Q(t) initially increases linearly and then gradually saturates, with a characteristic time scale determined by the barrier conductivity, storage capacity, and thickness. In Figure 3, this corresponds to the linear Q1(t) curve.

3.2. Diffusively Thin and Diffusively Thick Layer Regimes

The behavior of the backing layer depends on the ratio of layer thickness to characteristic diffusion length. A layer is diffusively thin when $l_k<<{\mu }_k$ and diffusively thick when the opposite inequality holds $l_k>{\mu }_k$, where the characteristic diffusion length is determined by (see Eq. (13):

$${\mu }_k={\displaystyle \frac{1}{\mathrm{Re}\left\{{\overline{\sigma }}_k\right\}}},$$

(31)

Physically, this distinction compares the layer thickness to the characteristic distance over which concentration disturbances can propagate during the observation time.

Considering that

$${\displaystyle \frac{sh(\sqrt{x})}{\sqrt{x}}}=1+{\displaystyle \frac{x}{3!}}+{\displaystyle \frac{x^2}{5!}}+\mathrm{...}=1+{\displaystyle \frac{x}{6}}+{\displaystyle \frac{x^2}{120}}+\mathrm{...},$$

(32)

$$ch(\sqrt{x})=1+{\displaystyle \frac{x}{2!}}+{\displaystyle \frac{x^2}{4!}}+\mathrm{...}=1+{\displaystyle \frac{x}{2}}+{\displaystyle \frac{x^2}{24}}+\mathrm{...}$$

(33)

$${\displaystyle \frac{th(\sqrt{x})}{\sqrt{x}}}={\displaystyle \frac{1+\frac{x}{3!}+\frac{x^2}{5!}+\mathrm{...}}{1+\frac{x}{2!}+\frac{x^2}{4!}+\mathrm{...}}}={\displaystyle \frac{1+\frac{x}{6}+\frac{x^2}{120}+\mathrm{...}}{1+\frac{x}{2}+\frac{x^2}{24}+\mathrm{...}}}$$

(34)

$${\displaystyle \frac{sh(\sqrt{x})}{\sqrt{x}}}\approx 1,$$

(35)

$$ch(\sqrt{x})\approx 1+{\displaystyle \frac{x}{2}},$$

(36)

$${\displaystyle \frac{th(\sqrt{x})}{\sqrt{x}}}\approx {\displaystyle \frac{1}{1+\frac{x}{2}}},$$

(37)

for geometrically thin layers satisfying $l_k<<{\mu }_k\sqrt{6}$the following approximations apply:

In contrast, for geometrically thick layers, $l_k>>{\mu }_k$

$$ch(\sqrt{x})\approx sh(\sqrt{x})$$

(38)

can be assumed.

3.3. Application to Epidermal Barrier and Backing Layer

For typical skin parameters [13,28,29] epidermis (VE) is always diffusively thin over the relevant time scale while dermis + hypodermis (backing layer) behaves as: a) thick at short times (for high harmonics) meaning distal boundaries irrelevant, b) thin at long times (for low harmonics) where distal boundaries influence Q(t).

In the short-time limit, the backing-layer diffusion impedance reduces to its characteristic value ${\overline{Z}}_{c2}={\overline{Z}}_2$ (Eqs. 25, 26), and Eq. (29) simplifies to

$$\overline{Q}(s)=A{n}_0{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{1}{{\overline{Z}}_{c1}sh({\overline{\sigma }}_1{l}_1)+{\overline{Z}}_{c2}ch({\overline{\sigma }}_1{l}_1)}}$$

(39)

Equation (39) applies to a backing layer of finite thickness in the quasi-static regime. At high frequencies (short-time limit), the backing layer becomes diffusively thin and the distal boundary conditions at the skin–systemic interface become explicitly visible through Eqs. (25) and (26).

3.4. Characteristic Transport Regimes

Assuming a geometrically thin epidermal barrier, four characteristic transport regimes can be identified.

1)

Backing layer as an ideal receiver

$$\overline{Q}(s)=n_{0}A{\displaystyle \frac{D_1{\gamma }_1}{l_1}}{\displaystyle \frac{1}{s^2}},$$

(40)

Substituting Eqs. (35) into Eq. (30), and using Eqs. (13) and (14), the spectral function of the cumulative amount becomes

Applying inverse Laplace transformation on Eq (40) yields [30]

$$Q(t)=n_{0}A{\displaystyle \frac{D_1{\gamma }_1}{l_1}}th(t),$$

(41)

As seen from Eqs. (40) and (41), in this case the cumulative drug amount is independent of both the backing-layer thickness and the distal boundary conditions. It depends exclusively on the barrier properties. The cumulative amount initially increases linearly and then gradually saturates, with a characteristic time scale determined by the barrier conductivity, storage capacity, and thickness.

2)

Diffusively thick backing layer (long-time limit)

Substituting Eqs. (35)–(38) into Eq. (39) and using Eqs. (13) and (14), the cumulative amount is described by

$$\overline{Q}(s)=n_{0}A{a}_2{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{\sqrt{s}}{s+b_1\sqrt{s}+b_2}},$$

(42)

with coefficients defined in

$$a_2={\displaystyle \frac{2{\gamma }_2{D}_1\sqrt{D_2}}{l_1^2}},b_1={\displaystyle \frac{2{\gamma }_2\sqrt{D_2}}{{\gamma }_1{l}_1}},b_2={\displaystyle \frac{2D_1}{l_1^2}}$$

(43)

Due to the presence of irrational poles and zeros, the inverse Laplace transform is evaluated via the substitution $p=\sqrt{s}$[31,32]

$$\overline{Q}(s)=n_{0}A{a}_2{\displaystyle \frac{1}{p^3}}{\displaystyle \frac{1}{p^2+b_{1}p+b_2}},$$

(44)

leading to

$$Q(t)=n_{0}A{a}_2\left[{\displaystyle \frac{K_1}{\sqrt{\pi t}}}+K_2+K_5{e}^{{\lambda }_1^{2}t}erfc\left({\lambda }_1\sqrt{t}\right)-K_6{e}^{{\lambda }_2^{2}t}erfc\left({\lambda }_2\sqrt{t}\right)\right],$$

(45)

with coefficients given by Eqs. (44) and (45)

$$K_5={\displaystyle \frac{K_3\left(-{\lambda }_1\right)+K_4}{{\lambda }_2-{\lambda }_1}},K_6={\displaystyle \frac{K_3\left(-{\lambda }_2\right)+K_4}{{\lambda }_2-{\lambda }_1}}$$

(46)

$$K_1=-{\displaystyle \frac{2b_1}{b_2}},K_2={\displaystyle \frac{2}{b_2}},K_3={\displaystyle \frac{2b_1}{b_2^2}},K_4={\displaystyle \frac{2\left(b_1^2-b_2\right)}{b_2^2}},{\lambda }_{1/2}={\displaystyle \frac{b_1\pm \sqrt{b_1^2-4b_2}}{2}},$$

(47)

In this regime, the cumulative drug amount exhibits a coupled dependence on both the barrier and backing-layer properties, while remaining insensitive to the distal boundary condition. Whether the distal boundary is an ideal sink or a reflecting interface becomes irrelevant when the backing layer is diffusively thick.

3)

Diffusively thin backing layer with ideal receiver (Dirichlet boundary) -short-time limit

Starting from Eqs. (25) and (29), and applying the expansions (Eqs. 35–38), the spectral function becomes

$$\overline{Q}(s)=n_{0}A{\displaystyle \frac{D_1{\gamma }_1}{l_1}}{\displaystyle \frac{1}{s^2}}{\displaystyle \frac{s+a_3}{s+b_3}},$$

(48)

$$a_3={\displaystyle \frac{D_2{\gamma }_2}{l_1{l}_2{\gamma }_1}},b_3={\displaystyle \frac{D_2{\gamma }_2{l}_1+D_1{\gamma }_1{l}_2}{{\gamma }_1{l}_2{l}_1^2}},$$

(49)

with parameters defined in Eq (49)

Inverse Laplace transformation performed by method of partial fraction [33] yields

$$Q(t)=n_{0}A{\displaystyle \frac{D_1{\gamma }_1}{l_1}}\left[{\displaystyle \frac{a_3}{b_3}}t+{\displaystyle \frac{a_3-b_3}{b_3^2}}\left(e^{-b_{3}t}-1\right)\right]h(t),$$

(50)

In this case, the cumulative drug amount doesn’t shows linear increase as it is case if the backing layer is ideal receiver (Eq 41). The slope and the characteristic time scale of increase depend on the transport properties of both the barrier and the backing layer (Eq 49).

4)

Diffusively thin backing layer with insulating distal boundary ((Neumann boundary)-short tim limit

Starting from Eqs. (26) and (29), and by applying Eqs (35)-(37), the corresponding spectral function becomes

$$\overline{Q}(s)=n_{0}A{\displaystyle \frac{2D_1{D}_2}{{\gamma }_1{l}_2{l}_1^2}}{\displaystyle \frac{1}{s}}{\displaystyle \frac{1}{s^2+b_{4}s+b_5}},$$

(51)

with coefficients

$$b_4={\displaystyle \frac{l_1^2{\gamma }_1{D}_2+2l_1{l}_2{\gamma }_2{D}_2+l_2^2{\gamma }_2{D}_2}{l_1^2{l}_2^2{\gamma }_1{\gamma }_2}},b_5={\displaystyle \frac{2D_1{D}_2}{l_1^2{l}_2^2{\gamma }_1{\gamma }_2}},$$

(52)

The inverse transformation of Eq (51) yields [34]

$$Q(t)=n_{0}A{\displaystyle \frac{2D_1{D}_2}{{\gamma }_1{l}_2{l}_1^2}}\left[{\displaystyle \frac{1}{{\alpha }^2+{\beta }^2}}+{\displaystyle \frac{1}{\beta \sqrt{{\alpha }^2+{\beta }^2}}}{e}^{-\alpha t}\sin (\beta t-\varphi )\right]h(t),$$

(53)

$$\alpha ={\displaystyle \frac{b_4}{2}},\beta =\sqrt{b_5-{\left({\displaystyle \frac{b_4}{2}}\right)}^2},\varphi =arctg\left({\displaystyle \frac{\beta }{-\alpha }}\right)$$

(54)

Here, the cumulative drug amount exhibits a markedly different short-time behavior compared to case (3), reflecting the decisive influence of the distal boundary condition when the backing layer is diffusively thin.

In Figure 3, we illustrate the cumulative drug amount predicted by the derived model using the parameters listed in Table 1 [13,28,29,35], assuming n₀=1 and A=1.

To quantitatively assess how the present approach compares with commonly used Robin-type boundary approximations, we additionally computed the cumulative drug amount obtained when the finite backing layer is replaced by an equivalent Robin boundary condition with h=D₂/l₂. The comparison, shown in Figure 4, demonstrates that Robin conditions can reproduce either the short-time or the long-time behavior, but fail to capture the full temporal dynamics governed by the true backing-layer diffusion impedance.

4. Summary

4.1. Regime Dependence of Cumulative Drug Transport

The temporal behavior of the cumulative drug amount Q(t) is determined by whether the backing layer behaves as diffusively thick or diffusively thin relative to the observation time.

When the backing layer is diffusively thick (short-time regime), the distal boundary condition does not influence the evolution of Q(t). Comparison of Eqs (41) and (45) shows that, in this regime, the diffusivity and storage properties of the backing layer determine both the shape and magnitude of the cumulative drug amount, while the skin–systemic interface remains effectively invisible (Figure 3a). The presence of an ideal sink (Dirichlet boundary condition) behind a thin barrier leads to a slower initial growth of Q(t).

As time progresses, the backing layer become diffusively thin (long-time regime). The situation then reverses: the backing layer no longer screens the influence of the distal boundary. As demonstrated by comparison of Eqs. (51) and (53), the distal boundary condition now plays a dominant role in shaping the temporal evolution of Q(t) (Figure 3).

Robin boundary conditions, frequently introduced to account for backing-layer effects, cannot fully reproduce this regime-dependent behavior (Figure 4). By reducing the entire backing layer to a constant effective surface resistance, they obscure the influence of both its geometry and the distal boundary, yielding results that are closest—but not equivalent—to those described by Eq. (53).

4.2. Practical Implications for Transdermal Drug Delivery

The analysis reveals that:

• The effectiveness of the epidermal barrier is not an intrinsic property, but is dynamically controlled by the transport properties of the backing layer.
• Early-time drug delivery is governed primarily by the diffusivity and thickness of the backing layer, while being insensitive to the distal boundary condition.
• Long-time cumulative uptake becomes highly sensitive to the distal boundary condition once the backing layer becomes diffusively thin.

These findings provide a physically grounded basis for designing impedance-controlled transdermal delivery systems without introducing phenomenological interfacial parameters.

5. Conclusions

In this work, we developed an analytical framework for transdermal drug delivery across multilayer skin membranes by modeling the system as a thin epidermal barrier dynamically coupled to a finite diffusive backing layer representing the dermis–hypodermis composite. The key conceptual advance is the treatment of deeper skin layers not as prescribed boundary conditions, but as active transport elements whose influence emerges naturally from their intrinsic conductive, capacitive, and geometric properties.

By introducing the concept of diffusion impedance, the interaction between the epidermal barrier and the backing layer is derived directly from the continuity of concentration and diffusive flux at internal interfaces. This eliminates the need for ad hoc interfacial conditions and provides a physically consistent description in which only true external boundary interactions are imposed. As a result, the effective transport behavior of the epidermis becomes a dynamic quantity governed by the response of the underlying layers rather than a fixed material parameter.

The analytical solutions reveal a clear distinction between transport regimes governed by macroscopic skin heterogeneity. In the short-time regime, when the backing layer behaves as diffusively thick, transport is controlled by backing-layer transport properties while distal boundaries remain effectively invisible. In the long-time regime, when the backing layer becomes diffusively thin, the influence of the distal boundary condition becomes dominant, leading to behavior that cannot be captured by classical boundary-condition-based models.

Within this framework, commonly used Dirichlet and Robin formulations arise naturally as limiting cases. An ideal sink corresponds to vanishing backing-layer impedance, reproducing the standard Dirichlet approximation. Robin-type conditions, while often employed to mimic backing-layer effects, are shown to be insufficient for capturing the regime-dependent impedance behavior, particularly when backing-layer geometry and distal boundary conditions play a significant role.

Although the present study focuses on classical diffusion and macroscopic layer heterogeneity, the formulation is inherently extensible. Because diffusion impedance is expressed through transport conductivity and characteristic impedance, the same analytical structure applies to anomalous, fractional, or memory-dependent diffusion models, in which these quantities acquire modified frequency dependence.

Overall, the proposed impedance-based framework provides a unified, physically grounded, and analytically tractable approach for describing transdermal transport across multilayer skin membranes. By directly linking cumulative drug uptake to backing-layer transport properties and boundary interactions, the model offers a robust basis for interpreting experimental data and for the rational design of impedance-controlled transdermal delivery systems.