A two-dimensional mathematical model of non-linear dual-sorption of percutaneous drug absorption

Background Certain drugs, for example scopolamine and timolol, show non-linear kinetic behavior during permeation process. This non-linear kinetic behavior is due to two mechanisms; the first mechanism being a simple dissolution producing mobile and freely diffusible molecules and the second being an adsorption process producing non-mobile molecules that do not participate in the diffusion process. When such 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, C, the donor concentration. The present paper studies the effect of the parameter value, C, when the region of the contact of the skin with drug, is a line segment on the skin surface. To confirm that dual-sorption process gives an explanation to non-linear kinetic behavior, the characteristic features that are used in one-dimensional models are (1) prolongation of half-life if the plot of flux versus time are straight lines soon after the vehicle removal, (2) the decrease in half-life with increase in donor concentration. This paper introduces another feature as a characteristic to confirm that dual-sorption model gives an explanation to the non-linear kinetic behavior of the drug. This new feature is "the prolongation of half-life is not a necessary feature if the plots of drug flux versus time is a non-linear curve, soon after the vehicle removal". Methods From biological point of view, a drug absorption model is said to be nonlinear if the sorption isotherm is non-linear. When a model is non-linear the relationship between lag-time and donor concentration is non-linear and the lag time decreases with increase in donor concentration. A two-dimensional dual-sorption model is developed for percutaneous absorption of a drug, which shows non-linear kinetic behavior in the permeation process. This model 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. The dual-sorption model is an initial/boundary value problem which consists of (1) one non-linear, two-dimensional, second-order parabolic equation, (2) boundary conditions, (3) one initial condition. Note that, the number of boundary conditions are, six and four, respectively, if the permeation process under consideration is, during the application of the vehicle and during the removal of the vehicle. Adopting the approach of method of lines, the initial/boundary value problem is transformed into an initial-value problem, which consists of (1) a system of non-linear ordinary differential equations, (2) one initial condition. The system of non-linear ordinary differential equations contains time-dependent non-homogeneous terms, if the permeation process under consideration is, during the application of the vehicle. To solve this initial-value problem, an eight-stage sequential algorithm which is second-order accurate, and requires only tri-diagonal solvers, is developed. Results Simulation of the numerical methods described is carried out with various values of the parameter C. The illustrations are given in the form of figures. The concentration profiles are viewed as parabolas along the mesh lines parallel to x-axis or y-axis. The flow rates in different subregions of the skin-region are studied. The shapes of the concentration profiles are examined before and after the steady-state concentration is reached. The concentration reaches steady-state when the flux reaches the steady state. The plots of flux versus time and cumulative amount of drug eliminated into the receptor cell versus time are given. Conclusion Based on the various values of the parameter, C, conclusions are drawn about (1) flow rate of the drug in different regions of the skin, (2) shape of the concentration profiles, (3) the time required to reach the steady-state value of the concentration, (4) concentration of the drug in different regions of the skin, when steady-state value of the concentration is reached, (5) the time required to reach the steady-state value of the flux, (6) time required to reach the steady-state value of the concentration of the drug, (7) half-life of the concentration of the drug and (8) lag-time. A comparison, between this two-dimensional model and the one-dimensional non-linear dual-sorption model that exists in the literature, is done based on (1) the shape of the concentration profiles at various time levels, (2) the time required to reach the steady-state value of the concentration, (3) lag-time and (4) half-life.


Background
A non-linear model takes into account the property of the skin to sorb and bind substances during the process of permeation in addition to the ordinary dissolution. Such a non-linear model can be used to explain the disparity between the steady-state diffusivity of the drug in the skin and the unsteady-state value computed from transient (time-lag) permeation experiments, see [1]. 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]. All these non-linear 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 non-linear 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 (1) to model the above situation mathematically (2) to obtain an efficient numerical method which can use a large time step compared to spacial discretization step; and also handles discontinuities between initial and boundary conditions in the mathematical model.

