A two-dimensional mathematical model of percutaneous drug absorption

Abstract

Background

When a drug is applied on the skin surface, the concentration of the drug accumulated in the skin and the amount of the drug eliminated into the blood vessel depend on the value of a parameter, r. The values of r depend on the amount of diffusion and the normalized skin-capillary clearence. It is defined as the ratio of the steady-state drug concentration at the skin-capillary boundary to that at the skin-surface in one-dimensional models. The present paper studies the effect of the parameter values, when the region of contact of the skin with the drug, is a line segment on the skin surface.

Methods

Though a simple one-dimensional model is often useful to describe percutaneous drug absorption, it may be better represented by multi-dimensional models. A two-dimensional mathematical model is developed for percutaneous absorption of a drug, which may be used when the diffusion of the drug in the direction parallel to the skin surface must be examined, as well as in the direction into the skin, examined in one-dimensional models. This model consists of a linear second-order parabolic equation with appropriate initial conditions and boundary conditions. These boundary conditions are of Dirichlet type, Neumann type or Robin type. A finite-difference method which maintains second-order accuracy in space along the boundary, is developed to solve the parabolic equation. Extrapolation in time is applied to improve the accuracy in time. Solution of the parabolic equation gives the concentration of the drug in the skin at a given time.

Results

Simulation of the numerical methods described is carried out with various values of the parameter r. The illustrations are given in the form of figures.

Conclusion

Based on the values of r, conclusions are drawn about (1) the flow rate of the drug, (2) the flux and the cumulative amount of drug eliminated into the receptor cell, (3) the steady-state value of the flux, (4) the time to reach the steady-state value of the flux and (5) the optimal value of r, which gives the maximum absorption of the drug. The paper gives valuable information which can be obtained by this two-dimensional model, that cannot be obtained with one-dimensional models. Thus this model improves upon the much simpler one-dimensional models. Some future directions of the work based on this model and the one-dimensional non-linear models that exist in the literature, are also discussed.

Background

A non-linear dual-sorption model of percutaneous drug absorption is described in [2]. Another mathematical model of percutaneous drug absorption and its exact solution is described in [3]. In [4], the non-linear percutaneous permeation kinetics of timolol is studied in vitro with human cadaver skin. A model for a suspension with a finite dissolution rate is solved numerically in [5]. In [6], a mathematical model is developed for percutaneous drug absorption with regular applications of the drug. All these models are one-dimensional and the region of contact of the skin with the drug, is a single point, say x = 0, where x measures the distance into the skin.

To the authors' knowledge, no two-dimensional mathematical model of percutaneous drug absorption exists in the literature, where the region of contact is a line segment, say x = 0, 0 ≤ y ≤ L _c . The purpose of this paper is to model the above situation mathematically and to obtain an efficient finite-difference method to solve it numerically. The usual loss of accuracy along the boundary, where the boundary conditions are of mixed type, is dealt with in a novel approach.

Methods

Model equations

It is an established fact that the stratum corneum, viable epidermis and upper dermis have distinctly different transport and partitioning characteristics, and diffusion through the dermis could be a critical determinant of the percutaneous absorption profiles at least for some kinds of drugs. The model is developed based on the classical view that a drug molecule is rapidly removed by the capillary once it is partitioned into the dermis. Therefore, the skin-capillary boundary may be, anatomically, interpreted as the stratum corneum-dermis boundary.

Let the drug be applied as an ointment on the skin-surface, say, {(0, y) : 0 ≤ y ≤ L _c }, at time t = 0. Let the thickness (distance between the skin-surface and skin-capillary boundary) of the skin be unity. Assuming that skin is an isotropic medium, that is, the diffusivity is the same in the x and y directions, the concentration c(x, y, t) of the drug is given by Fick's Second Law of Diffusion, namely

where L _d , L _u and L _c are positive real numbers. Assuming that the concentration at the uppermost epidermis {(0, y, t) : 0 ≤ y ≤ L _c } is maintained at unity, the boundary conditions for the model are given by

where

.

In the above notation, k _c is the normalized skin-capillary clearence, k _d is the diffusion parameter and r is the concentration at the skin-capillary boundary when the concentration at the uppermost epidermis is regarded as unity (see [3]). Assuming that there is no drug in the skin before the application, the initial distribution associated with the PDE is

c(x, y, 0) = 0, 0 ≤ x ≤ 1, -L _d ≤ y ≤ L _u     (8)

The flux J = J(t), through the skin to the receptor site, per unit area, is given by

