A two-dimensional mathematical model of percutaneous drug absorption

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.

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 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 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 approxi- 0  0  0  0  2   0  10  0  3   , ,  ,  ,  , , mates 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 ,..., ; , ,..., ; , At the grid point (x N+1 , y m ) the differential equation (1) can be written as where the vector C + is independent of time and is calculated from the equation This fully implicit algorithm defined by (34) is first-order accurate in time and second order accurate in space, and is L 0 -stable. The accuracy in time is improved by extrapolating, as in Lawson and Morris [7]. Let C 0 , C 1 and C 2 , defined by (35)

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.
Using the trapezoidal rule given in (39)  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.
(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.
Concentration profiles for r = 0. 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 ydirectional diffusion of drug molecules through the dermis may also be important and what is observed will be differ- ent 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.

Authors' Contributions
All the authors contributed equally to this work.
Concentration profiles for r = 0.8