Model equations (a) Concentration before the steady state is reached (Concentration during the application of the vehicle)
Let the drug be applied as an ointment on the skin-surface, say, {(0, y) : 0 ≤ y ≤ L c }, at time t = 0. The two-dimensional model developed in [6] takes into account the drug kinetics in the skin where the ointment is not directly applied. The one-dimensional non-linear dual-sorption model given in [2] postulates the total concentration of the drug, C T , in the skin is composed of two parts, (1) the mobile solute concentration C D , which is due to the mechanism of simple dissolution and is expressed as C D = K D C, in which K D is skin/receptor cell partition coefficient and C is donor-cell concentration, (2) the immobile solute concentration C I , which is due to the adsorption process, and is expressed as , in which is Langmuir's saturation constant and b is Langmuir's affinity constant. Hence the total concentration C T , in the skin is given by Based on [6] and [2], to describe the drug kinetics at any time t > 0, a two-dimensional non-linear dual-sorption model is developed. Let the thickness (distance between the skin-surface and skin-capillary boundary) of the skin be L s . Assuming that skin is an isotropic medium, that is, the diffusivity κ, is the same in the x and y directions, the (1) is non-singular as b > 0, > 0 and C D ≥ 0. Following [2], the concentration along the skin-receptor cell boundary may be considered to have the value zero. Assuming that the donor-cell concentration at the uppermost epidermis {(0, y, t) : 0 ≤ y ≤ L c , t > 0} is maintained at the value C, the boundary conditions for the model are given by 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 The cumulative amount of drug eliminated into the receptor cell per unit area at time τ (see [2]) is The initial/boundary-value problem given by (1) to (3) is named as "problem (PA)".

(b) Concentration after the steady state is reached (Concentration after the vehicle removal)
Let T s denote the time at which the steady-state is reached by the problem (PA). Suppose that at time t = T s , the vehicle is removed instantaneously from the skin. Then (1) remains unchanged, but the initial conditions and boundary conditions changes and consequently, (2a) to (2c) disappear from the mathematical model, leaving the PDE's The initial distribution associated with PDE (6) is (8) where C D (x, y, T s ) is computed by solving problem (PA). The initial/ boundary-value problem given by (6) to (8) is named as "problem (PB)".

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. At any time t > 0, the region associated with the skin is given by The region of the contact R c , of the drug and the skin surface is defined as R c = {(0, y) : 0 ≤ y ≤ L c }. Superimpose on R s a rectangular grid with mesh lengths, h > 0 and H > 0, respectively, in the x and y directions. Let k > 0 be a constant time step. Let N be a positive integer, and h = 1/(N + 1). Also it is assumed that L c = Q c H, where Q c , Q u and Q d are positive integers. 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, ...s, ...etc. with sk = T s . The values t j represent the time-levels for the problems (PA) and (PB) for j ≤ s and j ≥ s respectively. The finite-difference solution which approximates the solution C D (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 -R c -{(L s , y), -L d ≤ y ≤ L u }] × t > 0. Note that, for For the problem (PB), the differential equation (6), subject to the boundary conditions (7a) to (7d), is discretized at all grid points in the region For notational simplicity, denote

Discretization for problem (PA)
At any time-level t = t j , or , the {(Q d + Q u + 1) (N + 1) -(Q c + 1)} elements will be ordered in rows parallel to the x-axis and in vector form, and will be denoted by C D (t) and Q(t), where Q(t) is due to the boundary con- where

System of ordinary differential equations
Combining the discretizations (9) together with expressions for and given by (10) and (11) respectively, a system of ordinary differential equations is formed as . + ( )