(Fick's First Law of Diffusion), so that, from (7),

The cumulative amount of drug eliminated into the receptor cell per unit area at time τ (see [2]) is

Numerical methods

The approach to be adopted in obtaining a numerical solution, is the method of lines in which the initial/boundary-value problem to be solved is transformed into a first-order initial-value problem. The region associated with the skin is given by R _s = {(x, y, t) 0 ≤ x ≤ 1, -L _d ≤ y ≤ L _u , t > 0}. Let N be a positive integer, and H = 1/(N + 1). Also it is assumed that L _c = Q _C H, L _u = Q _u H, L _d = Q _d H, where Q _c , Q _u and Q _d are positive integers. Superimpose on R _s a square grid with mesh length H > 0 and let k > 0 be a constant time step. The grid points are given by (x _l , y _m , t _j ); x _l = lH, l = 0(1)(N + 1), y _m = mH, m = -Q _d (1)Q _u and t _j = jk, j = 0,1,2,.... The finite-difference solution which approximates the solution c(x, y, t) of the two-dimensional parabolic equation (1), is sought at each mesh point (x _l , y _m , t _j ) in the region [R _s -{0, y), 0 ≤ y ≤ L _c }] × t > 0. For notational simplicity, denote

At any time-level

or , the {(Q _d + Q _u + 1)(N + 2) - (Q _c + 1)} elements will be ordered in rows parallel to the x-axis and in vector form, and will be denoted by C(t) and Q, where Q is due to the boundary conditions and is independent of time. For W= C and Q, and or , denote

In the above notation n takes the value 0 if m ≠ 0(1)Q _c and takes the value 1 otherwise. The vector C(t) is a vector of unknowns and Q is given by by , for m = 0(1)Q _c . All other are zero. In the following discretizations, (x _l , y _m ) and c _l,m , respectively, denote (x _l ,y _m ,t _j ) and . At any time-level t = t _j , the differential equation (1) can be discretized as

If

and are respectively, suitable approximations to c _xx;l,m , and c _yy;l,m , the above discretization can be written as

Discretization at ( x ₀ , y _m ), m = - Q _d :-

The boundary condition (5) gives c _0,m-1= c _0,m+1. Discretizing the boundary condition (2) gives c _-1,m = c _1,m Substituting these values of c _0,m-1and c _-1,m in (13) gives

The above equation can be written in the form

where

Discretization at ( x ₀ , y _m ), l = 1(1) N , m = - Q _d :- Using the boundary condition (5), the required discretization is

where

Discretization at ( x ₀ , y _m ), m = - Q _d + 1, - Q _d + 2,...., -2, Q _c + 2, Q _c + 3,.., Q _u -1:- Using the the boundary conditions (2) and (4), the required discretization is

where

Discretization at ( x ₀ , y _m ), m = -1:- Using the boundary conditions (2) and (3), the required discretization is

where

Discretization at ( x ₁ , y _m ), m = 0,1,...., Q _c :- Using the boundary condition (3), the required discretization is

where

Discretization at ( x ₀ , y _m ), m = Q _c + 1:- Using the boundary conditions (3) and (4), the required discretization is

where

Discretization at ( x _l , y _m ), m = Q _u :- Using the boundary conditions (2) and (6), the required discretization is

where

Discretization at ( x _l , y _m ), l = 1(1) N , m = Q _u :- Using the boundary condition (6), the required discretization is

where

Discretization along the boundary x = 1, - L _d <= y <= L _u .

At the grid point (x _N+1, y _m ) the differential equation (1) can be written as

Let

c _xx;N+1,m = αc _N,m + βc _N+1,m + γc _xxx;N+1,m     (24)

Using the Taylor-series expansion for C _N,m about the point (x _N+1, y _m ) in (24), and then equating the constant terms, and the coefficients of c _xx;N+1,m and c _xxx;N+1,m , gives

α = 2/H², β = -2(1 + Hw)/H², γ = H/3.     (25)

Differentiating (1) with respect to x and (7) with respect to y and t, respectively, and then using the assumptions c _xyy;N+1,m = c _yyx;N+1,m and c _xt;N+1,m = c _tx;N+1,m , an expression for c _xxx;N+1,m is obtained as

c _xxx;N+1,m = -wc _t;N+1,m + wc _yy;N+1,m     (26)

Substituting (25) and (26) in (24) gives

c _xx;N+1m = αc _N,m + βc _N+1,m + γ{-wc _t;N+1,m + c _yy;N+1,m }     (27)

