- Email: [email protected]
Diffusional release of a dispersed solute from a spherical polymer matrix
journal of MEMBRANE SCIENCE ELSEVIER
Journal of Membrane Science 115 (1996) 171 - 178
Diffusional release of a dispersed solute from a spherical polymer matrix M.J. Abdekhodaie, Y.-L. Cheng * Department of Chemical Engineering and Applied Chemistry, University of Toronto, 200 College Street, Toronto, Ontario, Canada, M5S 1A4 Received 12 July 1995; revised 30 October 1995; accepted 18 December 1995
Abstract The exact analytical solution and an approximate solution based on the variational formulation method have been obtained for the release kinetics of a solute from a spherical polymeric matrix in which the initial solute loading A is above the solubility limit CS. Comparisons between these new solutions and previously reported approximate solutions based on either pseudo-steady state assumptions or the refined integral method indicate that the approximate solutions are very accurate when A / C s is much larger than unity. However, as drug loading decreases, deviations between the approximate solutions and the exact solution increase so that the approximate solutions are least accurate in the limit A / C s ~ 1. Keywords: Controlled release; Diffusion; Theory; Exact solution; Approximate solution
1. Introduction
been used to analyze the release kinetics of a dis-
The diffusional release of a solute from a polymeric matrix in which the initial drug loading A is greater than the solubility limit C s is a problem of interest in controlled release drug delivery. This
persed solute from polymeric matrices. So far two approximate solutions have been developed to analyze the release kinetics of dispersed solute from a polymeric matrix with spherical geometry [4,5]. In both of these solutions, perfect sink conditions were assumed, and matrix swelling of the matrix, the concentration dependence of the solute diffusion coefficient, and the external mass transfer resistance
problem can be formulated by considering a moving diffusional front separating the core containing undissolved drug and the partially extracted region, Moving boundary problems arise in a wide range of application areas, including diffusional release [1], heat transfer involving melting phase transition [2], and diffusion-controlled growth of particles [3]. A1though moving boundary problems have been widely studied, only a few mathematical descriptions have
* Corresponding author. Tel: (416) 978-5500; Fax: (416) 978-
8605; E-mail: y|[email protected],
were neglected. Higuchi [4] obtained an approximate solution based on the pseudo-steady state assumption that at any instant in time, steady state concentration profiles existed in the (time-dependent) partially extracted region. Since the diffusional front velocity should increase as solute loading decreases, the pseudo-steady-state assumption is expected to be least applicable at low solute loadings. Lee [5] used the refined integral method to obtain approximate
0376-7388/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved PII S 0 3 7 6 - 7 3 8 8 ( 9 6 ) 0 0 0 1 8 - X
172
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
analytical solutions for release from planar and spherical matrices containing dispersed drug. T h e approximate solution to the slab problem agreed well with the known exact solution obtained by Paul and McSpadden [6], with the largest deviation at drug loadings near the solubility limit. No exact solution for the spherical problem presently exists for comparisons to be made. Although Lee's solution is elegant and the result is easy to use, the refined integral method of solution is complex. It is therefore the primary objective of this study to develop an exact analytical solution for the release kinetics of dispersed solute from a sphere into an infinite sink. In addition, a simpler procedure for obtaining an approximate solution to this problem based on the variational formulation method will be presented. Comparisons between the exact solution and all known approximate solutions will be made.
Partially
/ExtractedRegionUnextracted / / ~
/
c
o
r
e
\ , . .~. . i ~ R
$
' r
Fig. 1. Schematic diagram of drug distribution in a partially extracted pellet: R, position of moving diffusion front; S, position of sphere surface.
matrix. Assuming perfect sink condition at all times, the initial and boundary conditions are: C[S, t] = 0
(2)
C [ R ( t ) , t] = C s 2. Analysis
(3)
OC
dR
D--=(A-Cs)-The diffusional release of a dispersed solute from a polymeric matrix with spherical geometry into a perfect sink fluid is considered. The diffusion coefficient is assumed to be independent of concentration, and solute diffusion is assumed to be the rate controlling step rather than polymer swelling or drug dissolution. Polymer swelling kinetics can be neglected in systems where little swelling occurs, or where swelling occurs rapidly relative to solute diffusion such that diffusional release essentially occurs from a pre-swollen sphere, The drug distribution in the sphere after a certain time of release is illustrated in Fig. 1 where S is the radius of the entire (possibly pre-swollen) pellet, R is the radius of the unextracted core and r is the polar radial position variable. The concentration is essentially constant in the unextracted core ( r < R), while the concentration in the partially extracted region, R < r < S, is a function of r and is determined by transient diffusion according to Fick's second law: Ot
r 2 Or ~
--~r ]
(1)
where D is the diffusion coefficient, and C is the concentration of dissolved solute concentration in the
at r = R
Or
(4)
dt
R(0) = S
(5)
where A is the initial solute loading per unit volume, and C~ is the solubility limit of the solute in the • matrix. 2.1. Exact solution
By using reduced dimensionless variables defined as: r
Dt
R
~:= ~
~"= -~-
F = -~
r[C~
0 = -~ ~ - ~ )
(6)
the defining equations can be reduced to their dimensionless forms: 00 020 (7) 0r 0~ 2 0[1, ~-] = 0 (8) 0[ F(~-), 'r] = F (9) 1 00
1
F (0) = 1
-
1
at ~ = F
(10)
(11)
173
M.J. Abdekhodaie, Y.-L. Cheng/ Journal of MembraneScience 115 (1996) 171-178
By defining K as: ( A ) K= ~ - 1 Eq. (10) becomes: 1 00 1 dF F 0~ F K d~"
(12)
at ~-- F
(13)
Eq. (7) can be reduced to an ordinary differential equation by using an appropriately defined combined variable U: (14)
1-~:
U=
x/~
such that the dimensionless concentration 0 becomes a total function of the combined variable:
In principle, Eq. (20) can be solved numerically with the initial condition on F g i v e n i n E q . ( l l ) . However, since Eq. (20) gives d I / d r ~ oo at ~"= 0, a value of F at 1- other than zero is needed as the initial condition for the numerical solution of Eq. (20). At very short times (small ~'), the position of the diffusion front is very close to the surface of the sphere and the partially extracted region is very thin. It is well-known that curvature effects are negligible in thin films such that thin cylindrical or spherical films behave similarly to thin planar films. It is therefore reasonable to expect that F(T) at early times for planar and spherical systems should be similar. Paul and McSpadden's [6] planar geometry solution gives: f ~ [ 12@TF]
[(1-/")
2
[ ]1-/"
The diffusion equation, Eq. (7), then becomes:
(21)
(U) d20 dU 2
dO dU
(16)
Integrating directly twice, we obtain: rf( U ] 0 = B t f ~ e ~ -~- ] + B2
(17)
where B l and B2 are constants of integration. By using boundary conditions (8) and (9), 0 becomes: [~] /" erf 0=
and hence: C -- =
Cs
SA M= = - -
1-
(19)
F] ev[ 2x/-T J
Using the last boundary condition, Eq. (13), gives: (1 - / " ) 2
exp
For any combination of solute loading and solubility, Eq. (21) can be used to calculate a value for iV at an arbitrarily small value of ~, (e.g. ~'= 10-6). These corresponding iV and ~- values can then be used as an initial condition (/"0, % ) i n combination with Eq. (20) to numerically solve for the time profile of the dimensionless diffusion front position in a sphere F(~-). The cumulative amount of solute released per unit area at any time t, M t, is related to the diffusional front position R(t) by a simple mass balance:
The total amount of drug, or the amount releasable at infinite time, per unit area, M~,, is given by:
/"erf[~] r
C~
4r
x/__~v/~ erf [_~___r ]
(23)
3
The fractional release at time t in terms of dimensionless quantities then becomes: M--L= (1 - / " 3) M~
1
-~
d /"
3 Cs F
(20)
A
S/~:erf(~_7~)d~:
(24)
M.J. Abdekhodaie, Y.-L Cheng /Journal of Membrane Science 115 (1996) 171-178
174
By substituting F(~-) profiles from the solution of Eq. (20) into Eq. (24) and integrating, the release kinetics can be obtained in terms of the time dependence of fractional release.
.
0.00
. 0.05 .
.
.0.15 .
0.10
0.20
0.25
2.5
3.0
1.0 ~
\
0.8 NCs = 2
2.2. Approximate solution using variational formulation method
0.6-
By taking Fick's second law as the Euler equation of the desired varational formulation and assuming a quadratic concentration profile, the following equations can be obtained to calculate the fractional release (see details in Appendix A): Mt - 1 - r
/
3 -(1
~r
Cs
/
\\ " \\ 0.4
1
~
0.5
1.0
1.5
2.0
(25)
where /" and a for each value of r can be obtained by solving the following equations:
= KF a( a - 3) (l-F):(1 ) (l-a) + F ~ ~-
Ncs=11
0.2
0.0
a
\
I
Ms + ~ - ( / ' + 1)
\
\
F
Fig. 2. F ( ~ ' ) with initial condition F 0 at ~" = ~'o o b t a i n e d f r o m P a u l - M c S p a d d e n p l a n a r g e o m e t r y solution for two values of
A / C s : (a) A / C s = 2 ; - - , ~'0 = 1 0 - 8 to 10 - 2 and . . . . . ro=0.1;(b) A / C s = l l ; - - , ~ - o = l O - 8 to 10 - 2 , - - ~'0=0.1
(26) pected, greater sensitivity was observed for the lower (27)
3. Results and discussion The exact solution for the fractional releases, MJM=, requires knowledge of the dimensionless moving front position F(~-). As described above, at early times, curvature effects can be neglected and F ( z ) from the known exact solution for the planar geometry problem can be adopted for the spherical problem. To use this approach, it is necessary to determine if the solution for _F(~-) is sensitive to the initial conditions selected. Fig. 2 shows the F(~-) profiles obtained from Eq. (20) using varying values of (/'0, %) pairs from Eq. (21) as the initial condition. These calculations were carried out for two drug loadings, A/Cs = 2 and A / C s = 11. The results indicate that for A / C s = 11, indicate that the F(~') profiles are insensitive to the initial condition selected for ~'o ranging from 10 -8 to 0.2. As ex-
drug l°ading ( A / C s = 2 ) " N° differences in I'('r) profiles were observed for ~'0 ranging from 10 -8 to 10 -2. For ~'0 = 0.1, however, a significantly different / ' ( 7 ) profile was calculated. At lower drug loadings, the diffusional front should move inwards at a higher velocity so that curvature effects should become significant at an earlier time. These results confirm the validity of our assumption that curvature effects can be neglected at early times. The results presented below were obtained using / ' values corresponding to ~-= 10 -6 from Eq. (21) as the initial conditions for Eq. (20). Fractional release profiles for various solute loading levels are shown in Fig. 3. The results show a non-linear dependence on r. As expected, as A / C s increases, the cumulative fractional release curves become increasingly sustained. In Fig. 4, fractional release is plotted as a function of z ( C J A ) . Interestingly, the release profiles collapse onto approximately the same curve. This normalized curve can be used to find fractional release profiles for any solute loading. To calculate fractional release based on the ap-
M.J. Abdekhodaie, Y.-L Cheng/Journal of Membrane Science 115 (1996) 171-178 NCs=2
7
11
21
0.8
1.0
1.0
F 0.5
0.5
175
101
NCs=7
0.6
0.0
Mt/M~
0.0
0.4
1.0-
i
~
~
0.1
0.2
0.3
0.0 0.4
[
0.0
0,3
0.6
0.9
1.2
0
1
2
3
4
1.0
0.2 r
0.0 0
I
I
I
1
2
3
0.5
0.5
NCs=I 1 4
0.0
0.0 0
Fig. 3. Fractional release versus ~- for a dispersed solute in a spherical polymer matrix at various solute loading levels based on the exact solution,
1
2
Fig. 5. Diffusion front position vs. time for a dispersed solute in a spherical polymer matrix at various solute loading levels: ,7, Higuchrs solution; (3, Lee's solution; [3, variational formulation solution; z~, exact solution.
proximate variational formulation solution, a direct elimination of F and " a " from Eqs. (25), (26) and (27) is not very convenient. It is more practical to correlate the fractional release with time through
NCs=2
3
101
1.0
0.8
0.6 Mr/M= 0.4
0.2-
o.0o.oo
'
o.04
'
o.o8
' o.12
o.16
x(C~A)
Fig. 4. Fractional release versus ~-(Cs / A ) for a dispersed solute in a spherical polymer matrix at various solute loading levels based on the exact solution,
their dependencies on these parameters. The pseudo-steady-state result given by Higuchi [4] and the solution presented by Lee [5] also require such a correlation through F . The dimensionless time dependence of the diffusion front position, F , and the fractional release, M i M e , are presented in Figs. 5 and 6, respectively. The exact solution and the variational formulation approximate solution presented in this work are shown along with Higuchi's and L e e ' s aproximate solutions. To compare these solutions, the deviation between the fractional release based on each of the approximate solutions and the exact solution at equally spaced values of r throughout the release duration were calculated, then averaged. The timeaveraged deviation at various drug loading levels are shown in Table 1. Increasing A / C , results in reduced deviations. For example, the deviations of the Higuchi, Lee and variational formulation approximate solutions from the exact solution for A / C ~ = 3 are 7.15%, 3.91% and 5.14%, respectively. The corresponding deviations for A / C ~ = 2 1 decline to 0.91%, 0.34% and 0.54%. In all three approximate solutions, an assumption had to be made regarding the concentration profile of the solute in the partially
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
176
1.o f
1.o f
M,/M~0.5
0.5 A/Cs=3
,
,
,
0.1
0.2
0.3
0.0 0.0
NCs=7
0.0 ¢ 0.4
1.0 ~
0.0
f
0.S NCs=I 1
0.0
i
0
0.3
, 0.6 •
, 0.9
1.2
1.0 ~
/
~/M~ 0.5
,
1
0.0
2
0
NCs=21
I
i
1
2
3
Fig. 6. Fractional release vs. time for a dispersed solute in a spherical polymer matrix at various solute loading levels: V, Higuchi's solution; ©, Lee's solution; [2, variational formulation
solution; z~, exact solution,
extracted region. In the Higuchi solution, a pseudo steady state profile was adopted, while a quadratic profile was used by Lee and in the present work. A parameter " a " , given in the appendix by Eq. (A12) was required in the quadratic profiles. This parameter was used to match the assumed profile to the position of the moving diffusion front. Although this parameter should vary with time, it was assumed to be time-independent in both approximate solutions. Since the moving front velocity increases as the solute loading decreases, the pseudo-steady-state assumption used in Higuchi's approximate solution, as Table 1 Comparison of average absolute percentage deviations of different approximate solutions from exact solution for a spherical polymer
matrix A/C S
Higuchi error (%)
well as the assumption of time-independence for the parameter " a " in the quadratic profiles, are less valid at low solute loadings, leading to greater inaccuracies at low loadings. Of the three approximate solutions, Lee's refined integral approximate solution was closest to the exact solution. The accuracy of refined integral and variational formulation solutions depends on the form of the assumed concentration profile. The assumption of a quadratic concentration profile resulted in a more accurate solution using the refined integral method for this particular problem. More accurate solutions can be achieved with improved assumed concentration profiles.
Lee
Variational
error (%)
formulation
solution
error (%)
3
7.1521
3.9082
5.1418
7 11 21
3.6999 1.6486 0.9072
0.7815 0.6361 0.3393
1.7284 1,0131 0,5379
4.
Conclusion
An exact moving boundary solution has been presented for the diffusional release of dispersed solutes from a spherical matrix into an infinite sink medium. This exact solution permits the validation of previously available approximate solutions presented by Higuchi [4] and Lee [5]. In addition, an alternative approximate solution based on a simple variation formulation method is presented. The approximate solutions have some advantages for describing release kinetics due to their relative simplicity. These solutions are very accurate when A / C s is in great excess of one, but in the limit of A ~ Cs, deviations from the exact solution increase.
5. List of symbols A initial solute loading per unit volume Bi constant of integration C solute concentration in the matrix at any time C s equilibrium solute solubility in the m a t r i x D solute diffusion coefficient in the matrix K initial solute loafing level defined by Eq. (12) M t amount of solute release per unit area M ~ total amount of solute release per unit area r spatial coordinate R position of the moving diffusion front S radius of the sphere t
release time
U combined variable defined by Eq. (14)
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) I71 -178
77 dimensionless
moving
Eq. (Al) dimensionless dimensionless
8 [
defined by
concentration defined by Eq. (6) spatial coordinate defined by Eq.
(6) T dimensionless cp dimensionless
time spatial coordinate
(Al) dimensionless
r
front position
moving
Eq. (6) !P dimensionless
defined by Eq.
front position
concentration
defined by
defined by Eq. (Al)
177
The next step is the assumption of a profile for Y# which approximates the concentration distribution in the partially extracted region. A suitable quadratic profile, as written below, generally gives excellent results. ?P=a(f)‘+b(t)+c
(A9
By using boundary conditions concentration profile becomes:
(A31 and (A41, the
(A10) Appendix A. Approximate solution by the variational formulation method Using reduced dimensionless S-r cp=-
Dt T=277=S
s
a?P
_-
variables defined as:
(All)
S
In addition, by introducing the assumed quadratic profile, Eq. (AlO), in the last boundary condition, Eq. (A5), the following equation can be obtained:
becomes:
a2!P
a=1K(I-q)g
q[V(r>,
aqo
(‘412)
aq2
TP[o, T] = 1 a!P -=
Eq. (A101 into Eq. (A8) and integrating
S-R
The general set of equations
a7-
Introducing gives:
r] =0
-K(l-n)z
v(O) = 0
atp=q
(‘45) (A61
By taking Fick’s second law as the Euler equation of the desired variational formulation [7], the following equation can be obtained:
(A? Since Eq. (A7) is true for any arbitrary time interval (7,’ TV) the integral relative to the time integration must vanish everywhere in this interval. By assuming a profile for 9, which is specified in the cp direction but depends on the unknown function q(r) in time, Eq. (A7) can be written as: (A81
By substituting Eq. (A12) into Eq. (Al 1) and using the initial condition, Eq. (A6), q versus r can be found. However, since the combination of these equations gives a transcendental differential equation that cannot be solved explicitly, a simplified method with another approximation can be used. Considering “a” to be independent on r in Eq. (A12) and solving this equation with the initial condition (A6) gives:
+2 i
1 77 --3
I
z-7
The combination
(~az-a+l~‘a-l)
(1 -a> K
(Al31
of Eqs. (Al 1) and (A121 gives:
=K(l_77)
(A14)
a(a-3)
The parameters “a” and 77 for each specified value of r can be obtained by solving Eqs. (A131 and (A14).To find fractional release, Eq. (AlO) can be substituted into Eq. (22) and divided by the theoretical amount of releasable solutee per unit volume of
178
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
the sphere given by Eq. (23). The resulting expression for fractional release becomes: ,17-)3
M~
a(r/
-- "r/
'r/--
( 7/-- 1 )
)]
2 -2 - 1
(A15)
To calculate fractional release, find r/ and " a " for each value of ~- should be found by solving Eqs. (A13) and (A14), then substitute these values into
Eq. (A15) to give
Mt/M~.
References [l] J. Crank, Mathematics of Diffusion, 2nd edn., Oxford University Press, Oxford, 1975, pp. 286-325.
[2] H.S. Carslaw and J.C. Jaegar, Conduction of Heat in Solids, 2nd edn., Oxford University Press, Oxford, 1959, pp. 282-296. [3] H. Reiss, J.R. Patel and K.A. Jackson, Approximate analytical solutions of diffusional boundary-value problems by the method of finite zone continuity, J. Appl. Phys., 48(12) (1977) 5274-5278. [4] T. Higuchi, Mechanism of sustained-action medication. Theoretical analysis of solid drugs dispersed in solid matrices, J. Pharm. Sci., 52 (1963) 1145-1149. [5] P.I. Lee, Diffusional release of a solute from a polymeric matrix. Approximate analytical solution, J. Membrane Sci., 7 (1980) 255-275. [6] D.R., Paul and S.K. McSpadden, Diffusional release of a solute from a polymer matrix, J. Membrane Sci., 1 (1976) 33-48. [7] R. Arpaci, Conduction Heat Transfer, Addison-Wesley, Ontario, 1966, pp. 435-477.
Journal of Membrane Science 115 (1996) 171 - 178
Diffusional release of a dispersed solute from a spherical polymer matrix M.J. Abdekhodaie, Y.-L. Cheng * Department of Chemical Engineering and Applied Chemistry, University of Toronto, 200 College Street, Toronto, Ontario, Canada, M5S 1A4 Received 12 July 1995; revised 30 October 1995; accepted 18 December 1995
Abstract The exact analytical solution and an approximate solution based on the variational formulation method have been obtained for the release kinetics of a solute from a spherical polymeric matrix in which the initial solute loading A is above the solubility limit CS. Comparisons between these new solutions and previously reported approximate solutions based on either pseudo-steady state assumptions or the refined integral method indicate that the approximate solutions are very accurate when A / C s is much larger than unity. However, as drug loading decreases, deviations between the approximate solutions and the exact solution increase so that the approximate solutions are least accurate in the limit A / C s ~ 1. Keywords: Controlled release; Diffusion; Theory; Exact solution; Approximate solution
1. Introduction
been used to analyze the release kinetics of a dis-
The diffusional release of a solute from a polymeric matrix in which the initial drug loading A is greater than the solubility limit C s is a problem of interest in controlled release drug delivery. This
persed solute from polymeric matrices. So far two approximate solutions have been developed to analyze the release kinetics of dispersed solute from a polymeric matrix with spherical geometry [4,5]. In both of these solutions, perfect sink conditions were assumed, and matrix swelling of the matrix, the concentration dependence of the solute diffusion coefficient, and the external mass transfer resistance
problem can be formulated by considering a moving diffusional front separating the core containing undissolved drug and the partially extracted region, Moving boundary problems arise in a wide range of application areas, including diffusional release [1], heat transfer involving melting phase transition [2], and diffusion-controlled growth of particles [3]. A1though moving boundary problems have been widely studied, only a few mathematical descriptions have
* Corresponding author. Tel: (416) 978-5500; Fax: (416) 978-
8605; E-mail: y|[email protected],
were neglected. Higuchi [4] obtained an approximate solution based on the pseudo-steady state assumption that at any instant in time, steady state concentration profiles existed in the (time-dependent) partially extracted region. Since the diffusional front velocity should increase as solute loading decreases, the pseudo-steady-state assumption is expected to be least applicable at low solute loadings. Lee [5] used the refined integral method to obtain approximate
0376-7388/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved PII S 0 3 7 6 - 7 3 8 8 ( 9 6 ) 0 0 0 1 8 - X
172
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
analytical solutions for release from planar and spherical matrices containing dispersed drug. T h e approximate solution to the slab problem agreed well with the known exact solution obtained by Paul and McSpadden [6], with the largest deviation at drug loadings near the solubility limit. No exact solution for the spherical problem presently exists for comparisons to be made. Although Lee's solution is elegant and the result is easy to use, the refined integral method of solution is complex. It is therefore the primary objective of this study to develop an exact analytical solution for the release kinetics of dispersed solute from a sphere into an infinite sink. In addition, a simpler procedure for obtaining an approximate solution to this problem based on the variational formulation method will be presented. Comparisons between the exact solution and all known approximate solutions will be made.
Partially
/ExtractedRegionUnextracted / / ~
/
c
o
r
e
\ , . .~. . i ~ R
$
' r
Fig. 1. Schematic diagram of drug distribution in a partially extracted pellet: R, position of moving diffusion front; S, position of sphere surface.
matrix. Assuming perfect sink condition at all times, the initial and boundary conditions are: C[S, t] = 0
(2)
C [ R ( t ) , t] = C s 2. Analysis
(3)
OC
dR
D--=(A-Cs)-The diffusional release of a dispersed solute from a polymeric matrix with spherical geometry into a perfect sink fluid is considered. The diffusion coefficient is assumed to be independent of concentration, and solute diffusion is assumed to be the rate controlling step rather than polymer swelling or drug dissolution. Polymer swelling kinetics can be neglected in systems where little swelling occurs, or where swelling occurs rapidly relative to solute diffusion such that diffusional release essentially occurs from a pre-swollen sphere, The drug distribution in the sphere after a certain time of release is illustrated in Fig. 1 where S is the radius of the entire (possibly pre-swollen) pellet, R is the radius of the unextracted core and r is the polar radial position variable. The concentration is essentially constant in the unextracted core ( r < R), while the concentration in the partially extracted region, R < r < S, is a function of r and is determined by transient diffusion according to Fick's second law: Ot
r 2 Or ~
--~r ]
(1)
where D is the diffusion coefficient, and C is the concentration of dissolved solute concentration in the
at r = R
Or
(4)
dt
R(0) = S
(5)
where A is the initial solute loading per unit volume, and C~ is the solubility limit of the solute in the • matrix. 2.1. Exact solution
By using reduced dimensionless variables defined as: r
Dt
R
~:= ~
~"= -~-
F = -~
r[C~
0 = -~ ~ - ~ )
(6)
the defining equations can be reduced to their dimensionless forms: 00 020 (7) 0r 0~ 2 0[1, ~-] = 0 (8) 0[ F(~-), 'r] = F (9) 1 00
1
F (0) = 1
-
1
at ~ = F
(10)
(11)
173
M.J. Abdekhodaie, Y.-L. Cheng/ Journal of MembraneScience 115 (1996) 171-178
By defining K as: ( A ) K= ~ - 1 Eq. (10) becomes: 1 00 1 dF F 0~ F K d~"
(12)
at ~-- F
(13)
Eq. (7) can be reduced to an ordinary differential equation by using an appropriately defined combined variable U: (14)
1-~:
U=
x/~
such that the dimensionless concentration 0 becomes a total function of the combined variable:
In principle, Eq. (20) can be solved numerically with the initial condition on F g i v e n i n E q . ( l l ) . However, since Eq. (20) gives d I / d r ~ oo at ~"= 0, a value of F at 1- other than zero is needed as the initial condition for the numerical solution of Eq. (20). At very short times (small ~'), the position of the diffusion front is very close to the surface of the sphere and the partially extracted region is very thin. It is well-known that curvature effects are negligible in thin films such that thin cylindrical or spherical films behave similarly to thin planar films. It is therefore reasonable to expect that F(T) at early times for planar and spherical systems should be similar. Paul and McSpadden's [6] planar geometry solution gives: f ~ [ 12@TF]
[(1-/")
2
[ ]1-/"
The diffusion equation, Eq. (7), then becomes:
(21)
(U) d20 dU 2
dO dU
(16)
Integrating directly twice, we obtain: rf( U ] 0 = B t f ~ e ~ -~- ] + B2
(17)
where B l and B2 are constants of integration. By using boundary conditions (8) and (9), 0 becomes: [~] /" erf 0=
and hence: C -- =
Cs
SA M= = - -
1-
(19)
F] ev[ 2x/-T J
Using the last boundary condition, Eq. (13), gives: (1 - / " ) 2
exp
For any combination of solute loading and solubility, Eq. (21) can be used to calculate a value for iV at an arbitrarily small value of ~, (e.g. ~'= 10-6). These corresponding iV and ~- values can then be used as an initial condition (/"0, % ) i n combination with Eq. (20) to numerically solve for the time profile of the dimensionless diffusion front position in a sphere F(~-). The cumulative amount of solute released per unit area at any time t, M t, is related to the diffusional front position R(t) by a simple mass balance:
The total amount of drug, or the amount releasable at infinite time, per unit area, M~,, is given by:
/"erf[~] r
C~
4r
x/__~v/~ erf [_~___r ]
(23)
3
The fractional release at time t in terms of dimensionless quantities then becomes: M--L= (1 - / " 3) M~
1
-~
d /"
3 Cs F
(20)
A
S/~:erf(~_7~)d~:
(24)
M.J. Abdekhodaie, Y.-L Cheng /Journal of Membrane Science 115 (1996) 171-178
174
By substituting F(~-) profiles from the solution of Eq. (20) into Eq. (24) and integrating, the release kinetics can be obtained in terms of the time dependence of fractional release.
.
0.00
. 0.05 .
.
.0.15 .
0.10
0.20
0.25
2.5
3.0
1.0 ~
\
0.8 NCs = 2
2.2. Approximate solution using variational formulation method
0.6-
By taking Fick's second law as the Euler equation of the desired varational formulation and assuming a quadratic concentration profile, the following equations can be obtained to calculate the fractional release (see details in Appendix A): Mt - 1 - r
/
3 -(1
~r
Cs
/
\\ " \\ 0.4
1
~
0.5
1.0
1.5
2.0
(25)
where /" and a for each value of r can be obtained by solving the following equations:
= KF a( a - 3) (l-F):(1 ) (l-a) + F ~ ~-
Ncs=11
0.2
0.0
a
\
I
Ms + ~ - ( / ' + 1)
\
\
F
Fig. 2. F ( ~ ' ) with initial condition F 0 at ~" = ~'o o b t a i n e d f r o m P a u l - M c S p a d d e n p l a n a r g e o m e t r y solution for two values of
A / C s : (a) A / C s = 2 ; - - , ~'0 = 1 0 - 8 to 10 - 2 and . . . . . ro=0.1;(b) A / C s = l l ; - - , ~ - o = l O - 8 to 10 - 2 , - - ~'0=0.1
(26) pected, greater sensitivity was observed for the lower (27)
3. Results and discussion The exact solution for the fractional releases, MJM=, requires knowledge of the dimensionless moving front position F(~-). As described above, at early times, curvature effects can be neglected and F ( z ) from the known exact solution for the planar geometry problem can be adopted for the spherical problem. To use this approach, it is necessary to determine if the solution for _F(~-) is sensitive to the initial conditions selected. Fig. 2 shows the F(~-) profiles obtained from Eq. (20) using varying values of (/'0, %) pairs from Eq. (21) as the initial condition. These calculations were carried out for two drug loadings, A/Cs = 2 and A / C s = 11. The results indicate that for A / C s = 11, indicate that the F(~') profiles are insensitive to the initial condition selected for ~'o ranging from 10 -8 to 0.2. As ex-
drug l°ading ( A / C s = 2 ) " N° differences in I'('r) profiles were observed for ~'0 ranging from 10 -8 to 10 -2. For ~'0 = 0.1, however, a significantly different / ' ( 7 ) profile was calculated. At lower drug loadings, the diffusional front should move inwards at a higher velocity so that curvature effects should become significant at an earlier time. These results confirm the validity of our assumption that curvature effects can be neglected at early times. The results presented below were obtained using / ' values corresponding to ~-= 10 -6 from Eq. (21) as the initial conditions for Eq. (20). Fractional release profiles for various solute loading levels are shown in Fig. 3. The results show a non-linear dependence on r. As expected, as A / C s increases, the cumulative fractional release curves become increasingly sustained. In Fig. 4, fractional release is plotted as a function of z ( C J A ) . Interestingly, the release profiles collapse onto approximately the same curve. This normalized curve can be used to find fractional release profiles for any solute loading. To calculate fractional release based on the ap-
M.J. Abdekhodaie, Y.-L Cheng/Journal of Membrane Science 115 (1996) 171-178 NCs=2
7
11
21
0.8
1.0
1.0
F 0.5
0.5
175
101
NCs=7
0.6
0.0
Mt/M~
0.0
0.4
1.0-
i
~
~
0.1
0.2
0.3
0.0 0.4
[
0.0
0,3
0.6
0.9
1.2
0
1
2
3
4
1.0
0.2 r
0.0 0
I
I
I
1
2
3
0.5
0.5
NCs=I 1 4
0.0
0.0 0
Fig. 3. Fractional release versus ~- for a dispersed solute in a spherical polymer matrix at various solute loading levels based on the exact solution,
1
2
Fig. 5. Diffusion front position vs. time for a dispersed solute in a spherical polymer matrix at various solute loading levels: ,7, Higuchrs solution; (3, Lee's solution; [3, variational formulation solution; z~, exact solution.
proximate variational formulation solution, a direct elimination of F and " a " from Eqs. (25), (26) and (27) is not very convenient. It is more practical to correlate the fractional release with time through
NCs=2
3
101
1.0
0.8
0.6 Mr/M= 0.4
0.2-
o.0o.oo
'
o.04
'
o.o8
' o.12
o.16
x(C~A)
Fig. 4. Fractional release versus ~-(Cs / A ) for a dispersed solute in a spherical polymer matrix at various solute loading levels based on the exact solution,
their dependencies on these parameters. The pseudo-steady-state result given by Higuchi [4] and the solution presented by Lee [5] also require such a correlation through F . The dimensionless time dependence of the diffusion front position, F , and the fractional release, M i M e , are presented in Figs. 5 and 6, respectively. The exact solution and the variational formulation approximate solution presented in this work are shown along with Higuchi's and L e e ' s aproximate solutions. To compare these solutions, the deviation between the fractional release based on each of the approximate solutions and the exact solution at equally spaced values of r throughout the release duration were calculated, then averaged. The timeaveraged deviation at various drug loading levels are shown in Table 1. Increasing A / C , results in reduced deviations. For example, the deviations of the Higuchi, Lee and variational formulation approximate solutions from the exact solution for A / C ~ = 3 are 7.15%, 3.91% and 5.14%, respectively. The corresponding deviations for A / C ~ = 2 1 decline to 0.91%, 0.34% and 0.54%. In all three approximate solutions, an assumption had to be made regarding the concentration profile of the solute in the partially
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
176
1.o f
1.o f
M,/M~0.5
0.5 A/Cs=3
,
,
,
0.1
0.2
0.3
0.0 0.0
NCs=7
0.0 ¢ 0.4
1.0 ~
0.0
f
0.S NCs=I 1
0.0
i
0
0.3
, 0.6 •
, 0.9
1.2
1.0 ~
/
~/M~ 0.5
,
1
0.0
2
0
NCs=21
I
i
1
2
3
Fig. 6. Fractional release vs. time for a dispersed solute in a spherical polymer matrix at various solute loading levels: V, Higuchi's solution; ©, Lee's solution; [2, variational formulation
solution; z~, exact solution,
extracted region. In the Higuchi solution, a pseudo steady state profile was adopted, while a quadratic profile was used by Lee and in the present work. A parameter " a " , given in the appendix by Eq. (A12) was required in the quadratic profiles. This parameter was used to match the assumed profile to the position of the moving diffusion front. Although this parameter should vary with time, it was assumed to be time-independent in both approximate solutions. Since the moving front velocity increases as the solute loading decreases, the pseudo-steady-state assumption used in Higuchi's approximate solution, as Table 1 Comparison of average absolute percentage deviations of different approximate solutions from exact solution for a spherical polymer
matrix A/C S
Higuchi error (%)
well as the assumption of time-independence for the parameter " a " in the quadratic profiles, are less valid at low solute loadings, leading to greater inaccuracies at low loadings. Of the three approximate solutions, Lee's refined integral approximate solution was closest to the exact solution. The accuracy of refined integral and variational formulation solutions depends on the form of the assumed concentration profile. The assumption of a quadratic concentration profile resulted in a more accurate solution using the refined integral method for this particular problem. More accurate solutions can be achieved with improved assumed concentration profiles.
Lee
Variational
error (%)
formulation
solution
error (%)
3
7.1521
3.9082
5.1418
7 11 21
3.6999 1.6486 0.9072
0.7815 0.6361 0.3393
1.7284 1,0131 0,5379
4.
Conclusion
An exact moving boundary solution has been presented for the diffusional release of dispersed solutes from a spherical matrix into an infinite sink medium. This exact solution permits the validation of previously available approximate solutions presented by Higuchi [4] and Lee [5]. In addition, an alternative approximate solution based on a simple variation formulation method is presented. The approximate solutions have some advantages for describing release kinetics due to their relative simplicity. These solutions are very accurate when A / C s is in great excess of one, but in the limit of A ~ Cs, deviations from the exact solution increase.
5. List of symbols A initial solute loading per unit volume Bi constant of integration C solute concentration in the matrix at any time C s equilibrium solute solubility in the m a t r i x D solute diffusion coefficient in the matrix K initial solute loafing level defined by Eq. (12) M t amount of solute release per unit area M ~ total amount of solute release per unit area r spatial coordinate R position of the moving diffusion front S radius of the sphere t
release time
U combined variable defined by Eq. (14)
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) I71 -178
77 dimensionless
moving
Eq. (Al) dimensionless dimensionless
8 [
defined by
concentration defined by Eq. (6) spatial coordinate defined by Eq.
(6) T dimensionless cp dimensionless
time spatial coordinate
(Al) dimensionless
r
front position
moving
Eq. (6) !P dimensionless
defined by Eq.
front position
concentration
defined by
defined by Eq. (Al)
177
The next step is the assumption of a profile for Y# which approximates the concentration distribution in the partially extracted region. A suitable quadratic profile, as written below, generally gives excellent results. ?P=a(f)‘+b(t)+c
(A9
By using boundary conditions concentration profile becomes:
(A31 and (A41, the
(A10) Appendix A. Approximate solution by the variational formulation method Using reduced dimensionless S-r cp=-
Dt T=277=S
s
a?P
_-
variables defined as:
(All)
S
In addition, by introducing the assumed quadratic profile, Eq. (AlO), in the last boundary condition, Eq. (A5), the following equation can be obtained:
becomes:
a2!P
a=1K(I-q)g
q[V(r>,
aqo
(‘412)
aq2
TP[o, T] = 1 a!P -=
Eq. (A101 into Eq. (A8) and integrating
S-R
The general set of equations
a7-
Introducing gives:
r] =0
-K(l-n)z
v(O) = 0
atp=q
(‘45) (A61
By taking Fick’s second law as the Euler equation of the desired variational formulation [7], the following equation can be obtained:
(A? Since Eq. (A7) is true for any arbitrary time interval (7,’ TV) the integral relative to the time integration must vanish everywhere in this interval. By assuming a profile for 9, which is specified in the cp direction but depends on the unknown function q(r) in time, Eq. (A7) can be written as: (A81
By substituting Eq. (A12) into Eq. (Al 1) and using the initial condition, Eq. (A6), q versus r can be found. However, since the combination of these equations gives a transcendental differential equation that cannot be solved explicitly, a simplified method with another approximation can be used. Considering “a” to be independent on r in Eq. (A12) and solving this equation with the initial condition (A6) gives:
+2 i
1 77 --3
I
z-7
The combination
(~az-a+l~‘a-l)
(1 -a> K
(Al31
of Eqs. (Al 1) and (A121 gives:
=K(l_77)
(A14)
a(a-3)
The parameters “a” and 77 for each specified value of r can be obtained by solving Eqs. (A131 and (A14).To find fractional release, Eq. (AlO) can be substituted into Eq. (22) and divided by the theoretical amount of releasable solutee per unit volume of
178
M.J. Abdekhodaie, Y.-L. Cheng / Journal of Membrane Science 115 (1996) 171-178
the sphere given by Eq. (23). The resulting expression for fractional release becomes: ,17-)3
M~
a(r/
-- "r/
'r/--
( 7/-- 1 )
)]
2 -2 - 1
(A15)
To calculate fractional release, find r/ and " a " for each value of ~- should be found by solving Eqs. (A13) and (A14), then substitute these values into
Eq. (A15) to give
Mt/M~.
References [l] J. Crank, Mathematics of Diffusion, 2nd edn., Oxford University Press, Oxford, 1975, pp. 286-325.
[2] H.S. Carslaw and J.C. Jaegar, Conduction of Heat in Solids, 2nd edn., Oxford University Press, Oxford, 1959, pp. 282-296. [3] H. Reiss, J.R. Patel and K.A. Jackson, Approximate analytical solutions of diffusional boundary-value problems by the method of finite zone continuity, J. Appl. Phys., 48(12) (1977) 5274-5278. [4] T. Higuchi, Mechanism of sustained-action medication. Theoretical analysis of solid drugs dispersed in solid matrices, J. Pharm. Sci., 52 (1963) 1145-1149. [5] P.I. Lee, Diffusional release of a solute from a polymeric matrix. Approximate analytical solution, J. Membrane Sci., 7 (1980) 255-275. [6] D.R., Paul and S.K. McSpadden, Diffusional release of a solute from a polymer matrix, J. Membrane Sci., 1 (1976) 33-48. [7] R. Arpaci, Conduction Heat Transfer, Addison-Wesley, Ontario, 1966, pp. 435-477.