 Research
 Open Access
 Published:
A twodimensional mathematical model of nonlinear dualsorption of percutaneous drug absorption
BioMedical Engineering OnLine volume 4, Article number: 40 (2005)
Abstract
Background
Certain drugs, for example scopolamine and timolol, show nonlinear kinetic behavior during permeation process. This nonlinear 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 nonmobile 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 dualsorption process gives an explanation to nonlinear kinetic behavior, the characteristic features that are used in onedimensional models are (1) prolongation of halflife if the plot of flux versus time are straight lines soon after the vehicle removal, (2) the decrease in halflife with increase in donor concentration. This paper introduces another feature as a characteristic to confirm that dualsorption model gives an explanation to the nonlinear kinetic behavior of the drug. This new feature is "the prolongation of halflife is not a necessary feature if the plots of drug flux versus time is a nonlinear 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 nonlinear. When a model is nonlinear the relationship between lagtime and donor concentration is nonlinear and the lag time decreases with increase in donor concentration. A twodimensional dualsorption model is developed for percutaneous absorption of a drug, which shows nonlinear 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 onedimensional models. The dualsorption model is an initial/boundary value problem which consists of (1) one nonlinear, twodimensional, secondorder 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 initialvalue problem, which consists of (1) a system of nonlinear ordinary differential equations, (2) one initial condition. The system of nonlinear ordinary differential equations contains timedependent nonhomogeneous terms, if the permeation process under consideration is, during the application of the vehicle. To solve this initialvalue problem, an eightstage sequential algorithm which is secondorder accurate, and requires only tridiagonal 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 xaxis or yaxis. The flow rates in different subregions of the skinregion are studied. The shapes of the concentration profiles are examined before and after the steadystate concentration is reached. The concentration reaches steadystate 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 steadystate value of the concentration, (4) concentration of the drug in different regions of the skin, when steadystate value of the concentration is reached, (5) the time required to reach the steadystate value of the flux, (6) time required to reach the steadystate value of the concentration of the drug, (7) halflife of the concentration of the drug and (8) lagtime.
A comparison, between this twodimensional model and the onedimensional nonlinear dualsorption 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 steadystate value of the concentration, (3) lagtime and (4) halflife.
Background
A nonlinear 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 nonlinear model can be used to explain the disparity between the steadystate diffusivity of the drug in the skin and the unsteadystate value computed from transient (timelag) permeation experiments, see [1]. A nonlinear dualsorption 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 nonlinear 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 nonlinear models are onedimensional 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 twodimensional nonlinear 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.
Methods
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 skinsurface, say, {(0, y) : 0 ≤ y ≤ L_{ c }}, at time t = 0. The twodimensional model developed in [6] takes into account the drug kinetics in the skin where the ointment is not directly applied. The onedimensional nonlinear dualsorption 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 donorcell 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
C_{ T }= C_{ D }+ C_{ I }.
Based on [6] and [2], to describe the drug kinetics at any time t > 0, a twodimensional nonlinear dualsorption model is developed. Let the thickness (distance between the skinsurface and skincapillary 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 drug concentration C_{ D }= C_{ D }(x, y, t) in the skin is governed by
where L_{ s }, L_{ d }, L_{ u }and L_{ c }are positive real numbers. Note that, (1) is nonsingular as b > 0, > 0 and C_{ D }≥ 0. Following [2], the concentration along the skinreceptor cell boundary may be considered to have the value zero. Assuming that the donorcell 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
C_{ D }(0, y, t) = K_{ D }C, 0 ≤ y ≤ L_{ c }, t > 0, (2b)
C_{ D }(L_{ s }, y, t) = 0, L_{ d }≤ y ≤ L_{ u }, t > 0, (2f)
Assuming that there is no drug in the skin before the application, the initial distribution associated with the PDE is
C_{ D }(x, y, 0) = 0,0 ≤ × ≤ L_{ s }, L_{ d }≤ y ≤ L_{ u } (3)
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/boundaryvalue 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 steadystate 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
C_{ D }(L_{ s }, y, t) = 0, L_{ d }≤ y ≤ L_{ u }, t > T_{ s } (7d)
The initial distribution associated with PDE (6) is
C_{ D }= C_{ D }(x, y, T_{ s }), 0 ≤ x ≤ L_{ s }, L_{ d }≤ y ≤ L_{ u }, (8)
where C_{ D }(x, y, T_{ s }) is computed by solving problem (PA). The initial/ boundaryvalue 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/ boundaryvalue problem to be solved is transformed into a firstorder initialvalue problem. At any time t > 0, the region associated with the skin is given by R_{ s } = {(x, y) : 0 ≤ x ≤ L_{ s },  L_{ d }≤ y ≤ L_{ u }}. 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, L_{ u }= Q_{ u }H, L_{ d }= Q_{ d }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 timelevels for the problems (PA) and (PB) for j ≤ s and j ≥ s respectively. The finitedifference solution which approximates the solution C_{ D }(x, y, t) of the twodimensional 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 the initialvalue problem (PB), R_{ c }= Φ (the empty set). For the problem (PA), the differential equation (1), subject to the boundary conditions (2a) to (2f), is discretized at all grid points in the region [R_{ s } R_{ c }{(L_{ s }, y), L_{ d }≤ y ≤ L_{ u }}] × [0 ≤ t ≤ T_{ s }]. 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 [R_{ s }{(L_{ s }, y), L_{ d }≤ y ≤ L_{ u }}] × [t > T_{ s }]. For notational simplicity, denote
Discretization for problem (PA)
At any timelevel t = t_{ j }, or , the {(Q_{ d }+ Q_{ u }+ 1) (N + 1)  (Q_{ c }+ 1)} elements will be ordered in rows parallel to the xaxis and in vector form, and will be denoted by C_{ D }(t) and Q(t), where Q(t) is due to the boundary condition (2b). For W= C_{ D }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_{ D }(t) is a vector of unknowns and the vector Q(t) is given by
All other
are zero. In the following discretizations, (x_{ l }, y_{ m }), p_{l,m}, q_{l,m}and C_{Dl,m}, respectively, denote (x_{ l }, y_{ m }, t_{ j }), , and . At any timelevel t = t_{ j }, the differential equation (1) can be discretized as
If
and are respectively, suitable approximations to C_{Dxx;l;m}and C_{Dyy;l;m}, the above discretization can be written as
At an interior grid point (x_{ l }, y_{ m }), and are given by
If (x_{ l }, y_{ m }) is a boundary point, using the boundary conditions (2a) to (2f), and are given by
where
A_{1} = 0, B_{1} = 2, if l = 0, m = Q_{ d }(1)(1), (Q_{ c }+ 1)(1)Q_{ u },
A_{1} = 0, B_{1} = 1, if l = 1, m = 0(1)Q_{ c },
A_{1} = 1, B_{1} = 1, if l = 1(1)N  1, m = Q_{ d }(1)Q_{ u },
A_{1} = 1, B_{1} = 0, if l = N, m = Q_{ d }(1)Q_{ u },
A_{2} = 0, B_{2} = 2, if l = 0(1)N, m = Q_{ d },
A_{2} = 2, B_{2} = 0, if l = 0(1)N, m = Q_{ u },
A_{2} = 1, B_{2} = 1, if l = 0, m = (Q_{ d }+ 1)(1)(2), (Q_{ c }+ 2)(1)(Q_{ u } 1),
A_{2} = 1, B_{2} = 0, if l = 0, m = 1
A_{2} = 0, B_{2} = 1, if l = 0, m = (Q_{ c }+ 1).
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
where
The value of m in E_{ d }and F_{ d }is Q_{ d }, and that in E_{ u }and F_{ u }is Q_{ c }+ 1. The abbreviations m 1, m 2, ... etc. represent m + 1, m + 2, ... etc. For m ≠ 0, 1, 2, ...Q_{ c }, E_{ m }and F_{ m }are square matrices of order N + 1 and are given by
For m = 0, 1, 2, ..., Q_{ c }, E_{ m }and F_{ m }are square matrices of order N and are given by
For the values of m = 1, Q_{ c }+ 1, and m = 0, Q_{ c }, the matrices G_{ m }are of orders (N + 1) × N and N × (N + 1), respectively, and are given by
For the values of m = 0, 1, ..., Q_{ c }, m ≠ 0, 1, 2, ...Q_{ c }, the matrices O_{ m }are the zero matrices of orders N and N + 1 respectively. The matrices O_{ d }and O_{ c }are zero matrices of orders (N + 1) × N and N × (N + 1), respectively.
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 onespace nonlinear dualsorption model. To the author's knowledge, there is no sequential algorithm in the literature to solve the twodimensional nonlinear parabolic equation (1). The aim of this section is to develop an eightstage sequential algorithm of tridiagonal solvers, to solve the system of nonlinear ordinary differential equations given by (12), which gives the solution of twospace nonlinear 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 twospace nonlinear equations. To achieve this, we proceed as follows.
Let Q(t) = R(t) + S(t), where the elements of the vectors 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 nonlinear 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 premultiplication, gives
Using the split form (15), and following [7] and [2], a fourstage sequential algorithm is developed to obtain , which is an approximation to C_{ D }(t + k) in (17). It is given by
(I  r_{1}kE) Z= C_{ D }(t) + r_{1}k R(t), (18a)
(I  r_{2}kE) V= Z+ r_{2}k R(t), (18b)
(I  r_{1}kF) Z= V+ r_{1}k S(t), (18c)
(I  r_{2}kF) (t + k) = Z + r_{2}k S(t). (18d)
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
(I  r_{1}kF) Z= C_{ D }(t) + r_{1}k S(t), (18e)
(I  r_{2}kF) V= Z+ r_{2}k S(t), (18f)
(I  r_{1}kE) Z= V+ r_{1}k R(t), (18g)
(I  r_{2}kE) (t + k) = Z + r_{2}k R(t). (18h)
Note that the approximations,
and , given by (18d) and (18h), are firstorder accurate in time and their linear combination defined by
is second order accurate in time. The eightstage 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 timelevel t = t_{ j }, j ≥ s, , the (Q_{ d }+ Q_{ u }+ 1) × (N + 1) elements, will be ordered in rows parallel to the xaxis 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
For l ≠ 0, m ≠ (Q_{ d }+ 1) (1) (Q_{ u } 1); and are given by (10) and for l = 0, m = (Q_{ d }+ 1)(1)(Q_{ u } 1),
and
Combining the discretizations given by (19), a system of differential equations is formed as
in which
In the above matrices, the value of m is Q_{ d }and m 1 = m + 1, m 2 = m + 2,... The matrices E_{ i }and F_{ i }, i = Q_{ d }(1)Q_{ u }, are square matrices of order N +1, and are given by (13). Proceeding as in problem (PA), the eightstage algorithm to solve the system of differential equations (20), subject to the initial conditions (8), is
(I  r_{1}kE) Z= C_{ D }(t), (21a)
(I  r_{2}kE) V= Z, (21b)
(I  r_{1}kF) Z= V, (21c)
(I  r_{2}kF) (t + k) = Z, (21d)
(I  r_{1}kF) Z= C_{ D }(t), (21e)
(I  r_{2}kF) V= Z, (21f)
(I  r_{1}kE) Z= V, (21g)
(I  r_{2}kE) (t + k) = Z, (21h)
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 eightstage algorithm developed in problem (PA).
Stability Analysis
The stability analysis of twostage sequential algorithm for the solution of onespace nonlinear secondorder parabolic equation is described in [2]. The stability analysis of fourstage sequential algorithm for the solution of twospace secondorder parabolic equation is described in [7]. In this section, following [7] and [2], it is to prove that the eightstage 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 (18e) to (18h), to evaluate (t + k) is given by
In (22) and (23), z = kλ where λ is an eigen value of the matrix E + F. Hence using (22) and (23) the symbol, S(z), of the method (18i) is
The method is L_{0} stable if
 S(z)  ≤ 1 (24)
S(z) → 0, as z → ∞

Proof of (24): Applying Brauer's theorem, it can be shown that all the eigen values λ of the matrix E + F are negative. Hence, z = kλ > 0.

Proof of (25):
From (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 parametervalues used are those used in [2]. For the experiments numbered n = 1, 2, 3, 4, and 5, the parameter C was given the values 4.1, 19.5, 43.1, 51.4 and 64.0 respectively. The rest of the parameter values are given by
L_{ s }= 0.004 cm, L_{ d }= 0.0320 cm, Lc = 0.128 cm, L_{ u }= 0.160 cm, k = 0.1 h, H = 0.0004, κ = 0.0000018 cm^{2}/h, = 5.0 mg/ml, K_{ D }= 1.1, h = 0.0002, N = 19, Q_{ c }= 320, Q_{ u }= 400, Q_{ d }= 80, b = 0.5599 ml/mg.
It was assumed in all numerical experiments that the drug was applied until steadystate 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 steadystate. 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 steadystate. 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, 28.0, 30.0, 32.0 and 34.0. In figure 3 the profile of concentration C_{ D }, are given at time t = = 128.7, the time at which concentration, C_{ D }, reaches the steady state for the first experiment. In figures 1, 2 and 3, the concentration profiles are viewed as parabolas along mesh lines x = x_{ l }, l = 0(1) N. The figures 1 and 2 show the flow rate in different subregions of the skin region R_{ s }, before the concentration reaches the steady state. Figure 3 shows the shape of the concentration profile when the concentration reaches steadystate.
In figure 4 and figure 5, for the value of C = 4.4, the concentration profiles at t = 0.5 and 3.0 are viewed as parabolas along the mesh lines y = y_{ m }, m = Q_{ d }(1)Q_{ u }. Figure 4 and figure 5 show the change in shape of the parabolas in the region away from the region of contact of drug and skin. In figure 6, corresponding to the value of C = 4.4, the concentration profiles at t = = 128.7 are shown along the mesh lines y = y_{ m }, m = Q_{ c }(1)Q_{ u }. The figure 6 shows the importance of twodimensional modelling.
Corresponding to each value of n = 1(1)5, the concentration profiles at t = are given in figures 7, 8, 9, 10 and 11, respectively. The above figures show the effect of the value of C on the steadystate concentration C_{ D }in different regions of R_{ s }.
In figure 12, the concentration profiles for the value of C = 4.4 are given at t = 128.7, 130.7, 132.7, 134.7, 136.7, and 138.7. This graph show the rate of decreasing of the concentration in different regions of R_{ s }after the removal of the vehicle.
The A_{ e }versus time and J versus time profiles are monitored for the time interval 0 ≤ t ≤ + 50 in each experiment.
The values J and A_{ e }are computed at each time step using the trapezoidal rule to approximate the integrals in (4) and (5). Thus, to secondorder accuracy, for j = 1, 2, 3,... etc,
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 step is used. Note that the conclusions obtained from all the five experiments are independent of the time step k. As an illustration, with k = 0.1 and k = 0.01, J versus time profiles are monitored for the first experiment (C = 4.4) and are given in figure 15. The time required to reach the steadystate value of the flux with k = 0.1 and k = 0.01 respectively are t = 128.7 and t = 130.37. When the numerical computation is done with k = 0.01, the difference in flux during the time interval [128.7 130.37] is 10^{10}, which is negligible. It is also noted that even if the drug is removed at t = 130.37 instead of t = 128.7, the conclusions drawn in next section remain valid.
Conclusion
Conclusion is presented in five parts.
Part 1: Conclusions based on concentration profiles till steadystate 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 y = y_{ m }, m = Q_{ d }(1)Q_{ u }. Divide the region of skin, R_{ s }, into two mutually exclusive regions R_{1} and R_{2} as follows.
R_{1} = R_{ c }× L_{ s }and R_{2} = R_{ s } R_{1}.
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 skincapillary boundary. Lateron time levels (from subplots (c) to (f) of figure 1), flow rate towards the skincapillary 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 skincapillary 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 lateron time levels, flow rate in the region R_{2} becomes larger than the flow rate in the region R_{1}, leading to a steadystate 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 = Q_{ d }(1)Q_{ u }. From figure 4 it is concluded that, during the initial timelevels, 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 onedimensional model of drug absorption was sufficient, the following results are expected from [2].

When steadystate 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 steadystate 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 twodimensional modelling.
5.1 From figure 4 and 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.
5.2 The concentration profiles along the mesh lines, y = y_{ m }, m = 5(1) Q_{ c } 5 are straight lines. These straight lines are shown in figure 6.
5.3 The concentration profiles along the mesh lines, y = y_{ m }, m = 0, 1, 2, 3, 4, Q_{ c } 4, Q_{ c } 3, Q_{ c } 2, Q_{ c } 1 and Q_{ c }are not straight lines, but they have the shape as shown in figure 6.
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 (that is, for all values of C).
Part 2:Conclusions based on concentration profiles when steadystate value of concentration is reached
The steadystate concentration profiles for all the five experiments are given in figure 7, figure 8, figure 9, figure 10 and figure 11. From these steadystate 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_{N1}, 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 steadystate is reached
The vehicle is removed when the concentration reaches steadystate. The concentration profiles are examined for ten more hours, after the vehicleremoval, in an interval of two hours. In figure 12, these concentration profiles are given for n = 1 (C = 4.4), at time levels t = = 128.7, t = 130.7, t = 132.7, t = 134.7, t = 136.7, and 138.7. From figure 12 the following conclusions are made.

(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 skinsurface is larger than that near to region of skincapillary 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.
Part 4:Conclusions based on figure 13(flux versus time) and figure 14(A_{e}versus time)
The graphs of J versus time and Ae versus time for the time interval 0 ≤ t ≤ ( + 50) are shown in figure 13. The following conclusions are drawn based on figure 13.

(1)
As the value of C increases, the steadystate value of flux increases.

(2)
As the value of C increases halflife decreases, during the time levels, soon after the vehicleremoval. There is no prolongation of halflife during lateron time, as stated in [2]. It is evident from figure 13 that, the plots of a drug flux versus time is a nonlinear curve, soon after the vehicle removal. Hence, in experimental studies even when the prolongation of halflife is not seen, the nonlinearity of drug flux versus time can be introduced as a characteristic to confirm that dualsorption model gives an explanation to nonlinear kinetic behavior of the drug.
The graphs of Ae versus time for the time interval 0 ≤ t ≤ + 50 are shown in figure 14.

(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, lagtime decreases. The lagtime is computed as tintercepts of the linear portion of graphs of Ae versus time. For various values of C, the lagtimes are compared with the lagtimes given in [2]. Let T_{L 1}denotes the lagtimes quoted in [2] as a result of onedimensional modelling of drugabsorption and T_{L 2}denotes the lagtimes obtained as a result of twodimensional modelling. Table 1 gives the values of T_{L 1}and T_{L 2}, corresponding to each value of C.
From Table 1 it is concluded that, for a particular value of C, the lagtime obtained as a result of twodimensional modelling is greater than the lagtime obtained as a result of onedimensional 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 lagtime obtained from the onedimensional 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 nonlinear 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 twodimensional nonlinear model is different from onedimensional nonlinear model, as the twodimensional model takes into account the drug kinetics at the site distant from the area where the ointment is directly applied.
References
 1.
Chandrasekaran SK, Michaels AS, Campbell PS, Shaw JE: Scopolamine permeation through human skin in vito. Am Inst Chem Eng J 1976, 22: 828–832.
 2.
Gumel AB, Kubota K, Twizell EH: A sequential algorithm for the nonlinear dualsorption model of percutaneous drug absorption. Mathematical Biosciences 1998, 152: 87–103. 10.1016/S00255564(98)100214
 3.
Kubota K, Ishizaki T: A calculation of percutaneous drug absorptionI. Theoretical. Comput Biol Med 1986, 16: 7–19. 10.1016/00104825(86)900582
 4.
Kubota K, Koyama E, Twizell EH: Dual sorption model for the nonlinear percutaneous permeation kinetics of timolol. J Pharmaceutical Sciences 1993, 82: 1205–1208.
 5.
Kubota K, Twizell EH, Maibach HI: Drug release from a suspension with a finite dissolution rate: Theory and its application to a betamethasone 17valerate patch. J Pharmaceutical Sciences 1994, 83: 1593–1599.
 6.
George K, Kubota K, Twizell EH: A twodimensional mathematical model of percutaneous drug absorption. BioMedical Engineering OnLine 2004, 3: 18. 10.1186/1475925X318
 7.
Gumel AB: Parellel and sequential algorithms for secondorder parabolic equations with applications. PhD thesis. Brunel University; 1993.
 8.
Twizell EH: Numerical Methods, with applications in the Biomedical Sciences. Ellis Horwood, Chichester and Wiley, Newyork; 1998.
Acknowledgements
The author thanks DST, Government of India, for the award of a BOYSCAST fellowship.
Author information
Affiliations
Corresponding author
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
George, K. A twodimensional mathematical model of nonlinear dualsorption of percutaneous drug absorption. BioMed Eng OnLine 4, 40 (2005). https://doi.org/10.1186/1475925X440
Received:
Accepted:
Published:
Keywords
 Concentration Profile
 Timolol
 Time Level
 Sequential Algorithm
 Drug Kinetic