Substituting (27) in (23) and then replacing c _tN+1,m by and using the approximation in the resulting equation, the required difference scheme along the boundary x = 1 is

where

At an interior grid point (x _l , y _m ), the parabolic equation (1) is discretized as

where

Combining the discretizations (15) to (22), (28) and (29), a system of ordinary differential equations is formed as

where

The value of m in B _d and D _d is -Q _d , and that in B _u and D _u is Q _c + 1. The abbreviations m 1, m 2, ... etc. represent m + 1, m + 2, ... etc. The matrices I _d and I _c are square matrices of order N + 1 and N + 2 respectively and are given by

For m = 0,1,2,...,Q _c , B _m , are square matrices of order N + 1 and each is given by

For m ≠ 0,1,2,...,Q _c , B _m are square matrices of order N + 2 and each is given by

For m = 0,1,...,Q _c , I _m and O _m are the unit square matrix and zero matrix of order N + 1. For m ≠ 0,1,2,...Q _c , I _m and O _m are the unit matrix and zero matrix of order N + 2. Solving the system of ordinary differential equations given by (30), subject to the initial condition(8), gives

C(t) = e^t(B+D) Q- (B + D)^-1 Q;     (31)

which satisfies the recurrence relation

C(t + k) = e^k(B+D) C(t) + (B + D)^-1 e^k(B+D) Q- (B + D)^-1 Q.     (32)

To obtain numerical solutions C(t + k), Padé approximants to the matrix exponential e^k(B+D)may be used in (32). A fully implicit method may be developed by first writing (32) in the form

C(t + k) = e^kDe^kB C(t) + (B + D)^-1 e^kDe^kB Q- (B + D)^-1 Q;     (33)

which has O(k²) error in time. Equation (33) may be approximated further, also with error O(k²), by introducing the (1,0) Padé approximants to the matrix exponential functions e^kD, e^kBand writing (33) in the split form

(I - kB)C* = C(t) + C⁺,     (34)

(I - kD)C(t + k) = C* - (I - kD)C⁺.

where the vector C⁺ is independent of time and is calculated from the equation

(B + D)C⁺ = Q.

This fully implicit algorithm defined by (34) is first-order accurate in time and second order accurate in space, and is L ₀-stable. The accuracy in time is improved by extrapolating, as in Lawson and Morris [7]. Let C⁰, C¹ and C², defined by (35), (36) and (37), respectively, be the three approximations to the solution of (30) at time level t + 2k. The Lawson-Morris algorithm takes the form

(I - kB)C* = C(t) + C⁺,

(I - kD)C(t + k) = C* - (I - kD)C⁺,

(I - kD)C** = C(t + k),

(I - kB)C⁰ (t + 2k) = C** - (I - kB)C⁺.     (35)

(I - 2kB)C* = C(t) + C⁺,

(I - 2kD)C¹ (t + 2k) = C* - (I - 2kD)C⁺.     (36)

(I - 2kD)C* = C(t + 2k) + C⁺,

(I - 2kB)C² (t + 2k) = C* - (I - 2kB)C⁺.     (37)

Then the linear combination of these three schemes (35), (36) and (37) defined by

is second-order accurate in time. That is, the first-order methods (35), (36) and (37) have been extrapolated to achieve second-order accuracy.

The values of the flux J(t _j ) and the cumulative amount of the drug eliminated into the receptor cell per unit area Ae(t _j ), at any time-level t = t _j = jk, j = 0,1,2, ....etc, are calculated using the trapezoidal rule to approximate the integrals in (10) and (11). Thus, for second-order accuracy,

and

Results and discussion

In order to examine the behaviour of the extrapolation method, a number of experiments, similar to those in one space dimension described in [3], are carried out for two space dimensions. The parameter r was given the values 0.00001, 0.0001, 0.001, 0.2, 0.4, 0.6 and 0.8. The space step h was given the value 0.1 so that N = 9 and the time step k is put equal to 0.001. Also, L _c = 0.5, L _u = 3 and L _d = 2.5 so that Q _c = 5, Q _u = 30 and Q _d = 25.

Using the extrapolation technique (38), corresponding to each value of r = 0.001, 0.2, 0.4, 0.6 and 0.8, the concentration profiles at t = 0.1, t = 0.5, t = 1.0 and t = 4.0 are given in Figures 1, 2, 3, 4, 5. These graphs show the effect of the value of the parameter r on the accumulation of the drug into the skin region.

Figure 1

