A Nonlinear Numerical Model of Percutaneous Drug Absorption K. KUBOTA Division of Clinical Pharmacology, Clinical Research Institute, National Medical Center, Tokyo 162, Japan AND
E. H. TWIZELL Department of Mathematics and Statistics, Brunei University, Uxbridge, Middlesex, England UBS 3PH Received 24 October 1990; revised 26 June 1991
ABSTRACT A nonlinear mathematical model developed by Chandrasekaran et al. is examined to monitor pharmacokinetic profiles in percutaneous drug absorption and is addressed to several associated problems that could occur in the data analysis of in vitro experiments. The formulation of the model gives rise to a nonlinear partial differential equation (PDE) of parabolic type, and a family of finite-difference methods is developed for the numerical solution of the associated initial/boundary-value problem. The value given to a parameter in this family determines the stability properties of the resulting method and whether the solution is obtained explicitly or implicitly. In the case of implicit members of the family it is seen that the solution of the nonlinear PDE is obtained by solving a linear algebraic system, the coefficient matrix of which is tridiagonal. The behaviors of two methods of the family are examined in a series of numerical experiments. Numerical differentiation and integration procedures are combined to monitor the cumulative amount of drug eliminated into the receptor cell per unit area as time increases. It is found that the use of the equation for the simple membrane model to estimate the permeability coefficient and lag time is warranted even if the system should be described by the dual-sorption model, provided cumulative amount versus time data collected for a sufficiently long time are used. However, being different from the behavior in the simple membrane model, the lag time, which can be estimated in this way, is dose-dependent and decreases with increasing donor cell concentration. On the other hand, the permeability coefficient in the dual-sorption model remains constant irrespective of the donor cell concentrations as in the simple membrane model.
1.
INTRODUCTION
Dose or concentration dependency for percutaneous drug absorption not a well-understood concept. Since the early reports of Chandrasekaran MATHEMATICAL
BIOSCIENCES
108:157-178(1992)
OElsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
is et 157
00255564/92/$05.00
158
K. KUBOTA AND E. H. TWIZELL
al. [l, 21 on the dose-dependent kinetics of transdermal scopolamine base, no paper, to our knowledge, has reported definitive evidence of the dependence of the pharmacokinetics on the dose of a drug in a study in vivo or on the concentration in a study in vitro, for percutaneous absorption of any drug other than scopolamine. Recently, Kubota and Yamada [5] reported on the dose-dependent in vitro transdermal kinetics of a finite dose of timololfree base. Although the detailed mechanism of dose-dependent transdermal timolol kinetics is not known, the findings of Kubota and Yamada [5] suggest that such dose dependency in percutaneous absorption kinetics may not be peculiar to scopolamine. Further studies of some of the features unique to dose- or concentration-dependent kinetics are thus warranted. In the present paper, the dual-mode sorption model of Chandrasekaran et al. [I, 21 is examined, with emphasis placed on the in vitro kinetics when the steady state is attained. It is known that J versus time profiles, where J represents flux from skin to perfusing water, cannot be distinguished from those for the simple diffusion membrane when only profiles for one drug concentration in the donor cell are analyzed. This suggests that it is preferable to use several different concentrations of a drug in the donor cell and to evaluate the presence or absence of the dose or concentration dependency of the percutaneous absorption kinetics, particularly those for lag time. The simple skin membrane model should be used for the analysis only after the dose or concentration dependency of the percutaneous absorption kinetics of the drug have been evaluated. It is noted also that the partition and diffusion coefficients are over- and underestimated, respectively, if the simple membrane model is incorrectly employed. The mathematical formulation of the dual-mode sorption model takes the form of a nonlinear initial/boundary-value problem of parabolic type for which the boundary conditions are functions. The derivation of this mathematical model is given in section 2. Later sections contain computational techniques for the numerical solution of the model equation and numerical results obtained using the associated algorithms. The numerical techniques developed are a family of finite-difference methods that obtain the solution of the nonlinear initial/boundary-value problem either explicitly or implicitly by solving a linear algebraic system at each time step. The results obtained are used to generate mobile (dissolved) and immobile solute concentration profiles and to compute the cumulative amount of drug eliminated into the receptor cell per unit area. A computer program written in BASIC and called MULTI (see Yamaoka et al. [14]) is used to estimate the partition parameter of the drug and the associated diffusion parameter. 2.
THE MODEL
EQUATIONS
The dual-sorption model for the in vitro experiment using excised human or animal skin is depicted in Figure 1. In this model, skin is considered to
PERCUTANEOUS
DRUG ABSORPTION
159
Skin
Donor Cell Diffusion Coefficient
D
I
I
Distance
Concentration,
total
dissolved
C
; KmC
x
;
C,
I KDC
I
immobile
Receptor Cell
CD
K I
CI
0
I
cf+aKDc
1 + aK,C I
I
FIG. 1. Dual-mode sorption model for in vitro percutaneous drug absorption. Total drug concentration C, in skin of thickness L is composed of dissolved (C,) and immobile (C,) concentrations. The skin has diffusion coefficient D for the dissolved (mobile) solute, and distance within the skin is represented by x (0 < x < L). At the donor cell-skin boundary (at x = O), C, is given as C&x = 0) = K,C and C, is given as C,(x = 0) = C;aK&/(l + aK,C), where K, is the skin/receptor cell partition coefficient for the dissolved (mobile) solute, C; is Langmuir’s capacity coefficient, a is Langmuir’s affinity constant, and C is the donor cell concentration. At this boundary, C, is given as C,(x = 0) = K,C, where K, is the apparent skin/receptor cell partition coefficient. At the skin-receptor cell boundary (at x = L), CT, C,, and C, are taken to be zero.
be a single membrane with thickness L and when the drug molecules reach the skin-receptor cell boundary they are removed rapidly, so the concentration of drug in the skin at this boundary may be assumed to have the value zero. The two mechanisms, dissolution and adsorption, are postulated by Chandrasekaran et al. [I, 21 as the processes of permeation through the skin; in these mechanisms, mobile molecules that are freely diffusible and nonmobile molecules that do not participate in the diffusion process are produced. Thus, the total concentration, C,, of a drug in the skin can be expressed as the sum of two components, that is to say, c,=c,+c,,
(1)
where C, and C, are the mobile (dissolved) and immobile solute concentrations, respectively (see also Paul [lo]). The value for C, on the uppermost epidermis, Co( x = 0) is given by CD( x = 0) = K,C,
(2)
160
K. KUBOTA AND E. H. TWIZELL
where C is the drug concentration in the solution in the donor cell and K, is a constant. In the classic steady-state percutaneous absorption experiment (Scheuplein and Blank [12]), an infinite supply of drug is effectively allowed for in the donor cell. Thus, C in Equation (2) can be regarded as a constant. It is noted that K, is the skin/vehicle partition coefficient for the mobile solute concentration and is equal to the effective skin/vehicle partition coefficient, K, [12] only when Co * C,. On the other hand, C, is given as c,=c;ac,/(l+ac,),
(3)
where C; is Langmuir’s saturation constant and II is the affinity constant. As the mobile solute concentration CD becomes large, C, tends to CJ, though it must be remembered that Co cannot exceed the solubility of the drug in skin; and as the affinity constant (I becomes large, C, tends to C; more rapidly with increasing Co. The flux J through the skin per unit area is given by the formula J=-Dax,
ac,
where D is the diffusion coefficient for the mobile solute through skin. The cumulative amount Ae eliminated into the receptor cell per unit area is given by the equation
'ac, x(x=
Ae=Ae(f)=j”J(x=L)dt=-D 0
s
and Fick’s second law of diffusion
for the dual-sorption
ac,
atzD
qdt
(5)
0
model is given as
a%, ax2
’
where
ac, ac, -=,+g[1+ at It follows from (6) and (7), therefore,
(I;j;D)2]%. that
o-0 (see also Frisch [4] and Paul [lo]).
PERCUTANEOUS
DRUG ABSORPTION
161
Chandrasekaran et al. [I] reported the values K, = 1.1, CT = 5 .O mg/mL, a = 0.509 mL/mg, and D = 5.0 x lo- lo cm2/s = O.OOOOOl8 cm2/h, while Michaels et al. [8] give L = 40pm = 0.004cm. Chandrasekaran et al. [1] also report the following set of values of C: (4.4,
19.5,
43.1,
the elements of which are measured
51.4,
in mg/mL
64.01,
(9)
and may be used in Equation
(2). Employing these values of K,, Cy, a, D, and L, the initial/boundaryvalue problem to be solved consists of the partial differential equation (PDE)
[It,1;;;DJ2]%= Ds;
0< x
O,
(10)
together with the initial distribution
c,( x,0) = 0; and the boundary
Oo t>o.
(12) (13)
takes a value from (9). The problem for C, using the family of numerical corresponding value of Cr. can then (I), which give
c,-c,[*+&]. 3.
ANALYSIS
USING SIMPLE
MEMBRANE
MODEL
One of the major interests in this paper is to elucidate the consequence of a case that is likely to occur in actual in vitro experiments. In such a case, the cumulative amount excreted into the receptor cell, which should be explained by the dual-sorption model, is inadvertently analyzed using the simple membrane model. The steady-state pharmacokinetics, observed when the donor cell concentration C is unchanged, can be characterized by (i) the steady-state flux,
162
K. KUBOTA AND E. H. TWIZELL
and (ii) the lag time TL, defined as the value for the t intersect of the asymptote of the Ae versus time curve. Of these two values, TL but not JS, reveals the nonlinear dependence on the donor cell concentration. This is because in the dual-sorption model the concentration gradient at the steady state of CD, dC,/ax, becomes a constant throughout the skin as in the simple membrane model and is given as - Co( x = O)/ L. J in (4) now becomes a constant, J,,, which is given by J,,,
J,, = DC,(x
= 0)/L
= DK,C/L.
(15)
Equation (15) designates that the steady-state flux Js, is proportional to the donor cell concentration C in the dual-sorption model, as in the simple membrane model. The permeability constant P, defined as
P= J,,/C, is thus simply described
as a constant P=
(16)
given by
DK,/L=
k”,L;,
(17)
where k> and L*d are the diffusion parameter and partition parameter defined by Okamoto et al. [9] in the dual-sorption model and given, respectively, as k; = D/L*,
(18)
L; = K,L.
(19)
On the other hand, according to the explicit expression for the lag time TL obtained by a general method for determining asymptotic steady-state solutions to Fick’s second law for permeation [lo], the TL in the dualsorption model is known to be dependent on the donor cell concentration [l, 2, lo] as l/2( OK&)*+
UK&-
(I+ UK&)
( aK&)3
m(l+
aK,C) I
Therefore, if the Ae-time data for more than two different values for the donor concentration C are available, it may not be difficult to see that the nonlinear model is required to analyze the data. However, it has not been determined whether any discrepancy occurs between the observed data and what the model predicts if the data of Ae
PERCUTANEOUS
163
DRUG ABSORPTION
versus time for one donor cell concentration are analyzed by a simple membrane model while the dual-sorption process prevails. In this paper, this problem is examined as follows; several Ae versus time plots are generated assuming that the dual-sorption process prevails in the in vitro diffusion cell experiment, and those data are fitted to the equation of the linear model given below by (see e.g., Crank [3, pp. 42-611 and Okamoto et al. [9])
r=1
kd,t-k--$e
(-k
(-l)‘exp
r%r*t) rz
1*(21)
Another question that has not been elucidated so far is whether the analysis of the Ae versus time data using the simple linear membrane model gives the correct J,, and TL values if the dual-sorption mechanism prevails. This might be a practical problem, particularly when enough of the Ae-time data after the steady state is reached are not available. In such a case, the Ae-time data are sometimes analyzed [9] using (21); that is, in general situations, when Ae is known or has been determined for a discrete set of values {t,, I,, . . . , t,}, the associated partition parameter L,, and diffusion parameter kdl can be estimated using formula (21) and the adaptation of the Levenberg-Marquardt algorithm given in [6], [7], and [ 131. The values TL,, P, , and Jss, are given using the estimated parameter values L,, and k,, in (21) as T,, = 1/6&j,, P, =
(22)
kd,Ld,3
(23)
and J,,, = P,C.
4.
NUMERICAL
SOLUTION
OF THE NONLINEAR
(24)
PDE
Consider the initial/boundary-value problem consisting of the nonlinear parabolic partial differential equation (lo), the initial distribution (1 I), and the boundary conditions (12) and (13). It will be assumed that the parameters C;, a, and K, have been determined and that C takes a discrete value from a set such as (9). Suppose that the interval 0 < x < L is divided into M + 1 subintervals each of width h so that (M + 1) h = L and the time variable t is descretized in steps of length 1. The open region R = [0 < x < L] x [t > 0] and its boundary dR consisting of the lines x = 0, x = L, and t = 0 have thus been covered by a rectangular mesh, the mesh having coordinates x, = mh
164
K. KUBOTA AND E. H. TWIZELL
and t,=nl,wherem=O,l,..., M+landn=0,1,2 ,.... Thetheoretical solution of the initial/boundary value problem {(lo), (1 l), (12), (13)) at the point ( x,, t,) is CD( x,, t,); the solution of an approximating numerical method at the same point will be denoted by 7;. It will be convenient to let aCJ*]-’ [see Eq. (lo)] and then p = DI/h* and q(C,) = [l+ C;a/(l+ and n = to let q~=[l+C;a/(l+ay~)*]-’ with m=1,2,...,M 0,1,2 )... . The derivatives in the PDE (10) will be approximated by the finitedifference replacements
and a2C,(W)
a2
= h-*[@{C,(x-
h,t+l)-2C,(x,t+l)
+C,(x+h,t+l)) +C&+
h,‘)}]
+(1-B){C,(x-h,t)-X&t) + C(h*),
(26)
in which 13 is a parameter 0 < 0 < 1. Evaluating q(C,) point (x,, t,) leads to the family of numerical methods -
pq;ey;::
+ (1+2pq;e)7:+’
-
in (10) at the base
p4;wz:
=pq~(l-e)y~_,+[1-2~q~(l-e)]~~+~q~(l-e)~~+,. The principal benefit gained by evaluating q(C,) at the point tion of any of the methods in (27) and resisting the temptation to at the advanced time level t = t,+ , is that the solution of nonlinear problem can be computed by solving a linear algebraic each time level. This linear algebraic system has the form A(“)r”+l
= B’“‘r”
+,$I),
(27) of applicaevaluate it the given0 system at
(28)
in which A(“) and B(“) are tridiagonal matrices of order M and the vector CC”)= [c’,“‘, 0, . . . ,OIT is also of order M with c\“) = pqyK,C; the boundary condition (12) is responsible for the absence of r3 from cr. The matrices A’“’ and B(“) are given in the Appendix. A fast, direct, well-known method for calculating P ’ + ’ from (28) is described in Twizell [13], [Ch. 11. The local truncation errors associated with the family of methods given
PERCUTANEOUS
DRUG ABSORPTION
16.5
by (27) are 0(/r* + 1) for 0 # l/2 and 0(/z* + 1*) for 8 = l/2. The von Neumann necessary condition for stability requires the localized restriction
where the subscript S refers to the spectral norm (see Twizell [ 13, Chs. 7 et seq.]). Clearly, any stability restriction will require modification at each time step because the PDE is nonlinear. The choice of which value of 0 to choose is influenced by a number of factors. Putting 8 = 0 allows the solution to be obtained explicitly because then A(“)=1 (the identity matrix of order M) in (28), but a linear stability analysis imposes the restriction pqL’ < l/2 in integrating from time t = t, to time t = t,,,, where qi’ = max,( qi 1. Clearly, this method of the family (27) is competitive only when the user can afford to use a small time step. Putting 8 = 1 gives a fully implicit method [B’“‘=I in (28)] that will be especially suitable for solutions that decay rapidly, and putting 0 = l/2 gives a more accurate implicit method that may not be suitable for use when there are discontinuities between initial conditions and boundary conditions, as indeed there are with the present problem [see (11) and (12)] if 1 is too large with respect to h (see Twizell [13, Ch. 81). Alternative methods of obtaining a numerical solution include the use of the so-called method of lines in which the parabolic initial/boundary-value problem ((lo), (1 l), (12), (13)) is transformed into a first-order initial-value problem. This approach is particularly useful when the PDE is linear and is described for a linear model of percutaneous drug absorption in Twizell [ 13, Ch. 81. In the case of a nonlinear PDE, the first-order initial-value problem must be integrated by a numerical method designed to solve such problems. These methods include Runge-Kutta methods, linear multistep methods employed in predictor-corrector combinations, hybrid methods, and multiderivative methods. Stability restrictions depend on the method being used and on the PDE itself and can be determined using the theory of weak stability for ordinary differential equations (see Twizell [ 13, Ch. 51). 5.
NUMERICAL
RESULTS
Using the values of K,, C;, a, D, and L reported in Section 2, the initial/boundary-value problem to be considered consists of the partial differential equation l+
2;;;;CDJZ
acLJ ac, 1 at=0.0000018
at
,
Oo,
= 0;
where C is obtained from (9). In order to examine the behavior of the family of computational methods given by (27), a number of numerical experiments were carried out. The parameter C was given the value 4.4 mg/mL first of all, the space step h was given the value 0.0002 cm, so that M = 19, the time step was put equal to 0.1 h, and the parameter 8 was given the values 1 and l/2. The solution y (an approximation to CD) of {(30), (31), (32), (33)) was then calculated for both values of 0 until the steady state was reached. It was found that the steady-state solution, the straight line joining the points (O,K,C) = (0, l.lC) and (L,O) = (0.004, 0) in the (x, Co) plane, was reached after about 170 time steps. The output of the computer program gave the values of y = Co, where x = 0, 0.0002, 0.0004,. . . ,0.004 cm, the corresponding values of C, being computed using (14). The numerical results for Co and C, are depicted in Figures 2 and 3, respectively, at times 0.5 h, 1.0 h, 2.0 h, 3.0 h, 4.5 h, 6 h, and 15 h (that is, after 5, 10, 20, 30, 45, 60, and 150 time steps). To the scales used in Figures 2 and 3 the corresponding curves for 8 = l/2 at the higher values of t are superimposed. An example of the closeness of the results for Co, for 0 = l/2 and 0 = 1, at time t = 4.5 h is revealed in Table 1. Using the fully implicit method (0 = l), the solution is seen to decay more quickly as x increases than when the implicit method (0 = l/2) is used, as was noted earlier. At the lower values of t the differences between the values of Co were more marked, as can be seen in Figure 4, where the distributions for f3 = l/2 and 0 = 1 at time t = 0.5 h (after five time steps) are depicted. The theoretical solution of the initial/boundary-value problem is not known, so it is not possible to say which of the two values of 8 gives the “better” set of results. Some indication may be obtained by refining the mesh. For instance, with 0 = 1, reducing the space step h by a factor of 3 and the time step 1 by a factor of 10 reduces the local truncation error by a factor of about 10. However, the storage and CPU times will be increased substantially by such moves, and the user may well be satisfied with results obtained on a coarser mesh.
PERCUTANEOUS
167
DRUG ABSORPTION
CD(mg/ml)
0
8
24
16
32
40
x(w) FIG.
2. Values of C,
using 9 = 1, I = 0.1 h, h = 2 pm, M = 19, and C = 4.4
mg/mL.
The numerical experiments were then repeated for the other values of C listed in (9). There was little appreciable difference between the distributions computed using any of the five values of C, with either 8 = l/2 or 1, at any time t. Once again the differences were greater at the lower values of t, although only those for the least and greatest values of C (4.4 and 64.0) were markedly different. These distributions of CD at time t = 0.5 h (after five time steps), normalized as Co /K,C, are depicted in Figure 5. The differences between the distributions CD are due to the nonlinear dose response of the model; when C = 64.0 mg/mL, the distribution precedes that of C = 4.4 mg/mL because the small fraction of drug that has diffused
168
K. KUBOTA AND E. H. TWIZELL CThg/ml)
9
1
0
8
24
16
32
40
XC!Jd FIG.
3. Values of C,
using 0 = 1, I = 0.1 h, h = 2 wm, M= 19, and C = 4.4
mg/mL.
into the deeper portion is immobilized when C = 64.0 mg/mL as shown in (1) and (3). This difference would be reflected eventually by the dosedependent behavior of the lag time in (20). The value of Ae at time t can be estimated from Equation (5) by employing a numerical integration formula. First of all, however, X, /ax must be estimated at x = L using a numerical differentiation formula. The numerical method used to solve the model PDE (10) incurred an 0( h*) error in approximating the space derivative by (26). There is no point, therefore, in using a numerical differentiation formula that enjoys better than second-order accuracy to estimate 0” = NJ,< L, nl)/ax. To this end, the second-order backward difference formula P”=(r~-,-4y~++3y~+,)/2h+O(h2) can be used (see Twizell [13; p. 921); note yl;fr+, = 0 from (13).
(34)
PERCUTANEOUS
169
DRUG ABSORPTION TABLE 1
Values of C, (mg/mL) at Time t = 4.5 h Using (27) with 19= l/2 and 1 and C = 4.4 mg/mL
e=i
e=1/2
x(e) 0
4.840 4.285 3.735 3.195 2.671 2.166 1.685 1.229 0.799 0.393 0.000
4 8 12 16 20 24 28 32 36 40
4.840 4.280 3.726 3.183 2.656 2.150 1.669 1.215 0.789 0.387 0.000
It is an easy matter to adapt the computer program to compute 6C,(L,n1)/6x at each time step, following which Ae can be computed using the trapezoidal rule Ae(nl)=-+Dl
/3°+2nc1/3q+pn I
(see
Twizell
6.
RESULTS METHOD
q=l
[13, p. 991). OF THE
ANALYSIS
BY THE
1
+0(1*)
SIMPLE
(35)
MEMBRANE
The Ae versus time plots produced by (35) for C = 4.4 (triangles) and 64.0 (circles) mg/mL, taking h = 0.0002 cm and I= 0.1 h, are shown in Figure 6a. The solid lines in Figure 6a are the Ae(t)-time curves given by (21) fitted to the hourly Ae versus time plots (up to 24 h) generated by (35). The nonlinear least squares computer program written in BASIC and called MULTI (see Yamaoka et al. [14]) is used. The plots generated by Equation (35) are well explained by the “incorrect” model equation (2 1). To illuminate that Js, and P are dependent not on C but on lag time TL (20), or that T,., (22) is dependent on the donor cell concentration, Ae(t)/C versus time plots are shown in Figure 6b. In Table 2, the values of k,, and L,, estimated by fitting the hourly Ae(t) versus time until t, (t, = 4, 6, 8, 10, 12, 18, and 24 h), generated by (35), to Equation (21) are shown. The TL, (22), P, (23), and J,,, (24) values are also shown in Table 2 and compared to TL (20), P (16), and J,, (15) in the dual-sorption model, respectively When C = 4.4 mg/mL, the numerical values for TL, P and Js, were 3.33 h, 0.495~ 1O-3 cm/h, and
K. KUBOTA AND E. H. TWIZELL CD(mg/ml)
0
8
24
16
32
40
*(urn)
FIG. 4. Values of C, after five time steps using I = 0.1 h, h = 2 pm, M = 19, and C = 4.4 mg/mL.
2.18 pg/(cm’*h), respectively, and when C = 64.0 mg/mL those values were 1.81 h, 0.495 X 10e3 cm/h, and 31.7 pg/(cm2*h), respectively. As is shown in Table 2, the values for P,, JSS,, and TL, (particularly that for TL,) particularly when are overestimated and the k,, value is underestimated, the Ae versus time plots collected for the short-time duration containing no or few data after the steady state is attained are fitted to (21). However, when the Ae(t) versus time plots collected for the long time duration are used, those values for TL, (22), P, (23), and J,,, (24) approach TL (20), P(16), and J,, (15), respectively, even if the “incorrect” model equation (21) is used. For instance, in Table 2, the error of the estimated TL, value (22) is within 20% of the true TL value (22) when the Ae versus time plots
PERCUTANEOUS
171
DRUG ABSORPTION
c,/ CKp 1.0
0.8
-
0.6
-
0.4
-
0.2
-
0
I
I
8
16
I
24
32
40
x (urn)
FIG. 5. Values of C, (normalized to C,/K,C) at time t = 0.5 h using 0 = 1, I = 0.1 h, h = 2 pm, and M = 19. Values of C are in mg/mL.
collected for more than twice the estimated TL, value are analyzed. When the Ae versus time plots collected for more than four times the estimated TL, value are analyzed, the error is within 10% of the true value. 7.
DISCUSSION
The primary conceptual parameter value estimated in steady-state percutaneous drug absorption kinetics is the permeability constant P defined by Scheuplein [II] and Scheuplein and Blank [12] as in (16). Based on the simple membrane model, there are two alternative approaches for estimating diffusion and partition parameters. In the first (approach A), the diffusion parameter is estimated from the lag time T, [I].
172
K. KUBOTA AND E. H. TWIZELL 800
cy -
600
E 3J 2
400
$ 200
0
6
12 Time
18
24
(h)
@) FIG. 6. (a) Profiles of cumulative amount excreted into receptor cell, Ae, plotted versus time, generated in the dual-sorption model where donor cell concentrations C are 4.4 (A) and 64 (0) mg/mL. The solid lines represent the best fits predicted by Equation (21) for the simple skin membrane model. (b) Ae/C versus time profiles, symbols and solid lines represent the same profiles as in (a).
As given concentration
in (20), the lag time TL is dependent C. Defining the lag time prolongation
on the donor cell factor S to be the
quantity
S=l+6C~al/2(aK,C)2+aK,C-(1+aK,C)In(l+aK,C) I ( aKLG3
7 (36)
PERCUTANEOUS
DRUG ABSORPTION
173
TABLE 2 Values of k,, , L,,, P,, J,, and TLl Estimated When Hourly Ae- t Plots up to Time f,(h) Generated by the Dual-Sorption Model Are Fitted to Equation (2 1) Describing the Simple Diffusion Model.
tz (h)
Qa @-‘I
Plb
&I,=
JSS,'
(cm x 10-3)
(cm/h x 10e3)
C = 4.4 mg/mL 4 0.0356 6 0.0395 8 0.0422 10 0.0441 12 0.0454 18 0.0474 24 0.0482 _e 0.0507f
20.7 15.4 13.1 12.0 11.4 10.6 10.3 9.179
0.736 0.606 0.553 0.530 0.517 0.503 0.499 0.495h
C = 64 mg/mL 4 0.0781 6 0.0843 8 0.0871 10 0.0885 12 0.0894 18 0.0906 24 0.0911 _e 0.0953f
7.14 6.13 5.81 5.66 5.58 5.48 5.44 5.20g
0.557 0.517 0.505 0.501 0.499 0.497 0.496 0.495h
[~gI(cm*.h)l 3.24 2.67 2.44 2.33 2.28 2.21 2.20 2.18 35.7 33.1 32.3 32.1 31.9 31.8 31.7 31.7
TLld (b) 4.68 4.22 3.95 3.78 3.67 3.52 3.46 3.29’ 2.13 1.98 1.91 1.88 1.87 1.84 1.83 1.75’
a Estimated from (21).bCalculated as k,,L,, (23). ‘Calculated as P,C (24). d Calculated as 1/6k,, (22). eIdeal case where P, (23) and TL, (22) are estimated exactly as P (17) and T, (20) in the dual-sorption model, respectively. f Estimated as 1/6T,. ‘Estimated as P/k,. hGiven in (17). i Given in (20).
a new conceptual
parameter,
TL, defined as
T;= T,/S=
L2/6D=
1/6k;,
may be introduced. Because S approaches unity as C approaches infinity, Ti is the fictive lag time obtained when the donor cell concentration is assumed to be infinite. In this section, a mechanism is outlined for estimating the parameter values, using the simple membrane model, from the data where the dualsorption process prevails. First, though the lag time TL is dose-dependent, the profiles of Ae versus time in Figure 6, observed for dual-mode sorption models, are difficult to distinguish from those for simple membrane models whenever
174
K. KUBOTA
AND E. H. TWIZELL
results for only one concentration of a drug in the donor cell are analyzed. It is thus quite likely in such cases that the diffusion and partition coefficients will be estimated using a simple membrane model even when a dosedependent model is required. The apparent diffusion parameter k,,, obtained when the simple membrane model is appropriate, is given, using S defined in (36), as
k,,= k’,lS, following equation
which the apparent partition
parameter
(38) L,, can be defined by the
L,, = L*dS. Second, the results shown in Table 2 indicate that the values for TL, (22), P, (23), and J,,, (24) of the dual-sorption model are well tracked as the values in the dual-sorption model, T, (20), P (17), and Js, (15)) respectively, even when the incorrect model equation (21) is employed, provided that Ae versus time plots collected for a sufficiently long time are analyzed. Thus, it is warranted to use the ‘ ‘ incorrect’ ’ model equation (2 1) to estimate the TL (20), P(17), and Js, (15) values in the dual-sorption model from Ae versus time data obtained for each donor cell concentration. This feature is of practical significance because TL, P, and Js, are sometimes estimated by fitting the Ae versus time data to Equation (21) (see Okamoto et al. [9]) even if the data after the steady state is attained are available. This procedure is useful because the actual Jss values in the experiments are usually fluctuating, and in some cases it is difficult to find the time after which J,, can be considered to have been attained. However, it is notable that if few or no Ae-time data after the steady state is attained are available, Ld,, P,, TL, and J,,, are generally overestimated and k,, is underestimated (Table 2). Third, if the TL, /S values [which should be TL in (37)] for several donor concentrations are of a similar magnitude to one another, it suggests that the dual-sorption process is the likely mechanism causing the nonlinear pharmacokinetic profiles. This is because the values of the lag-time prolongation factor S in (36) are determined by a, C;, and K, obtained from the sorption isotherm experiment, and the consistency of TL /S values for different donor cell concentrations indicates that two nonlinear features (observed in the sorption isotherm and diffusion experiments) are explained as two aspects of a dual-sorption model. In the second method (approach B) to estimate the diffusion and partition parameters based on the single membrane model, the partition parameter is estimated from the sorption isotherm experiment as reported by Scheuplein
PERCUTANEOUS
DRUG ABSORPTION
175
[ 1 I 1. In situations where a dual-sorption model is appropriate, when equilibrium is attained in the sorption isotherm experiment, the total concentration of drug, C,, is constant throughout the skin (the epidermal sheet) and is given by the equation
(40) where C is the concentration Defining R by
of the drug in the solution
R=l+Cfa/(l+aK,C), the apparent partition
coefficient K,
The corresponding
partition
K,
and defining
kd2 as k,
parameter
= P/L,?
(41)
is given as
= CT/C=
Ld, = K,L
at equilibrium.
(42)
K,R. Ld2 is given by
= K,RL
= L*,R
(43)
it follows that
kd2 =
k”,lR,
(44)
where L>, ki, and P are given in (19), (18), and (17), respectively. Though the above expressions for Ld2 and kd2 are dependent on the concentration of solution at equilibrium, C, it is quite likely that only one concentration at equilibrium will be employed to estimate the partition coefficient K,, if it is not realized that the dual-sorption process occurs in the experiment. As is seen in Table 3, when the parameter values for scopolamine are employed (see (9) and [2]), the values of R and S are of the same order but are not identical. Thus the parameter values estimated by
TABLE 3 Values of R and S Corresponding to C in (9) [2] c
R
s
4.4 19.5 43.1 51.4 64.0
1.735 1.214 1.101 1.085 1.069
2.221 1.481 1.255 1.219 1.181
176
K. KUBOTA AND E. H. TWIZELL
approaches A and B would be of similar magnitude to each other if the concentration employed at equilibrium is of the same magnitude as that in the diffusion experiment. Thus, in this case, the fact that the dual-sorption process prevails may again be easily missed. In conclusion, it is important to employ various donor cell concentrations in the diffusion experiment and various concentrations at equilibrium in the sorption isotherm experiment to know whether the dual-sorption process occurs in the experiment. If a single concentration is employed in the diffusion and sorption isotherm experiments, an important feature, that the process should be described by the dual-sorption model but not the simple diffusion model, might not easily be recognized. An early priority is the ability to estimate the parameters CT and a, first used in Equation (3). Chandrasekaran et al. [l, 21 reported values of C, for the set of values of C given in (9). The parameters Ci and a may be estimated using the data pairs (C, Cr.) with Equation (40); the Levenberg-Marquardt algorithm for doing so is given in Levenberg [6], Marquardt [7], and Twizell [ 13, Ch. 21. The research for this paper was begun while one of us (E. H. T.) was a visitor to the National Medical Center (NMC), Tokyo. The substantial grant received from The Royal Society of London that made the visit possible, and the help of Dr. Takashi Ishizaki (NMC) in organizing the visit, are gratefully acknowledged. We are very grateful to two anonymous referees for their many valuable comments and to Mrs. B. Yates and Ms. M. Demmar of Brunei University for their help in preparing the typescript. APPENDIX.
THE MATRICES
A’“’ AND B(“)
The matrices A’“’ and B(“) are given by
where a:=
1+2pqG8,
b$=
-pq;e,
PERCUTANEOUS
DRUG ABSORPTION
177
and
0,” B(n)
0
Pz”
= 9
0
KG
44-I Kl
I
PZG1 44
with
a;=1-2pq;(l-e),
The meanings
of p, q;
P;=Pq;(l-q.
and ~9are described
in the main body of the paper.
REFERENCES 1 S. K. Chandrasekaran, A. S. Michaels, P. S. Campbell, and J. E. Shaw, Scopolamine permeation through human skin in vitro, Am. Inst. Chem. Eng. J. 21:828-832 (1976). 2 S. K. Chandrasekaran, W. Bayne, and J. E. Shaw, Pharmacokinetics of drug permeation through human skin, J. Pharm. Sci. 67: 1370- 1374 (1978). 3 J. Crank, The Mathematics of Diffusion. 2nd ed., Clarendon Press, Oxford, 1975. 4 H. L. Frisch, The time lag in diffusion, J. Phys. Chem. 61:93-95 (1957). 5 K. Kubota and T. Yamada, Finite dose percutaneous drug-absorption: theory and its application to an in vitro timolol permeation, J. Pharm. Sci. 79:1015-1019 (1990). 6 K. Levenberg, A method for the solution of certain nonlinear problems in least squares, Q. Appl. Math. 2:164-168 (1944). 7 D. W. Marquardt, An algorithm for least squares estimation of nonlinear parameters, SZAMJ. 11:431-441 (1963). 8 A. S. Michaels, S. K. Chandrasekaran, and J. E. Shaw, Drug permeation through human skin: theory and in vitro experimental measurement, Am. Inst. Chem. Eng. J. 21:985-996 (1975). 9 H. Okamoto, M. Hashida, and H. Sezaki, Structure-activity relationship of I-Alkylor 1-alkenylazacycloalkanone derivatives as percutaneous penetration enhancers, J. Pharm. Sci. 77~418-424 (1988). 10 D. R. Paul, Effect of immobilizing absorption on the diffusion time lag, J. Polym. Sci. (Part A-2) 7:1811-1818 (1969). 11 R. J. Scheuplein, Mechanism of percutaneous absorption. I. Routes of penetration and the influence of solubility, J. Invest. Dermatol. 45:334-346 (1965).
178 12
K. KUBOTA AND E. H. TWIZELL
R. J. Scheuplein and I. H. Blank, Permeability of the skin, Physiol. Rev. 51:702-747 (1971). 13 E. H. Twizell, Numerical Methods, with Applications in the Biomedical Sciences, Ellis Horwood, Chichester and Wiley, New York, 1988. 14 K. Yamaoka, Y. Tanigawara, T. Nakagawa, and T. Uno, A pharmacokinetic analysis program (MULTI) for a microcomputer, J. Pharmacobio-Dyn. 4:879-885 (1981).
E. H. TWIZELL Department of Mathematics and Statistics, Brunei University, Uxbridge, Middlesex, England UBS 3PH Received 24 October 1990; revised 26 June 1991
ABSTRACT A nonlinear mathematical model developed by Chandrasekaran et al. is examined to monitor pharmacokinetic profiles in percutaneous drug absorption and is addressed to several associated problems that could occur in the data analysis of in vitro experiments. The formulation of the model gives rise to a nonlinear partial differential equation (PDE) of parabolic type, and a family of finite-difference methods is developed for the numerical solution of the associated initial/boundary-value problem. The value given to a parameter in this family determines the stability properties of the resulting method and whether the solution is obtained explicitly or implicitly. In the case of implicit members of the family it is seen that the solution of the nonlinear PDE is obtained by solving a linear algebraic system, the coefficient matrix of which is tridiagonal. The behaviors of two methods of the family are examined in a series of numerical experiments. Numerical differentiation and integration procedures are combined to monitor the cumulative amount of drug eliminated into the receptor cell per unit area as time increases. It is found that the use of the equation for the simple membrane model to estimate the permeability coefficient and lag time is warranted even if the system should be described by the dual-sorption model, provided cumulative amount versus time data collected for a sufficiently long time are used. However, being different from the behavior in the simple membrane model, the lag time, which can be estimated in this way, is dose-dependent and decreases with increasing donor cell concentration. On the other hand, the permeability coefficient in the dual-sorption model remains constant irrespective of the donor cell concentrations as in the simple membrane model.
1.
INTRODUCTION
Dose or concentration dependency for percutaneous drug absorption not a well-understood concept. Since the early reports of Chandrasekaran MATHEMATICAL
BIOSCIENCES
108:157-178(1992)
OElsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
is et 157
00255564/92/$05.00
158
K. KUBOTA AND E. H. TWIZELL
al. [l, 21 on the dose-dependent kinetics of transdermal scopolamine base, no paper, to our knowledge, has reported definitive evidence of the dependence of the pharmacokinetics on the dose of a drug in a study in vivo or on the concentration in a study in vitro, for percutaneous absorption of any drug other than scopolamine. Recently, Kubota and Yamada [5] reported on the dose-dependent in vitro transdermal kinetics of a finite dose of timololfree base. Although the detailed mechanism of dose-dependent transdermal timolol kinetics is not known, the findings of Kubota and Yamada [5] suggest that such dose dependency in percutaneous absorption kinetics may not be peculiar to scopolamine. Further studies of some of the features unique to dose- or concentration-dependent kinetics are thus warranted. In the present paper, the dual-mode sorption model of Chandrasekaran et al. [I, 21 is examined, with emphasis placed on the in vitro kinetics when the steady state is attained. It is known that J versus time profiles, where J represents flux from skin to perfusing water, cannot be distinguished from those for the simple diffusion membrane when only profiles for one drug concentration in the donor cell are analyzed. This suggests that it is preferable to use several different concentrations of a drug in the donor cell and to evaluate the presence or absence of the dose or concentration dependency of the percutaneous absorption kinetics, particularly those for lag time. The simple skin membrane model should be used for the analysis only after the dose or concentration dependency of the percutaneous absorption kinetics of the drug have been evaluated. It is noted also that the partition and diffusion coefficients are over- and underestimated, respectively, if the simple membrane model is incorrectly employed. The mathematical formulation of the dual-mode sorption model takes the form of a nonlinear initial/boundary-value problem of parabolic type for which the boundary conditions are functions. The derivation of this mathematical model is given in section 2. Later sections contain computational techniques for the numerical solution of the model equation and numerical results obtained using the associated algorithms. The numerical techniques developed are a family of finite-difference methods that obtain the solution of the nonlinear initial/boundary-value problem either explicitly or implicitly by solving a linear algebraic system at each time step. The results obtained are used to generate mobile (dissolved) and immobile solute concentration profiles and to compute the cumulative amount of drug eliminated into the receptor cell per unit area. A computer program written in BASIC and called MULTI (see Yamaoka et al. [14]) is used to estimate the partition parameter of the drug and the associated diffusion parameter. 2.
THE MODEL
EQUATIONS
The dual-sorption model for the in vitro experiment using excised human or animal skin is depicted in Figure 1. In this model, skin is considered to
PERCUTANEOUS
DRUG ABSORPTION
159
Skin
Donor Cell Diffusion Coefficient
D
I
I
Distance
Concentration,
total
dissolved
C
; KmC
x
;
C,
I KDC
I
immobile
Receptor Cell
CD
K I
CI
0
I
cf+aKDc
1 + aK,C I
I
FIG. 1. Dual-mode sorption model for in vitro percutaneous drug absorption. Total drug concentration C, in skin of thickness L is composed of dissolved (C,) and immobile (C,) concentrations. The skin has diffusion coefficient D for the dissolved (mobile) solute, and distance within the skin is represented by x (0 < x < L). At the donor cell-skin boundary (at x = O), C, is given as C&x = 0) = K,C and C, is given as C,(x = 0) = C;aK&/(l + aK,C), where K, is the skin/receptor cell partition coefficient for the dissolved (mobile) solute, C; is Langmuir’s capacity coefficient, a is Langmuir’s affinity constant, and C is the donor cell concentration. At this boundary, C, is given as C,(x = 0) = K,C, where K, is the apparent skin/receptor cell partition coefficient. At the skin-receptor cell boundary (at x = L), CT, C,, and C, are taken to be zero.
be a single membrane with thickness L and when the drug molecules reach the skin-receptor cell boundary they are removed rapidly, so the concentration of drug in the skin at this boundary may be assumed to have the value zero. The two mechanisms, dissolution and adsorption, are postulated by Chandrasekaran et al. [I, 21 as the processes of permeation through the skin; in these mechanisms, mobile molecules that are freely diffusible and nonmobile molecules that do not participate in the diffusion process are produced. Thus, the total concentration, C,, of a drug in the skin can be expressed as the sum of two components, that is to say, c,=c,+c,,
(1)
where C, and C, are the mobile (dissolved) and immobile solute concentrations, respectively (see also Paul [lo]). The value for C, on the uppermost epidermis, Co( x = 0) is given by CD( x = 0) = K,C,
(2)
160
K. KUBOTA AND E. H. TWIZELL
where C is the drug concentration in the solution in the donor cell and K, is a constant. In the classic steady-state percutaneous absorption experiment (Scheuplein and Blank [12]), an infinite supply of drug is effectively allowed for in the donor cell. Thus, C in Equation (2) can be regarded as a constant. It is noted that K, is the skin/vehicle partition coefficient for the mobile solute concentration and is equal to the effective skin/vehicle partition coefficient, K, [12] only when Co * C,. On the other hand, C, is given as c,=c;ac,/(l+ac,),
(3)
where C; is Langmuir’s saturation constant and II is the affinity constant. As the mobile solute concentration CD becomes large, C, tends to CJ, though it must be remembered that Co cannot exceed the solubility of the drug in skin; and as the affinity constant (I becomes large, C, tends to C; more rapidly with increasing Co. The flux J through the skin per unit area is given by the formula J=-Dax,
ac,
where D is the diffusion coefficient for the mobile solute through skin. The cumulative amount Ae eliminated into the receptor cell per unit area is given by the equation
'ac, x(x=
Ae=Ae(f)=j”J(x=L)dt=-D 0
s
and Fick’s second law of diffusion
for the dual-sorption
ac,
atzD
qdt
(5)
0
model is given as
a%, ax2
’
where
ac, ac, -=,+g[1+ at It follows from (6) and (7), therefore,
(I;j;D)2]%. that
o-0 (see also Frisch [4] and Paul [lo]).
PERCUTANEOUS
DRUG ABSORPTION
161
Chandrasekaran et al. [I] reported the values K, = 1.1, CT = 5 .O mg/mL, a = 0.509 mL/mg, and D = 5.0 x lo- lo cm2/s = O.OOOOOl8 cm2/h, while Michaels et al. [8] give L = 40pm = 0.004cm. Chandrasekaran et al. [1] also report the following set of values of C: (4.4,
19.5,
43.1,
the elements of which are measured
51.4,
in mg/mL
64.01,
(9)
and may be used in Equation
(2). Employing these values of K,, Cy, a, D, and L, the initial/boundaryvalue problem to be solved consists of the partial differential equation (PDE)
[It,1;;;DJ2]%= Ds;
0< x
O,
(10)
together with the initial distribution
c,( x,0) = 0; and the boundary
Oo t>o.
(12) (13)
takes a value from (9). The problem for C, using the family of numerical corresponding value of Cr. can then (I), which give
c,-c,[*+&]. 3.
ANALYSIS
USING SIMPLE
MEMBRANE
MODEL
One of the major interests in this paper is to elucidate the consequence of a case that is likely to occur in actual in vitro experiments. In such a case, the cumulative amount excreted into the receptor cell, which should be explained by the dual-sorption model, is inadvertently analyzed using the simple membrane model. The steady-state pharmacokinetics, observed when the donor cell concentration C is unchanged, can be characterized by (i) the steady-state flux,
162
K. KUBOTA AND E. H. TWIZELL
and (ii) the lag time TL, defined as the value for the t intersect of the asymptote of the Ae versus time curve. Of these two values, TL but not JS, reveals the nonlinear dependence on the donor cell concentration. This is because in the dual-sorption model the concentration gradient at the steady state of CD, dC,/ax, becomes a constant throughout the skin as in the simple membrane model and is given as - Co( x = O)/ L. J in (4) now becomes a constant, J,,, which is given by J,,,
J,, = DC,(x
= 0)/L
= DK,C/L.
(15)
Equation (15) designates that the steady-state flux Js, is proportional to the donor cell concentration C in the dual-sorption model, as in the simple membrane model. The permeability constant P, defined as
P= J,,/C, is thus simply described
as a constant P=
(16)
given by
DK,/L=
k”,L;,
(17)
where k> and L*d are the diffusion parameter and partition parameter defined by Okamoto et al. [9] in the dual-sorption model and given, respectively, as k; = D/L*,
(18)
L; = K,L.
(19)
On the other hand, according to the explicit expression for the lag time TL obtained by a general method for determining asymptotic steady-state solutions to Fick’s second law for permeation [lo], the TL in the dualsorption model is known to be dependent on the donor cell concentration [l, 2, lo] as l/2( OK&)*+
UK&-
(I+ UK&)
( aK&)3
m(l+
aK,C) I
Therefore, if the Ae-time data for more than two different values for the donor concentration C are available, it may not be difficult to see that the nonlinear model is required to analyze the data. However, it has not been determined whether any discrepancy occurs between the observed data and what the model predicts if the data of Ae
PERCUTANEOUS
163
DRUG ABSORPTION
versus time for one donor cell concentration are analyzed by a simple membrane model while the dual-sorption process prevails. In this paper, this problem is examined as follows; several Ae versus time plots are generated assuming that the dual-sorption process prevails in the in vitro diffusion cell experiment, and those data are fitted to the equation of the linear model given below by (see e.g., Crank [3, pp. 42-611 and Okamoto et al. [9])
r=1
kd,t-k--$e
(-k
(-l)‘exp
r%r*t) rz
1*(21)
Another question that has not been elucidated so far is whether the analysis of the Ae versus time data using the simple linear membrane model gives the correct J,, and TL values if the dual-sorption mechanism prevails. This might be a practical problem, particularly when enough of the Ae-time data after the steady state is reached are not available. In such a case, the Ae-time data are sometimes analyzed [9] using (21); that is, in general situations, when Ae is known or has been determined for a discrete set of values {t,, I,, . . . , t,}, the associated partition parameter L,, and diffusion parameter kdl can be estimated using formula (21) and the adaptation of the Levenberg-Marquardt algorithm given in [6], [7], and [ 131. The values TL,, P, , and Jss, are given using the estimated parameter values L,, and k,, in (21) as T,, = 1/6&j,, P, =
(22)
kd,Ld,3
(23)
and J,,, = P,C.
4.
NUMERICAL
SOLUTION
OF THE NONLINEAR
(24)
PDE
Consider the initial/boundary-value problem consisting of the nonlinear parabolic partial differential equation (lo), the initial distribution (1 I), and the boundary conditions (12) and (13). It will be assumed that the parameters C;, a, and K, have been determined and that C takes a discrete value from a set such as (9). Suppose that the interval 0 < x < L is divided into M + 1 subintervals each of width h so that (M + 1) h = L and the time variable t is descretized in steps of length 1. The open region R = [0 < x < L] x [t > 0] and its boundary dR consisting of the lines x = 0, x = L, and t = 0 have thus been covered by a rectangular mesh, the mesh having coordinates x, = mh
164
K. KUBOTA AND E. H. TWIZELL
and t,=nl,wherem=O,l,..., M+landn=0,1,2 ,.... Thetheoretical solution of the initial/boundary value problem {(lo), (1 l), (12), (13)) at the point ( x,, t,) is CD( x,, t,); the solution of an approximating numerical method at the same point will be denoted by 7;. It will be convenient to let aCJ*]-’ [see Eq. (lo)] and then p = DI/h* and q(C,) = [l+ C;a/(l+ and n = to let q~=[l+C;a/(l+ay~)*]-’ with m=1,2,...,M 0,1,2 )... . The derivatives in the PDE (10) will be approximated by the finitedifference replacements
and a2C,(W)
a2
= h-*[@{C,(x-
h,t+l)-2C,(x,t+l)
+C,(x+h,t+l)) +C&+
h,‘)}]
+(1-B){C,(x-h,t)-X&t) + C(h*),
(26)
in which 13 is a parameter 0 < 0 < 1. Evaluating q(C,) point (x,, t,) leads to the family of numerical methods -
pq;ey;::
+ (1+2pq;e)7:+’
-
in (10) at the base
p4;wz:
=pq~(l-e)y~_,+[1-2~q~(l-e)]~~+~q~(l-e)~~+,. The principal benefit gained by evaluating q(C,) at the point tion of any of the methods in (27) and resisting the temptation to at the advanced time level t = t,+ , is that the solution of nonlinear problem can be computed by solving a linear algebraic each time level. This linear algebraic system has the form A(“)r”+l
= B’“‘r”
+,$I),
(27) of applicaevaluate it the given0 system at
(28)
in which A(“) and B(“) are tridiagonal matrices of order M and the vector CC”)= [c’,“‘, 0, . . . ,OIT is also of order M with c\“) = pqyK,C; the boundary condition (12) is responsible for the absence of r3 from cr. The matrices A’“’ and B(“) are given in the Appendix. A fast, direct, well-known method for calculating P ’ + ’ from (28) is described in Twizell [13], [Ch. 11. The local truncation errors associated with the family of methods given
PERCUTANEOUS
DRUG ABSORPTION
16.5
by (27) are 0(/r* + 1) for 0 # l/2 and 0(/z* + 1*) for 8 = l/2. The von Neumann necessary condition for stability requires the localized restriction
where the subscript S refers to the spectral norm (see Twizell [ 13, Chs. 7 et seq.]). Clearly, any stability restriction will require modification at each time step because the PDE is nonlinear. The choice of which value of 0 to choose is influenced by a number of factors. Putting 8 = 0 allows the solution to be obtained explicitly because then A(“)=1 (the identity matrix of order M) in (28), but a linear stability analysis imposes the restriction pqL’ < l/2 in integrating from time t = t, to time t = t,,,, where qi’ = max,( qi 1. Clearly, this method of the family (27) is competitive only when the user can afford to use a small time step. Putting 8 = 1 gives a fully implicit method [B’“‘=I in (28)] that will be especially suitable for solutions that decay rapidly, and putting 0 = l/2 gives a more accurate implicit method that may not be suitable for use when there are discontinuities between initial conditions and boundary conditions, as indeed there are with the present problem [see (11) and (12)] if 1 is too large with respect to h (see Twizell [13, Ch. 81). Alternative methods of obtaining a numerical solution include the use of the so-called method of lines in which the parabolic initial/boundary-value problem ((lo), (1 l), (12), (13)) is transformed into a first-order initial-value problem. This approach is particularly useful when the PDE is linear and is described for a linear model of percutaneous drug absorption in Twizell [ 13, Ch. 81. In the case of a nonlinear PDE, the first-order initial-value problem must be integrated by a numerical method designed to solve such problems. These methods include Runge-Kutta methods, linear multistep methods employed in predictor-corrector combinations, hybrid methods, and multiderivative methods. Stability restrictions depend on the method being used and on the PDE itself and can be determined using the theory of weak stability for ordinary differential equations (see Twizell [ 13, Ch. 51). 5.
NUMERICAL
RESULTS
Using the values of K,, C;, a, D, and L reported in Section 2, the initial/boundary-value problem to be considered consists of the partial differential equation l+
2;;;;CDJZ
acLJ ac, 1 at=0.0000018
at
,
Oo,
= 0;
where C is obtained from (9). In order to examine the behavior of the family of computational methods given by (27), a number of numerical experiments were carried out. The parameter C was given the value 4.4 mg/mL first of all, the space step h was given the value 0.0002 cm, so that M = 19, the time step was put equal to 0.1 h, and the parameter 8 was given the values 1 and l/2. The solution y (an approximation to CD) of {(30), (31), (32), (33)) was then calculated for both values of 0 until the steady state was reached. It was found that the steady-state solution, the straight line joining the points (O,K,C) = (0, l.lC) and (L,O) = (0.004, 0) in the (x, Co) plane, was reached after about 170 time steps. The output of the computer program gave the values of y = Co, where x = 0, 0.0002, 0.0004,. . . ,0.004 cm, the corresponding values of C, being computed using (14). The numerical results for Co and C, are depicted in Figures 2 and 3, respectively, at times 0.5 h, 1.0 h, 2.0 h, 3.0 h, 4.5 h, 6 h, and 15 h (that is, after 5, 10, 20, 30, 45, 60, and 150 time steps). To the scales used in Figures 2 and 3 the corresponding curves for 8 = l/2 at the higher values of t are superimposed. An example of the closeness of the results for Co, for 0 = l/2 and 0 = 1, at time t = 4.5 h is revealed in Table 1. Using the fully implicit method (0 = l), the solution is seen to decay more quickly as x increases than when the implicit method (0 = l/2) is used, as was noted earlier. At the lower values of t the differences between the values of Co were more marked, as can be seen in Figure 4, where the distributions for f3 = l/2 and 0 = 1 at time t = 0.5 h (after five time steps) are depicted. The theoretical solution of the initial/boundary-value problem is not known, so it is not possible to say which of the two values of 8 gives the “better” set of results. Some indication may be obtained by refining the mesh. For instance, with 0 = 1, reducing the space step h by a factor of 3 and the time step 1 by a factor of 10 reduces the local truncation error by a factor of about 10. However, the storage and CPU times will be increased substantially by such moves, and the user may well be satisfied with results obtained on a coarser mesh.
PERCUTANEOUS
167
DRUG ABSORPTION
CD(mg/ml)
0
8
24
16
32
40
x(w) FIG.
2. Values of C,
using 9 = 1, I = 0.1 h, h = 2 pm, M = 19, and C = 4.4
mg/mL.
The numerical experiments were then repeated for the other values of C listed in (9). There was little appreciable difference between the distributions computed using any of the five values of C, with either 8 = l/2 or 1, at any time t. Once again the differences were greater at the lower values of t, although only those for the least and greatest values of C (4.4 and 64.0) were markedly different. These distributions of CD at time t = 0.5 h (after five time steps), normalized as Co /K,C, are depicted in Figure 5. The differences between the distributions CD are due to the nonlinear dose response of the model; when C = 64.0 mg/mL, the distribution precedes that of C = 4.4 mg/mL because the small fraction of drug that has diffused
168
K. KUBOTA AND E. H. TWIZELL CThg/ml)
9
1
0
8
24
16
32
40
XC!Jd FIG.
3. Values of C,
using 0 = 1, I = 0.1 h, h = 2 wm, M= 19, and C = 4.4
mg/mL.
into the deeper portion is immobilized when C = 64.0 mg/mL as shown in (1) and (3). This difference would be reflected eventually by the dosedependent behavior of the lag time in (20). The value of Ae at time t can be estimated from Equation (5) by employing a numerical integration formula. First of all, however, X, /ax must be estimated at x = L using a numerical differentiation formula. The numerical method used to solve the model PDE (10) incurred an 0( h*) error in approximating the space derivative by (26). There is no point, therefore, in using a numerical differentiation formula that enjoys better than second-order accuracy to estimate 0” = NJ,< L, nl)/ax. To this end, the second-order backward difference formula P”=(r~-,-4y~++3y~+,)/2h+O(h2) can be used (see Twizell [13; p. 921); note yl;fr+, = 0 from (13).
(34)
PERCUTANEOUS
169
DRUG ABSORPTION TABLE 1
Values of C, (mg/mL) at Time t = 4.5 h Using (27) with 19= l/2 and 1 and C = 4.4 mg/mL
e=i
e=1/2
x(e) 0
4.840 4.285 3.735 3.195 2.671 2.166 1.685 1.229 0.799 0.393 0.000
4 8 12 16 20 24 28 32 36 40
4.840 4.280 3.726 3.183 2.656 2.150 1.669 1.215 0.789 0.387 0.000
It is an easy matter to adapt the computer program to compute 6C,(L,n1)/6x at each time step, following which Ae can be computed using the trapezoidal rule Ae(nl)=-+Dl
/3°+2nc1/3q+pn I
(see
Twizell
6.
RESULTS METHOD
q=l
[13, p. 991). OF THE
ANALYSIS
BY THE
1
+0(1*)
SIMPLE
(35)
MEMBRANE
The Ae versus time plots produced by (35) for C = 4.4 (triangles) and 64.0 (circles) mg/mL, taking h = 0.0002 cm and I= 0.1 h, are shown in Figure 6a. The solid lines in Figure 6a are the Ae(t)-time curves given by (21) fitted to the hourly Ae versus time plots (up to 24 h) generated by (35). The nonlinear least squares computer program written in BASIC and called MULTI (see Yamaoka et al. [14]) is used. The plots generated by Equation (35) are well explained by the “incorrect” model equation (2 1). To illuminate that Js, and P are dependent not on C but on lag time TL (20), or that T,., (22) is dependent on the donor cell concentration, Ae(t)/C versus time plots are shown in Figure 6b. In Table 2, the values of k,, and L,, estimated by fitting the hourly Ae(t) versus time until t, (t, = 4, 6, 8, 10, 12, 18, and 24 h), generated by (35), to Equation (21) are shown. The TL, (22), P, (23), and J,,, (24) values are also shown in Table 2 and compared to TL (20), P (16), and J,, (15) in the dual-sorption model, respectively When C = 4.4 mg/mL, the numerical values for TL, P and Js, were 3.33 h, 0.495~ 1O-3 cm/h, and
K. KUBOTA AND E. H. TWIZELL CD(mg/ml)
0
8
24
16
32
40
*(urn)
FIG. 4. Values of C, after five time steps using I = 0.1 h, h = 2 pm, M = 19, and C = 4.4 mg/mL.
2.18 pg/(cm’*h), respectively, and when C = 64.0 mg/mL those values were 1.81 h, 0.495 X 10e3 cm/h, and 31.7 pg/(cm2*h), respectively. As is shown in Table 2, the values for P,, JSS,, and TL, (particularly that for TL,) particularly when are overestimated and the k,, value is underestimated, the Ae versus time plots collected for the short-time duration containing no or few data after the steady state is attained are fitted to (21). However, when the Ae(t) versus time plots collected for the long time duration are used, those values for TL, (22), P, (23), and J,,, (24) approach TL (20), P(16), and J,, (15), respectively, even if the “incorrect” model equation (21) is used. For instance, in Table 2, the error of the estimated TL, value (22) is within 20% of the true TL value (22) when the Ae versus time plots
PERCUTANEOUS
171
DRUG ABSORPTION
c,/ CKp 1.0
0.8
-
0.6
-
0.4
-
0.2
-
0
I
I
8
16
I
24
32
40
x (urn)
FIG. 5. Values of C, (normalized to C,/K,C) at time t = 0.5 h using 0 = 1, I = 0.1 h, h = 2 pm, and M = 19. Values of C are in mg/mL.
collected for more than twice the estimated TL, value are analyzed. When the Ae versus time plots collected for more than four times the estimated TL, value are analyzed, the error is within 10% of the true value. 7.
DISCUSSION
The primary conceptual parameter value estimated in steady-state percutaneous drug absorption kinetics is the permeability constant P defined by Scheuplein [II] and Scheuplein and Blank [12] as in (16). Based on the simple membrane model, there are two alternative approaches for estimating diffusion and partition parameters. In the first (approach A), the diffusion parameter is estimated from the lag time T, [I].
172
K. KUBOTA AND E. H. TWIZELL 800
cy -
600
E 3J 2
400
$ 200
0
6
12 Time
18
24
(h)
@) FIG. 6. (a) Profiles of cumulative amount excreted into receptor cell, Ae, plotted versus time, generated in the dual-sorption model where donor cell concentrations C are 4.4 (A) and 64 (0) mg/mL. The solid lines represent the best fits predicted by Equation (21) for the simple skin membrane model. (b) Ae/C versus time profiles, symbols and solid lines represent the same profiles as in (a).
As given concentration
in (20), the lag time TL is dependent C. Defining the lag time prolongation
on the donor cell factor S to be the
quantity
S=l+6C~al/2(aK,C)2+aK,C-(1+aK,C)In(l+aK,C) I ( aKLG3
7 (36)
PERCUTANEOUS
DRUG ABSORPTION
173
TABLE 2 Values of k,, , L,,, P,, J,, and TLl Estimated When Hourly Ae- t Plots up to Time f,(h) Generated by the Dual-Sorption Model Are Fitted to Equation (2 1) Describing the Simple Diffusion Model.
tz (h)
Qa @-‘I
Plb
&I,=
JSS,'
(cm x 10-3)
(cm/h x 10e3)
C = 4.4 mg/mL 4 0.0356 6 0.0395 8 0.0422 10 0.0441 12 0.0454 18 0.0474 24 0.0482 _e 0.0507f
20.7 15.4 13.1 12.0 11.4 10.6 10.3 9.179
0.736 0.606 0.553 0.530 0.517 0.503 0.499 0.495h
C = 64 mg/mL 4 0.0781 6 0.0843 8 0.0871 10 0.0885 12 0.0894 18 0.0906 24 0.0911 _e 0.0953f
7.14 6.13 5.81 5.66 5.58 5.48 5.44 5.20g
0.557 0.517 0.505 0.501 0.499 0.497 0.496 0.495h
[~gI(cm*.h)l 3.24 2.67 2.44 2.33 2.28 2.21 2.20 2.18 35.7 33.1 32.3 32.1 31.9 31.8 31.7 31.7
TLld (b) 4.68 4.22 3.95 3.78 3.67 3.52 3.46 3.29’ 2.13 1.98 1.91 1.88 1.87 1.84 1.83 1.75’
a Estimated from (21).bCalculated as k,,L,, (23). ‘Calculated as P,C (24). d Calculated as 1/6k,, (22). eIdeal case where P, (23) and TL, (22) are estimated exactly as P (17) and T, (20) in the dual-sorption model, respectively. f Estimated as 1/6T,. ‘Estimated as P/k,. hGiven in (17). i Given in (20).
a new conceptual
parameter,
TL, defined as
T;= T,/S=
L2/6D=
1/6k;,
may be introduced. Because S approaches unity as C approaches infinity, Ti is the fictive lag time obtained when the donor cell concentration is assumed to be infinite. In this section, a mechanism is outlined for estimating the parameter values, using the simple membrane model, from the data where the dualsorption process prevails. First, though the lag time TL is dose-dependent, the profiles of Ae versus time in Figure 6, observed for dual-mode sorption models, are difficult to distinguish from those for simple membrane models whenever
174
K. KUBOTA
AND E. H. TWIZELL
results for only one concentration of a drug in the donor cell are analyzed. It is thus quite likely in such cases that the diffusion and partition coefficients will be estimated using a simple membrane model even when a dosedependent model is required. The apparent diffusion parameter k,,, obtained when the simple membrane model is appropriate, is given, using S defined in (36), as
k,,= k’,lS, following equation
which the apparent partition
parameter
(38) L,, can be defined by the
L,, = L*dS. Second, the results shown in Table 2 indicate that the values for TL, (22), P, (23), and J,,, (24) of the dual-sorption model are well tracked as the values in the dual-sorption model, T, (20), P (17), and Js, (15)) respectively, even when the incorrect model equation (21) is employed, provided that Ae versus time plots collected for a sufficiently long time are analyzed. Thus, it is warranted to use the ‘ ‘ incorrect’ ’ model equation (2 1) to estimate the TL (20), P(17), and Js, (15) values in the dual-sorption model from Ae versus time data obtained for each donor cell concentration. This feature is of practical significance because TL, P, and Js, are sometimes estimated by fitting the Ae versus time data to Equation (21) (see Okamoto et al. [9]) even if the data after the steady state is attained are available. This procedure is useful because the actual Jss values in the experiments are usually fluctuating, and in some cases it is difficult to find the time after which J,, can be considered to have been attained. However, it is notable that if few or no Ae-time data after the steady state is attained are available, Ld,, P,, TL, and J,,, are generally overestimated and k,, is underestimated (Table 2). Third, if the TL, /S values [which should be TL in (37)] for several donor concentrations are of a similar magnitude to one another, it suggests that the dual-sorption process is the likely mechanism causing the nonlinear pharmacokinetic profiles. This is because the values of the lag-time prolongation factor S in (36) are determined by a, C;, and K, obtained from the sorption isotherm experiment, and the consistency of TL /S values for different donor cell concentrations indicates that two nonlinear features (observed in the sorption isotherm and diffusion experiments) are explained as two aspects of a dual-sorption model. In the second method (approach B) to estimate the diffusion and partition parameters based on the single membrane model, the partition parameter is estimated from the sorption isotherm experiment as reported by Scheuplein
PERCUTANEOUS
DRUG ABSORPTION
175
[ 1 I 1. In situations where a dual-sorption model is appropriate, when equilibrium is attained in the sorption isotherm experiment, the total concentration of drug, C,, is constant throughout the skin (the epidermal sheet) and is given by the equation
(40) where C is the concentration Defining R by
of the drug in the solution
R=l+Cfa/(l+aK,C), the apparent partition
coefficient K,
The corresponding
partition
K,
and defining
kd2 as k,
parameter
= P/L,?
(41)
is given as
= CT/C=
Ld, = K,L
at equilibrium.
(42)
K,R. Ld2 is given by
= K,RL
= L*,R
(43)
it follows that
kd2 =
k”,lR,
(44)
where L>, ki, and P are given in (19), (18), and (17), respectively. Though the above expressions for Ld2 and kd2 are dependent on the concentration of solution at equilibrium, C, it is quite likely that only one concentration at equilibrium will be employed to estimate the partition coefficient K,, if it is not realized that the dual-sorption process occurs in the experiment. As is seen in Table 3, when the parameter values for scopolamine are employed (see (9) and [2]), the values of R and S are of the same order but are not identical. Thus the parameter values estimated by
TABLE 3 Values of R and S Corresponding to C in (9) [2] c
R
s
4.4 19.5 43.1 51.4 64.0
1.735 1.214 1.101 1.085 1.069
2.221 1.481 1.255 1.219 1.181
176
K. KUBOTA AND E. H. TWIZELL
approaches A and B would be of similar magnitude to each other if the concentration employed at equilibrium is of the same magnitude as that in the diffusion experiment. Thus, in this case, the fact that the dual-sorption process prevails may again be easily missed. In conclusion, it is important to employ various donor cell concentrations in the diffusion experiment and various concentrations at equilibrium in the sorption isotherm experiment to know whether the dual-sorption process occurs in the experiment. If a single concentration is employed in the diffusion and sorption isotherm experiments, an important feature, that the process should be described by the dual-sorption model but not the simple diffusion model, might not easily be recognized. An early priority is the ability to estimate the parameters CT and a, first used in Equation (3). Chandrasekaran et al. [l, 21 reported values of C, for the set of values of C given in (9). The parameters Ci and a may be estimated using the data pairs (C, Cr.) with Equation (40); the Levenberg-Marquardt algorithm for doing so is given in Levenberg [6], Marquardt [7], and Twizell [ 13, Ch. 21. The research for this paper was begun while one of us (E. H. T.) was a visitor to the National Medical Center (NMC), Tokyo. The substantial grant received from The Royal Society of London that made the visit possible, and the help of Dr. Takashi Ishizaki (NMC) in organizing the visit, are gratefully acknowledged. We are very grateful to two anonymous referees for their many valuable comments and to Mrs. B. Yates and Ms. M. Demmar of Brunei University for their help in preparing the typescript. APPENDIX.
THE MATRICES
A’“’ AND B(“)
The matrices A’“’ and B(“) are given by
where a:=
1+2pqG8,
b$=
-pq;e,
PERCUTANEOUS
DRUG ABSORPTION
177
and
0,” B(n)
0
Pz”
= 9
0
KG
44-I Kl
I
PZG1 44
with
a;=1-2pq;(l-e),
The meanings
of p, q;
P;=Pq;(l-q.
and ~9are described
in the main body of the paper.
REFERENCES 1 S. K. Chandrasekaran, A. S. Michaels, P. S. Campbell, and J. E. Shaw, Scopolamine permeation through human skin in vitro, Am. Inst. Chem. Eng. J. 21:828-832 (1976). 2 S. K. Chandrasekaran, W. Bayne, and J. E. Shaw, Pharmacokinetics of drug permeation through human skin, J. Pharm. Sci. 67: 1370- 1374 (1978). 3 J. Crank, The Mathematics of Diffusion. 2nd ed., Clarendon Press, Oxford, 1975. 4 H. L. Frisch, The time lag in diffusion, J. Phys. Chem. 61:93-95 (1957). 5 K. Kubota and T. Yamada, Finite dose percutaneous drug-absorption: theory and its application to an in vitro timolol permeation, J. Pharm. Sci. 79:1015-1019 (1990). 6 K. Levenberg, A method for the solution of certain nonlinear problems in least squares, Q. Appl. Math. 2:164-168 (1944). 7 D. W. Marquardt, An algorithm for least squares estimation of nonlinear parameters, SZAMJ. 11:431-441 (1963). 8 A. S. Michaels, S. K. Chandrasekaran, and J. E. Shaw, Drug permeation through human skin: theory and in vitro experimental measurement, Am. Inst. Chem. Eng. J. 21:985-996 (1975). 9 H. Okamoto, M. Hashida, and H. Sezaki, Structure-activity relationship of I-Alkylor 1-alkenylazacycloalkanone derivatives as percutaneous penetration enhancers, J. Pharm. Sci. 77~418-424 (1988). 10 D. R. Paul, Effect of immobilizing absorption on the diffusion time lag, J. Polym. Sci. (Part A-2) 7:1811-1818 (1969). 11 R. J. Scheuplein, Mechanism of percutaneous absorption. I. Routes of penetration and the influence of solubility, J. Invest. Dermatol. 45:334-346 (1965).
178 12
K. KUBOTA AND E. H. TWIZELL
R. J. Scheuplein and I. H. Blank, Permeability of the skin, Physiol. Rev. 51:702-747 (1971). 13 E. H. Twizell, Numerical Methods, with Applications in the Biomedical Sciences, Ellis Horwood, Chichester and Wiley, New York, 1988. 14 K. Yamaoka, Y. Tanigawara, T. Nakagawa, and T. Uno, A pharmacokinetic analysis program (MULTI) for a microcomputer, J. Pharmacobio-Dyn. 4:879-885 (1981).