Recurrence relation and its implementation via sequential algorithm
In [7], an eight stage sequential algorithm is described to solve two space linear parabolic equations. In [2], a sequential algorithm of two tridiagonal solvers is described to solve one-space non-linear dual-sorption model. To the author's knowledge, there is no sequential algorithm in the literature to solve the two-dimensional non-linear parabolic equation (1). The aim of this section is to develop an eight-stage sequential algorithm of tridiagonal solvers, to solve the system of non-linear ordinary differential equations given by (12), which gives the solution of two-space non-linear parabolic equation (1). The main idea used is, to rewrite the system of equations (12) in such a way that, the numerical techniques described in [7] and [2] can be extended to two-space non-linear equations. To achieve this, we proceed as follows.

R(t) and S(t) are, respectively, denoted by and ,
and are defined by Note that the vector, R(t), is the contribution of C Dxx to the vector Q(t). The vector, S(t), is the contribution of C Dyy to the vector Q(t). The system of non-linear ordinary differential equations given by (12) can be written as It is known (see [2,8]) that the system of ordinary differential equations given by (12), subject to the initial condition (3), satisfies the recurrence relation where D* = diag{d/dt}. The exponential term in the recurrence relation (16) will be approximated by its (2, 0) Padé approximant to give which, following pre-multiplication, gives .
Using the split form (15), and following [7] and [2], a four-stage sequential algorithm is developed to obtain , which is an approximation to C D (t + k) in (17). It is given by where (see [8]) Another solution, , which is an approximation to C D (t + k) in (17), is obtained by interchanging the matrices E and F; and the vectors R and S; in (18a) to (18d). It is given by Note that the approximations, and , given by (18d) and (18h), are first-order accurate in time and their linear combination defined by is second order accurate in time. The eight-stage algorithm, defined by (18a) to (18i) is, L 0 stable and uses only tridiagonal solver.

Initial/boundary value problem (PB)
The "problem (PB)" is modelled by the equations (6) to (8). At any time-level t = t j , j ≥ s, , the (Q d + Q u + 1) × (N + 1) elements, will be ordered in rows parallel to the x-axis and in vector form, and will be denoted as C D (t), with At a grid point x l , y m , the required discretization is Note that the sequential algorithm defined by (21a) to (21i) is obtained from the equations (18a) to (18i), replacing the vectors R(t) and S(t), by the zero vector of length (Q d + Q u + 1) × (N + 1). This shows the efficiency of the splitting of the vector Q(t), and the eight-stage algorithm developed in problem (PA).

Stability Analysis
The stability analysis of two-stage sequential algorithm for the solution of one-space non-linear second-order parabolic equation is described in [2]. The stability analysis of four-stage sequential algorithm for the solution of twospace second-order parabolic equation is described in [7]. In this section, following [7] and [2], it is to prove that the eight-stage sequential algorithms described by (18a) to (18i) and (21a) to (21i) are L 0 stable.
The amplification matrix R*(k(E + F)) of the method defined by (18a) to (18d) is given by neglecting the higher order terms O (-k 3 ).
Hence the symbol R*(-z), of the method defined by (18a) to (18d), to evaluate (t + k) is given by Similarly the symbol R + (--z), of the method defined by  (24) and (25), it is concluded that the numerical method defined by (18a) to (18i) is L 0 stable. Similarly the numerical method defined by (21a) to (21i) is L 0 stable. Since the method is L 0 stable, the discontinuities around y = 0 and y = L c are not propagated.

Results and discussion
In order to examine the behavior of the recurrence relations (18a) to (18i) and (21a) to (21i), a series of five numerical experiments, similar to those described in [2], are carried out for two space dimensions. In these experiments the parameter-values used are those used in [2]. It was assumed in all numerical experiments that the drug was applied until steady-state concentration profile is reached. Let denotes the time in hours, at which the concentration profiles for the n th (n = 1, 2, 3, 4 and 5) experiment reaches the steady-state. Suppose that at time t = , the vehicle is removed instantaneously from the skin. The pattern of the concentration profile is observed for fifty hours more, after the vehicle is removed at the steady-state. For the n th (n = 1(1)5) experiment, the values of C D are computed (1) using the sequential algorithm (18a) to (18i), in the time interval 0 <t ≤ , (2) using the sequential algorithm (21a) to (21i), in the time interval <t ≤ + 50. Following [2], for the value of n = 1 (C = 4.4), the profiles of concentration C D at t = 0.5. 1.0, 2.0, 3.0, 4.5 and 6.0 are given in figure 1. In figure 2, the profiles of concentration C D are given at t = 24.0, 26.0, The values J and A e are computed at each time step using the trape-zoidal rule to approximate the integrals in (4) and (5). Thus, to second-order accuracy, for j = 1, 2, 3,... etc,