Concentration profiles for r = 0.001

Figure 2

Concentration profiles for r = 0.2

Figure 3

Concentration profiles for r = 0.4

Figure 4

Concentration profiles for r = 0.6

Figure 5

Concentration profiles for r = 0.8

Using the trapezoidal rule given in (39) and (40), the graphs of J versus time and Ae versus time, respectively, are plotted in Figures 6 and 7. These graphs show the effect of the value of r on the the flux and the cumulative amount of drug eliminated into the receptor cell.

Figure 6

J versus time for various values of r

Figure 7

Ae versus time for various values of r

Some authors prefer to retain the dimensions of quantities, but the full range of information can be given in a simple manner when dimensions are eliminated. The authors believe that, for this particular study, it is more appropriate not to use dimensions.

Conclusion

The value of r is not easy to control and the major features relating to it include lipophilicity and the molecular weight of the drug. Various drugs with various lipophilicity and molecular weight will have different r-values, so that the profiles will be different between drugs. The value of r will be just one of the factors in selecting the best drug and best vehicle (ointment) design. Various values of r have been used to illustrate what will be predicted by the model. Based on the profiles given, the following conclusions are drawn.

From Figures 1, 2, 3, 4, 5, the following conclusions are made.

(1) For any value of r, the flow rate towards the skin-capillary boundary is larger than the flow rate near to the skin surface. (2) For larger values of r, the difference between the concentration near to the skin surface and the concentration near to the skin-capillary boundary decreases at a faster rate, compared to smaller values of r.

From Figures 6 and 7, however, the following conclusions are made.

(1) As the value of r increases, the flux and the cumulative amount of drug eliminated into the receptor cell, at a particular time, decrease. (2) As the value of r decreases, the steady-state value of the flux increases. (3) The time required to reach the steady-state value of the flux is smaller for smaller values of r. (4) To the scale used, the graphs for r = 0.00001, 0.0001 and 0.001 were found to be coincident, and it is concluded that the value of r = 0.001 is the optimal value of r which gives the maximum absorption of the drug into the blood stream.

The skin has a complicated structure while the model presented is simple. Though the complicacy itself does not always preclude the successful use of a simple model, if the major assumptions made in this model do not hold, what will be observed may differ from what the model predicts. For example, as to the x-axis, the skin consists of stratum corneum and dermis. The use of a simple model for the x-direction may be warranted because the major barrier is normally assumed to be provided by the stratum corneum, [1].

Considering the diffusion of a drug, however, where the diffusion through dermis is a critical component, the y-directional diffusion of drug molecules through the dermis may also be important and what is observed will be different from what the model in this study predicts. The model presented in this study may therefore be used to examine the mechanism of diffusion in the y-direction.

As shown in the non-linear model in [2], which assumes that the drug molecules are either dissolved (mobile) or immobile inside the skin, the drug permeation profiles are considerably different from those described in the linear model. In this study, it was shown that the drug permeation described by the two-dimensional linear model is different from that by the one-dimensional linear model. The percutaneous drug absorption described by a non-linear two-dimensional model will be studied in future work and further comparison based on [2] will then be warranted.

The mathematical model considered in this paper does not take into account the complex, in vivo physiology of skin. The model is developed by simplifying (or ignoring some complexity of) the skin structure. Several models in the literature, used to describe the kinetic behaviour of a drug, have in fact made this simplification. These models do not give any insight into the drug kinetics at the site distant from the area where the ointment is directly applied. As the title of this paper suggests, the model described has taken into account the above issue in the form of y-directional diffusion. That is, the time-course of the drug concentration (which may be best measured in an in vitro study), at the site distant from the area where the ointment is directly applied, is examined in the present paper. No attempt has been made to relate the model to experimental data, because the authors have no access to experimental data which can be related to the two-dimensional model, where the issue of the y-directional diffusion is examined. The authors hope that publication of this model will stimulate researchers to carry out relevant experimental studies, based on diffusion in the direction parallel to the skin surface.

This two-dimensional model improves upon the much simpler models in one-dimension, by giving information of the time-course of drug kinetics in the skin, where the ointment is not directly applied, assuming that the process can be described by a simple homogeneous membrane. If what is observed is different from what this model describes, this may give some clue to elucidate the mechanism of y-directional diffusion of the drug (for example, the diffusivity in the stratum corneum may be different between x- and y- directions or the y-directional diffusion through the dermis, once it is permeated through the stratum corneum, may change). Those features cannot be known from one-dimensional models.