R k E F I r kE I r kE I r kF I r kF
where and In literature see [2], the drug absorption experiments are monitored for a large time interval. Hence it is desirable to use a numerical method which gives acceptable results with larger time steps compared to spacial discretization. Since the numerical method used in this section are L 0 stable, it is expected to give accurate results when a large time

Conclusion
Conclusion is presented in five parts.

Part 1:-Conclusions based on concentration profiles till steady-state concentration, C D , is reached
Concentration profiles at various time levels are examined until steady state is reached. Graphically these concentration profiles are represented as parabolas along mesh lines x = x l , l = 0(1) N or parabolas along mesh lines Concentration profiles for C = 4.4, when steady-state is reached Figure 7 Concentration profiles for C = 4.4, when steady-state is reached.
Concentration profiles for C = 19.5, when steady-state is reached Figure 8 Concentration profiles for C = 19.5, when steady-state is reached. Concentration profiles for C = 43.1, when steady-state is reached.
Concentration profiles for C = 51.4, when steady-state is reached Figure 10 Concentration profiles for C = 51.4, when steady-state is reached.
Divide the region of skin, R s , into two mutually exclusive regions R 1 and R 2 as follows.
The following conclusions are drawn based on drug kinetics in the regions R 1 , R 2 and R s .
(1) In figure 1 the concentration profiles are viewed as parabolas along mesh lines x = x l , l = 0(1)N. Consider the flow rate of concentration of the drug in the region R 1 . From the subplots (a) to (c) of figure 1, it is concluded that during the initial time levels, flow rate towards the skin surface is larger than flow rate towards the skin-capillary boundary. Later-on time levels (from subplots (c) to (f) of figure 1), flow rate towards the skin-capillary boundary becomes larger than flow rate near to the skin surface.
(2) In figure 2 the concentration profiles are viewed as parabolas along mesh lines x = x l , l = 0(1)N.
From figure 1 and figure 2 it is concluded that in the region R 2 , the flow rate near to the skin surface is larger than flow rate towards the skin-capillary boundary.
(3) Consider the flow rate in the whole region R s . From figure 1 and figure 2 it is concluded that, during the initial time levels, flow rate in the region R 1 is larger than the flow rate in the region R 2 . But later-on time levels, flow rate in the region R 2 becomes larger than the flow rate in the region R 1 , leading to a steady-state concentration profile as given in figure 3. In figure 3 the concentration profiles are viewed as parabolas along mesh lines x = x l , l = 0(1)N.
(4) At any time level t = t j , consider the concentration of the drug along the mesh line x = x l , l = 0(1)N. From figure  1, 2 and 3, it is concluded that the concentration of the drug at a mesh point (x l , y m , t j ) ∈ R 1 , m = 0(1)Q c is larger than the concentration of the drug at a mesh point (x l , y m , t j ) ∈ R 2 , m ≠ 0(1)Q c .
(5) In figure 4 and figure 5, the concentration profiles are viewed as parabolas along mesh lines y = y m , m = - figure 4 it is concluded that, during the initial time-levels, the shape of the parabolas in the whole region R s retains the same shape as the shape of the concentration profiles given in page 97 of [2]. If a one-dimensional model of drug absorption was sufficient, the following results are expected from [2].
• When steady-state concentration reaches, the shape of the concentration profiles viewed along the mesh lines y = y m , m = -Q d (1)Q u is same as the shape of the parabolas given in page 97 of [2].
• When steady-state concentration reaches, concentration profiles along the mesh lines y = y m , m = 0(1)Q c are straight lines.
Contrary to this result, the following results are obtained from two-dimensional modelling. figure 5 it is concluded that, as time increases, the parabolas in the region R 2 have a convex shape, whereas the parabolas given in page 97 of [2] have a concave shape. Concentration profiles for C = 64.0, when steady-state is reached Figure 11 Concentration profiles for C = 64.0, when steady-state is reached.

Part 2:-Conclusions based on concentration profiles when steady-state value of concentration is reached
The steady-state concentration profiles for all the five experiments are given in figure 7, figure 8, figure 9, figure 10 and figure 11. From these steady-state concentration profiles, the following conclusions are drawn.
(1) As the value of C increases, the concentration at any point (x l , y m , t j ) increases.
(2) As the value of C increases, the difference between the concentrations of the drug at mesh points (x l , y m , t j ) ∈ R 1 , m = 0(1)Q c and (x l , y m , t j ) ∈ R 2 , m ≠ 0(1)Qc increases.
(3) As the value of C increases, the difference between the concentrations of the drug at mesh points (x 0 , y m , t j ) and (x N-1 , y m , t j ) increases for m = -Q d (1)Q c . This increase in difference is more prominent in the region R 1 than in the region R 2 .

Part 3:-Conclusions based on concentration profiles after the removal of the vehicle when steady-state is reached
The vehicle is removed when the concentration reaches steady-state. The concentration profiles are examined for ten more hours, after the vehicle-removal, in an interval of two hours. In figure 12, these concentration profiles are T s 1 (1) After the removal of the vehicle, rate of decrease in the region R 1 is larger than that in R 2 .
(2) After the removal of the vehicle, rate of decrease near to the region of skin-surface is larger than that near to region of skin-capillary boundary. Even though the illustration of figures is done only for one experiment (with n = 1, C = 4.4), that is only for the first experiment, the conclusions drawn are true for all other experiments also. (1) As the value of C increases, the steady-state value of flux increases.
(2) As the value of C increases half-life decreases, during the time levels, soon after the vehicle-removal. There is no prolongation of half-life during later-on time, as stated in [2]. It is evident from figure 13 that, the plots of a drug flux versus time is a non-linear curve, soon after the vehicle removal. Hence, in experimental studies even when the prolongation of half-life is not seen, the nonlinearity of drug flux versus time can be introduced as a characteristic to confirm that dual-sorption model gives an explanation to non-linear kinetic behavior of the drug.    (3) As the value of C increases, at any time level t = t j , the value of Ae(t) increases.
(4) As the value of C increases, lag-time decreases. The lagtime is computed as t-intercepts of the linear portion of graphs of Ae versus time. For various values of C, the lagtimes are compared with the lag-times given in [2]. Let T L1 denotes the lag-times quoted in [2] as a result of onedimensional modelling of drug-absorption and T L2 denotes the lag-times obtained as a result of two-dimensional modelling. Table 1 gives the values of T L1 and T L2 , corresponding to each value of C.
From Table 1 it is concluded that, for a particular value of C, the lag-time obtained as a result of two-dimensional modelling is greater than the lag-time obtained as a result of one-dimensional modelling. For a linear model, lag time is independent of donor concentration C. The linear model corresponding to (1) is obtained by substituting = 0. Its lag time is computed as 1.56 which is greater than 1.44, the lag-time obtained from the one-dimensional linear model (substitute = 0 in (8) of [2]).
It is an established fact that the drug permeation profiles of a nonlinear model is different from that of a linear model, as a non-linear model assumes that the drug molecules are either dissolved or immobile inside the skin. In this article it is shown that the drug permeation described by the two-dimensional non-linear model is different from one-dimensional non-linear model, as the two-dimensional model takes into account the drug kinetics at the site distant from the area where the ointment is directly applied.