- Email: [email protected]
Analysis of controlled release in disordered structures: a percolation model
> ~ ~ j o u r n a l of MEMBRANE SCIENCE
v
ELSEVIER
Journal of Membrane Science 113 (1996) 21-30
Analysis of controlled release in disordered structures" a percolation model Alessandra Adrover a,,, Massimiliano Giona b Mario Grassi c '~Dipartimento di lngegneria Chimica, Universitil di Roma "La Sapien=a", Via Eudossiana 18, 00184 Roma, Itah' h Dipartimento di lngegneria Chimica, Universitd* di Cagliari, Piaz=a d'Armi, 09123 Cagliari, ltal3 • Dipartimento di lngegneria Chimiea, Universitd di Trieste, Pia:=a Europa I, 34127 Trieste, Italy Accepted 9 August 1995
Abstract
The analysis of controlled-release kinetics in disordered models (percolation clusters) of polymer matrices is developed in detail. It is shown that the "anomalous diffusion" experimentally found in the release kinetics from swollen gels, M , / M ~ ~ t", n > 1 / 2 , cannot possibly be related to hindered diffusive motion in disordered media, i.e. is not a consequence of the disordered structure of the matrix. A kinetic model is proposed to account for an exponent n greater than 1/2. This is based on a two-phase kinetics which takes into account both field effects (deriving from potential interaction between solute and polymeric matrix) and entrapping effects due to geometric constraints. Kevwords. Release kinetics: Percolation: Anomalous transport: Two-phase models
1. Introduction In recent years much work has been focused on the development of controlled release devices [ 1 - 4 ] . These drug-delivery systems provide release rates which are determined by the physical properties of the device itself. There exist different kinds of drug delivery systems, amongst which polymeric matrices are widely used. The rate-determining step to produce a specific dosage form is represented by solute diffusion inside the polymeric matrix, depending on the structural properties of the gel. Gel structures can be subdivided into two classes: chemical gels, which form a lattice that is stable also in the presence of a swelling agent; and physical gels, which may undergo structural deformation in time (dissolution). * Corresponding author. 0376-7388/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved S S D 1 0 3 7 6 - 7 3 8 8 ( 95 ) 0 0 2 2 0 - 0
Two different drug-release mechanisms occur in hydrophilic cross-linked polymer systems depending on the macromolecular state of the system, which can be in a glassy and rubbery state during the initial stage of the swelling agent diffusion or in a rubbery state after thermodynamic equilibrium is attained between the network and the surrounding dissolution medium ( usually water) [5]. In the study of swelling controlled systems, the overall solute release depends on the thickness of the swollen layer and on the velocity of the swelling interface. In this situation, the polymeric matrix structure changes in time and solute release develops simultaneously with the propagation of the swelling interface [6,7]. In the case of solute diffusion through swollen gels, the mechanism of transport depends on the relative size of the diffusing species, on the mesh formed by the macromolecular chains [8,9] and on the potential effects related to solute/polymeric
22
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
interactions. Morphological features associated with reduced macromolecular chain mobility and potential barriers to solute diffusion can be responsible of the decrease in the solute diffusivity. Assuming that field effects are negligible, the diffusion coefficient D of the solute molecules in liquid-filled porous matrices deviates from the corresponding coefficient in the dilute bulk solutions DE (Einstein diffusivity) by a factor G, called the accessibility factor, which depends exclusively on steric properties (size exclusion effect, SEE) [10]
D= DEG,
(1)
An approach to this kind of problem can be based on the simulation of transport in disordered structures representing the polymeric matrix: transport parameters are evaluated from reasonable models of the porenetwork [ 11,12]. Besides the hindered effects related to the porosity and to the crosslinking of the matrix, the key problem in the kinetics of controlled release from swollen gels is represented by the anomalous release behaviour from polymeric devices, which generally show, for intermediate time scales, an empiric power-law behaviour [6,13,14],
M,/M~~f', M,/M~ <0.5,
(2)
with an exponent n usually greater than 1/2 (pseudo non-Fickian diffusion). M,/M= is the fraction of solute released up to time t. A value of n > 1/2 is contrary to the classical release models of Fickian diffusion [ 3 ]. In order to understand the origin of the "anomalous" exponent n, we analyze the release kinetics from disordered models of polymeric matrices by taking percolation clusters as the paradigm of a disordered structure [ 15,16]. Percolation lattices are simple and general models of porous structures with spatially uncorrelated disorder [16]. We focus attention on swollen or non-swellable devices in the form of slabs and develop a percolation model of this kind of drugdelivery system in order to evaluate the influence of the complex and disordered nature of the polymeric matrices on the drug-diffusion coefficient and on the exponent n. The cases of lag and burst operating conditions [ 13 ] are analyzed. We show by means of numerical simulations that the transport mechanism in disordered structures is not responsible for the anomalies found in sorption kinet-
ics. In order to explain the apparent power-law release (Eq. (2) with an exponent n greater that 1/2), we propose a simple kinetic model which considers the entrapping effects of solute molecules inside the polymeric matrix, deriving from field interactions or geometric hindrances, in the form of a two-phase kinetics. The numerical results on percolation lattices prove to be in agreement with Eq. (2). The article is organized as follows. First, we discuss the simulation method for analyzing transport phenomena in disordered lattice models. Second, we analyze two typical release conditions (lag and burst effect) in percolation structures in order to evaluate G and n. We then present the two-phase kinetic model, showing that the entrapping kinetics may be responsible for the anomalous behaviour represented by Eq. (2). Finally, we analyze the differences between anomalous transport in fractal and disordered structures (in the sense discussed by Havlin and Ben-Avraham [ 17] ) and the notion of anomalous transport in controlled release systems in order to clarify this problem, which is not only a matter of semantics.
2. Percolation models and numerical simulations A percolation lattice is a spatially uncorrelated random distribution of empty and occupied sites characterized by a percolation probability p representing the frequency of empty sites (site percolation). In the application of percolation models to transport in polymeric matrices, the occupied sites are the elements of the matrix and the empty sites are the pores through which solute particles diffuse. From percolation theory [ 16], it is known that there exists a critical percolation probability, Pc, above which transport is allowed. For values of P>Pc there exists a connected cluster of empty sites (called the infinite cluster) all over the lattice. In two dimensional square lattices Pc = 0.593. For each value of p greater than Pc, in addition to the infinite cluster there exists a family of finite subclusters of various lattice sizes. In order to generate the infinite connected cluster alone, Leath's growth algorithm (spreading percolation) can be adopted [18]. As extensively discussed by many authors [ 15 ], percolation lattices constitute a "model for all seasons" to represent a disordered structure with uncorrelated spatial disorder.
A. Adrover et al. /Journal of Membrane Science 113 (1996) 21 30 The classical approach to the analysis of transport in disordered media and fractals is based on Monte Carlo simulation of random walks, which is grounded on the equivalence between Brownian motion and diffusion [ 19]. The random-walk approach is particularly suitable for analysis of the statistical properties associated with transport (e.g. the scaling of the mean square displacement with time) [ 17] but is of less utility if we want to solve specific transport boundary-value problems, as in the case of sorption kinetic simulations. For this reason, we propose a simple and somewhat elementary finite-difference numerical method for parabolic and elliptic differential equations in disordered lattice models. This satisfies the conservation principle [20] and furnishes promising results in the analysis of linear and nonlinear transport models on fractals and on disordered structures [21,22]. Let us consider the problem of diffusion in two dimensions,
Oc ~-[O~'c 02c~
(3)
where c is the concentration of the transported entity. A generic explicit algorithm for solving Eq. (3) takes the form (n+
%
I ~ --
.(n)
-~i
+fl'
~
jo,)
hk~j,
(4)
h,h } ~ I( i,]
where ¢3' = ~ A t / A x and JJ,'2Lo is the local flux at time in) n, Jhk~"~,j= (~hk'~"~--Gj)/Ax.. Let M be a disordered medium on a square lattice, (id) the lattice coordinates of a generic point and l ( i j ) the set of the first four nearest neighbours of (id). The medium M can be represented by means of its characteristic function
c o( n + l ) =c~j( n )
_~_ [3 ¢
~_,
23
(n} xIKk)J),k~,
(h,k)El(i,j)
'"'+[3
= ¢_!]
E
x(h,k)[d,i"-,,, I
~6)
{ h,k ~ E l ( i , . ] )
where [3 = DA t / A x z. It is important to observe that Eq. (6) satisfies the conservation principle, i.e.
~x(i,j) (i,j}
~
x(h,k)[Gk-c,] =0,
(7}
(It,k) ~ : l ( i , j )
which is a key condition lbr the applicability of a finitedifference scheme deriving from mass conservation. Of course, from the stability analysis of explicit schemes,/3 < 0.5. It should also be noted that the same approach makes it possible to develop more sophisticated numerical techniques (e.g. implicit alternating direction) in a straightforward way. Moreover it is easy to extend the same numerical technique to the case of multiphase transport problems (see section 4).
3. Lag and burst effects: evaluation of the accessibility factor Let us consider the controlled release of a solute from a slab of a swollen polymeric matrix of thickness L, total volume Vo, and void fraction e. The drug is released into a recovery volume V. Drug concentration in the slab is initially equal to zero and in x = 0 the slab is placed in contact with a solution in which drug concentration is kept constant and equal to co (infinite reservoir volume). In this case the polymeric matrix acts as a membrane. In a one-dimensional approach, the balance equation reads as
&,
32c
X( i,j )
X( i,j ) =
l, ( i , j ) ~ M 0, otherwise.
(5)
The initial and boundary conditions describing the lageffect experiment are [ 13 ] : c(x,0) = 0 , 0 < x < L
Eq. (4) is still valid in a disordered structure if care is taken that Eq. (4) is defined only in a point (id') ~ M , (n) i.e. for x ( i J ) = 1, and the local fluxes l~'hk~j are identically equal to zero if x ( i j ) = 0. According to these observations, the simplest explicit algorithm for the analysis of transport in a disordered lattice structure reads
[3c
(9)
hD Oc)
L G=,"
10)
where h = Voe/V. The latter boundary condition takes into account the fact that the recovery volume V is not infinite and consequently stationary conditions will be asymptotically reached leading to a finite value of M-~.
A. Adroveret al. / Journalof MembraneScience 113 (1996)21-30
24
Mt/Moo
From the solution of the differential scheme ( 8 ) - (10), the expression of the fraction of solute M,/M= released up to time t attains the form: y ' B,,exp(
Mr~M== 1
m=l
-
0.75
A2mDt)
~
a
(11) 0.50
~_Bm m=l
/ /
0.25
(A~L2 +h2)c°S(AmL) h,,COS(hmL)=h/L. Bin- L[(A2mL2+h2) + h ] '
a) p-o.95 b) p.o.a5 c) p-0.75 d) p-0.65
0.00 0
(12) For t tending to infinity the leading term in Eq. ( 11 ) corresponds to the smallest eigenvalue hj (fundamental mode) and consequently: (1 -mr~m=) ~ e x p ( - h 2 1 D t ) , t>> 1/(DA2).
(13)
The slope of the curve In (1 -M,/M=) vs. t therefore enables us to estimate the effective diffusion coefficient D. Eq. ( ! 3) is the typical exponential relaxation associated with diffusive motion. Simulation of the lag effect on percolation lattices can be simply performed by following the approach outlined in section 2 and imposing the boundary and initial conditions ( 9 ) - ( 1 0 ) on the lattice. Simulations were performed on a two-dimensional square lattice of size 5 0 × 150 for different values of the percolation probability p ( e = p ) , p > Pc and by choosing h = 0. lp. The size of the lattice is large enough for us to neglect finite-size effects. The choice of a larger size in the direction orthogonal to transport makes it possible to treat the phenomenon as one-dimensional by evaluating the concentration as the average over the direction orthogonal to transport. Fig. 1 shows the resulting sorption curves for different values of the percolation probability p. If s(p,h) is the slope of the curve ( 1 - M , / M = ) vs. t in a lognormal plot (see Fig. 2), then for large t we have
200
400
t
Fig. 1. M,/M~ vs. t for the lag-effect simulations [Eqs. (8)-( 10)] for different values of the percolation probabilityp: L = 1 a.u., D = 1 a.u..
taking into account its internal porosity. Conversely, we can define the accessibility factor GE(p) ( " E " stands for Euclidean in contrast to " p " , percolative) by taking the slab volume Vo as the reference volume and therefore by considering for the factor h the expression h = h E - Vo/V, since the effects of the porosity are encompassed in the accessibility factor. The relationship between G" and G E is therefore given by GE(p) = GP(p)p.
(15)
Fig. 3 shows the comparison of the values of GE(p) obtained by applying Eq. (14) to sorption kinetic simulations with the values of GE(p) obtained from Monte Carlo simulations [ 12]. The figure shows a fairly good level of agreement. It should also be observed that while s(p,h) and Al are functions o f p and h, the accessibility factor is only a function of p, i.e. it depends exclusively on the steric properties of percolation lattices and does not depend on the boundary condition. The burst con1-( Mr/M** ) 10 o ~
~
al °-O.g5 b) p-o.a5
E\\\
c,o.o, d,°-o°,
s(p,h) GP(p) =
AE(h,L ),
(14)
where GP(p) is the value of the accessibility factor for a fixed value p of the percolation probability evaluated by taking as reference volume the effective volume of the slab Vo'p in which the solute diffuses. The superscript " p " indicates simply that we assume the effective volume of the slab as the reference volume by
1o0[
\ 100\
200 t 300 Fig. 2. Log-normal plot of 1- M,IM=vs. t for the data of Fig. 1.
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30 GE
2oMC
10
~
burst effect
•
,oo,-,c,
05 i
i i
0.0
.
.
.
.
.
s(p) rfl/4L2.
GP(p) -
0 ql
.
.
.
00
.
.
.
.
.
.
.
0 5 (p-pc)/(1-Pc) 10
Fig 3 Comparison of the values of GE(p) obtained from Monte Carlo simulation [ 12] with the values of GE(p) for lag and bursteffect s i m u l a t i o n .
dition represents another physical situation often encountered in controlled release systems: the matrix is initially loaded up to a given solute concentration Co. This observation is confirmed by simulation results of controlled release in the burst condition [13], The transport equation describing this phenomenon is still given by Eq. (8), to be solved with the following initial and boundary conditions:
c(x,O)=co, O
(20)
In the case of burst-effect conditions, the accessibility factor is evaluated with respect to the effective slab volume allowed to the solute diffusion, i.e. with the above notations it is given by GP(p). The corresponding values of GE(p) are shown in Fig. 3, in agreement with the lag-effect simulation experiments. The second important question that can be addressed by means of percolative models is the evaluation of the exponent n appearing in Eq. (2) from sorption experiments on disordered lattice structures. By considering, for example, the burst effect (a similar result can be obtained from lag conditions), simulations on percolation lattices give a value of n equal to 0.5 for p at a distance from the critical threshold. This is evident from the slope of the curves ofM,/M~ ( forM,/Ms < 0.5 ) plotted 1-( Mr/Moo ) 10° %
\'\
-
\\
~\ '\
(16) (17)
The latter boundary condition comes from the assumption of perfect sink conditions in x = L. The solute fraction released up to time t attains the form
a
~
'
"\ e
0~xC, _ o = 0 , c(L,t) =0.
25
a) p-100 b) p-0.90 d
'c
c) p-080
e) p-065
d) 1:}-0.70
10-2 0 1 2 3 4 t 5 Fig 4 Log-normal plot of I -M,/M~ vs t for the bursteffect sim ulations for different values of the percolation probability p: L = I a u,D=l a u oMt/M=
M,/M=
1
8 77"
_
=
1 I
(2m
-
1 )2
exp
( 2m--1)27r2Dt] 4L=
]
(18)
'° I 1O- 1
and
In( 1 - M,/M~) ~ - ~ _ D t .
a) p-l.00 b) p-0.90 c) p-0,80 d) p-070 e) p-0.65
(19)
Again, if s(p) is the slope of the curve ( 1 - M , / M = ) vs t in a log-normal plot (see Fig. 4), for large values of t we have
/~/~//~
102
........ ................................. 10 -4 10 -a 10 "a 10 -~ 100 t 101 Fig 5 Log-log plot of M,/M= vs t for a burst-effect simulation for different values of the percolation probability p
26
A. Adroveret al. / Journal of MembraneScience 113 (1996) 21-30
in log-log scale, Fig. 5. Only forp = 0.65 does the slope change slightly and attain the apparent value of n = 0.38 _+0.03. This result can be explained by considering that for p tending toward criticality (Pc = 0.593) the highly heterogeneous nature of the porous matrix emerges, as appears from anomalous characteristics in diffusive motion (n4: 1/2). However, the resulting exponent n is always less than 1/2. An analysis of this topic is developed in section 5. To sum up, the apparent deviation from Fickian behaviour experimentally observed in many controlled-release devices cannot be explained by means of anomalous transport effects in disordered matrices. Consequently, the physical origins of the anomalies in M,/M=-kinetics should be sought in other controlling mechanisms.
4. Two-phase kinetic model The above results indicate that an exponent n greater than 1/2 is not related to anomalous diffusion effects in disordered media. Deviation from n = 1/2 can be due to different phenomena: concentration dependentdiffusion coefficient, erosion of the gel matrix, convective contribution to mass transport, etc. In the case of non-erodible matrices, under the assumption of a practical concentration-independent diffusion coefficient and negligible convective effects, a simple way of explaining the phenomenological Eq. (2) is to regard it as a consequence of interaction effects between solute and polymeric matrix. This can be modeled by assuming the presence, in the release matrix, of two phases: a gel phase ("entrapping phase") and a sol phase in the liquid solution. Two-phase models in the analysis of controlled release have been considered by other authors [ 13,23 ] in order to describe the effects of interfacial transfer kinetics on release rates for those systems which show significant free energy barriers for transport across liquid-liquid interfaces [ 13 ] or to model the release-condition above the solubility threshold [ 23 ]. In the proposed model we assume that a fraction es of solute molecules is initially entrapped in the gel phase and irreversibly released in the sol phase in which it diffuses with an effective diffusion coefficient D. The phenomenon is close to an adsorption and irreversible desorption first-order kinetics with a kinetic coefficient
k. The dynamic equations describing the phenomenon are: Oc 02c Oq o t = Dox 2 0 t
(21)
Oq --= -kq, Ot
(22)
where c(x,t) and q(x,t) are drug concentrations in the gel and sol phases respectively. Eqs. (21) and (22) are to be solved with the initial and boundary conditions (burst effect): c(x,O) -----(1--E~)coO
(23)
0c
(24)
Xx=
=0, c(L,t) =0. 0
In the simulations on percolation lattices, the stationary phase is represented by the first nearest neighbouring site of each site representing the polymeric matrix (see Fig. 6) and is therefore implicitly a function of position since it is defined at the matrix boundary sites. Fig. 7 shows the simulation results of the two-phase kinetics on the percolation cluster for different values of the kinetic constant k ( in arbitrary units), keeping constant the values of the initial partition coefficient es and of the percolation probability p (es = 1, p = 0.80). Fig. 8 shows the behaviour of the exponent n as a function of k for different values of p. For high values of k the internal release kinetics is so fast that the presence of
I
polymericmatrix~
sites
matrixboundary
sites
Fig. 6. Schematicrepresentationof two-phasesimulation in percolation lattices. In the simulationsthe concentrationq is definedonly on the nearest neighbors of a matrix site (matrix boundary site) while c is defined in all the sites not belonging to the polymeric matrix.
A. Adrover et al. /Journal of Membrane Science 113 (1996) 21-30
27 ]
100 Mt/M.__~
n /
1.2 !m ,i /
0.8-
10-; f bc
[°
~' ,//
.I k.5o0 b) k-2oo
e ~
c) k-lO0
d) k-50
e) k-lO
f) k-1 ,
¢
i m_~ 0.4 i
!
102L .... , ,,, , , , ,, , 10-2 10-1 10 o t 101 Fig. 7. Log-log plot of M,/M= vs. t for the two-phase percolation model for different values of k (a.u.) and fixed values o f p and e,
i
0.0
(p = 0.80, e+= 1).
~
0.00
0.25
.
.
.
0.50
.
0.75
1.00
6~
Fig. 10. Plot of n vs. ~+ for a fixed value of p and k (p = 0.80, k = 1 1.2
a.u.). a
\
n
1.0
c
0.8~ 0.6 0.4 0.2
a) p-0.90 b) p-0.80 c) p-0.70 d) p-0.65 10 -1
10 o
101
102
10 a k
104
Fig. 8. Plot o f n vs. k (a.u.) for a fixed value o f p and ~+ (p=0.80, am=l).
100 F
Mt/M°°
i
10 -1 ~s.O. 1 i i
10 -2 !
//
/
r-
/
[ 10-3~
/
/ /'
rC+.lO ....
i
F
'
10 -+ . . . . . . . . . . . . . . . . . . . . . . . 10-3 10"2 10-1 10° t
,ll
101
Fig. 9. Log-log plot of M,/M= vs. t for the two-phase percolation model for different values of e, and fixed values ofp and k (p = 0.80, k = l a.u.).
the gel phase is practically negligible and n tends to the Fickian value 0.5. The influence of the initial partition coefficient es of the solute molecules between the two phases is evident from Fig. 9, which shows the simulation results of the two-phase model for different values of a,, keeping constant the percolation probability and the kinetic constant (p = 0.80, k = 1 a.u. ). Evaluation of the slopes of the sorption curves of Fig. 9 for intermediate time scales reveals that n is a monotonically increasing function of es (see Fig. 10). Therefore, we may conclude that a simple kinetics model which takes into account the internal release kinetics between a sol and a gel phase is able to predict a release exponent n greater than 1/2. As the simulations show, the value of n depends on the kinetic parameters introduced in the model. It should be noted in passing that the predictions of the model proposed have been checked against experimental results on controlled release kinetics by developing a detailed data analysis. In [25] it is shown that the predictions of the model (21 ) - ( 2 2 ) are in agreement with experimental results on the release rate of theophylline from water-swollen scleroglucan matrices.
5. " A n o m a l o u s transport" in fractal m e d i a and in controlled release s y s t e m s
The empirical relation (2) is usually regarded as the deviation from Fickian behaviour in controlled release
28
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
systems. Since deviations from Fick's law are also observed in diffusion through disordered fractal structures, and since polymeric matrices can be regarded as disordered media, it is useful to analyze the differences between these two "anomalous transport phenomena" critically in the light of the results obtained in the previous sections. In the theory of dynamic phenomena on fractals, anomalous transport means deviation from Fick's law in the specific sense that the mean square displacement of a Brownian particle, (r2(t)), deviates from a linear behaviour with time and follows a scaling law,
(r2(t)) - t ~,
(25)
with an exponent/3 less than unity. Eq. (25) is valid asymptotically, i.e. for t ~ w. For some classes of deterministic fractals it is possible to obtain explicit expressions for the exponent/3 by means of renormalization group techniques (decimation methods) [ 17]. The propagation of a diffusion front in a fractal medium does not display Gaussian behaviour but has a stretched Gaussian form. For more details on all these topics see [ 17,24]. In Eq. (25) /3 is a universal exponent characteristic of the fractal structure of the medium. Eq. (25) implies that the heterogeneity at all the length scales of a fractal structure induces a characteristic slowing down (/3 < 1 ) in the diffusional propagation for all the time scales. This effect can be interpreted as the appearance of correlations (for uncorrelated processes/3 is strictly equal to unity) induced exclusively by geometric and topological disorder, although particle dynamics at microscopic scales is still Brownian. Anomalous transport in swollen polymers (the case of simultaneously swelling front propagation is not considered here since the phenomenology is slightly different) also shows an exponent n which deviates from the classical value n = 1/2 (which is strictly the counterpart of/3 = 1 for the scaling of the mean square displacement), but the physical origin of this deviation is completely different from the case of diffusion in disordered or fractal media. The scaling law (2) is not valid asymptotically, as in the case of Eq. (2), and is not a manifestation of universal physical properties, as in the case of the exponent/3. The simulations on percolation lattices developed in the previous sections have shown that for these disordered structures near criticality the exponent n in pure diffusion is definitely less than 1/2, as should be
expected since disorder tends to slow down diffusive motion. Consequently, an "anomalous" exponent n greater than 1/ 2 (which implies an enhancement of the release) cannot possibly be related to non-Fickian diffusion as in the case of fractal structures. The phenomenological origin of this exponent should be sought in the presence of correlations in the diffusive motion deriving from entrapping or field effects (which ultimately have the same physical origin). For these reasons, we have introduced in section 5 the simplest physical model consistent with these assumptions, and the simulation results confirm that such a model can predict an apparent scaling exponent n greater than 1 / 2. The simulation results confirm that the exponent n is not universal but depends on the kinetic constant k. To sum up, the appearance of a pseudo "non-Fickian" exponent n > 1/2 can be regarded as the superposition of regular Fickian transport and of kinetic effects described macroscopically with simple twophase model, but whose physical origin depends ultimately on the specific interactions between solute molecules and polymeric matrices.
6. Conclusions and discussion The expression "anomalous transport" usually refers to two different transport situations: the pseudo non-Fickian behaviour in controlled release kinetics (appearing as an exponent n > 1 / 2 in the sorption curve scaling with time) and the theory of transport in disordered and fractal structures, which are characterized by an exponent/3 less than unity in the scaling of the mean square displacement with time. It is intuitive to suppose that the pseudo non-Fickian behaviour observed in many controlled release experiments could be attributed to the disordered structure of the polymeric matrix. The simulations of the sorption kinetics in disordered models (represented by percolation clusters) in the case of both lag and burst conditions have clearly shown that an exponent n greater than 1/2 cannot be related to the disordered geometric and topological structure of the matrix. Indeed, in the case of percolation clusters near criticality, we observed an exponent n less than 1/2, as is to be expected from the scaling theory of transport in fractal structures [ 17,24].
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
This result is fairly important in the case of drug diffusion through swollen polymeric matrices (since in this case the anomalies cannot be attributed to the simultaneous propagation of the gel front) because it indicates that the origin of the apparent non-Fickian behaviour should be attributed to other kinetic mechanisms rather than to the microscopic disordered structure. We therefore propose a very simple kinetic model by including the effects associated with the solute/ polymeric matrix interactions. In macroscopic models, the interactions can be described by means of a twophase kinetic model, as developed in section 4. It is important to observe that the entrapping effects expressed macroscopically by means of two-phase models may be related to different physical mechanisms (see e.g. the different approaches in the analysis of usual adsorption phenomena): the presence of potential energy fields related to solute/matrix interactions, the entrapping effects related to hindered constraints and the physical adsorption of the solute molecules on the polymeric matrix. The analysis of these important aspects involves detailed microscopic modeling of the specific solute/matrix interactions, which goes beyond the scope of investigations of macroscopic models. However, in a macroscopic perspective it is important to observe that the simple kinetic model we propose is able to account for the fundamental features of the apparent non-Fickian release. Moreover, experimental results on the release of theophylline from water swollen scleroglucan matrices [25] confirm the results predicted by the two-phase model presented here. A final comment is related to the simulation method adopted in dealing with transport in disordered lattice models. The numerical method discussed, which has furnished good results in the simulation of Brownian motion on fractals [ 20], represents an improvement on Monte Carlo simulation and can be applied to a great variety of transport phenomena of interest in controlled-release modelling (e.g. the analysis of potential fields on release kinetics from polymeric materials). Of course, if the structural parameters of the polymeric matrix are known, this technique can be applied to a microscopic lattice model of the polymeric structure by constructing a lattice simulator of the matrix through the reconstruction method of porous solids proposed by Adler et al. [26].
29
7. List of symbols Bm e Co
D G g~ Gp h J k L M M,/M~ n
P p¢. q
(r~(t)) s(p,h ) t V
Vo X
see Eq. ( 11 ) drug concentration initial drug concentration effective diffusion coefficient see Eq. 3 Euclidean accessibility factor percolative accessibility factor
Voe/V local flux kinetic coefficient gel slab thickness disordered medium fractional solute release anomalous exponent [ see Eq. ( 2 ) ] percolation probability percolation threshold drug concentration in the gel phase mean square displacement see Eq. (14) time recovery volume gel slab total volume spatial coordinate
7.1. Greek letters
[3
universal exponent characteristics of a fractal medium DAt/Ax
A1
X( ij)
void fraction initial drug partition coefficient between the sol and the gel phase fundamental mode [see Eq. (13) ] characteristic function of the porous medium M
Acknowledgements We should like to thank A.R. Giona for constant encouragement and useful discussions. One of the authors (M.G.) would like to thank M. Liauw for useful discussions.
30
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
References [ 1] R.W. Baker and H.K. Lonsdale, Controlled release: mechanism and rates, in A.C. Tanquary and R.E. Lacey (Eds.), Controlled Release of Biologically Active Agents, Plenum Press, New York, NY, 1974, p 73. [2] R.S. Langer and N.A. Peppas, Chemical and physical structure of polymers as carriers for controlled release of bioactive agents: A review, J. Macrom. Sci., Rev. Macromol. Chem. Phys., C23 (1983) 61. [ 3 ] J. Crank, Mathematics of Diffusion, 2nd ed., Clarendon Press, Oxford, 1975, [4] L.S. Langer and D.L. Wise (Eds.), Medical Application of Controlled Release: Application and Evaluations, Vol. I-II, CRC Press, Boca Raton, FL, 1984. [5] P.I. Lee, Kinetics of drug release from hydrogel matrices, J. Controlled Release, 2 (1985) 277. [ 6 ] P.L. Ritger and N.A. Peppas, A simple equation for description of solute release II. Fickian and anomalous release from swellable devices, J. Controlled Release, 5 (1987) 37. [7] J.S. Vrentas, C.M. Jarzebski and J.L. Duda, A Deborah number for diffusion in polymer-solvent systems, AIChE J., 21 (1975) 894. [8] C.T. Reinhart and N.A. Peppas, Solute diffusion in swollen membrane. Part II. Influence of crosslinking on diffusive properties, J. Membrane. Sci., 18 (1984) 227. [9] R.A. MacDonald, Modelling macromolecular diffusion through a porous medium, J. Membrane Sci., 68 (1992) 93. [ 10] W.M. Deen, Hindered transport of large molecule in liquidfilled pores, AIChE J., 33 1409. [ 11 ] R.A. Siegel, J. Kost and R. Langer, Mechanistic studies of macromolecular drug release from macroporous polymers I. Experiments and preliminary theory concerning completeness of drug release, J. Controlled Release, 8 (1989) 223. [ 12] M. Giona, A. Adrover and D. Spera, Percolation analysis of hindered diffusion, Chem. Biochem. Eng. Q., 7 (1993) 199.
[ 13] J. Hadgraft, Calculation of drug release rates from controlled release devices, the slab, Int. J. Pharm., 2 (1979) 177. [ 14] N.R. Vyavahare, M.G. Kulkarni and R.A. Mashelkar, Zero order release from swollen hydrogels, J. Membrane Sci., 54 (1990) 221. [15] D. Deutscher, R. Zallen and J. Adler (Eds.), Percolation Structures and Processes, Adam Hilger, Bristol, 1993. [16] D. Stauffer and A. Aharony, Introduction to Percolation Theory, Taylor and Francis, London, 1992. [ 17 ] S. Havlin and D. Ben-Avraham, Diffusion in disordered media, Adv. Phys., 36 (1987) 695. [18] T. Vicksek, Fractal Growth Phenomena, World Science, Singapore, 1989. [19] N. Wax (Ed.), Selected Papers on Noise and Stochastic Processes, Dover, New York, 1954. [20] M. Giona, A. Adrover and A.R. Giona, Finite-difference methods in the analysis of transport phenomena in disordered structures, in Fractals in Natural and Applied Sciences, M.M. Novak (Ed.), North-Holland, 1994. [ 21 ] M. Giona and A. Viola, The pressure~tiffusion equation in the Darcy regime: scaling analysis and anomalies on fractals, in M.M. Novak (Ed.), Fractals in Natural and Applied Sciences, North-Holland, 1994. [ 22 ] M. Giona, A. Adrover and A.R. Giona, Convection~:liffusion transport in disordered structures: numerical analysis based on the exit-time equation, Chem. Eng. Sci., in press. [23] R.S. Harland, C. Dubernet, J.P. Benoit and N.A. Peppas, A model of dissolution-controlled, diffusional drug release from non-swellable polymeric microspheres, J. Controlled Release, 7 (1988) 207. [24] A. Bunde and S. Havlin, Fractals and Disordered Systems, Springer Verlag, Berlin, 1992. [25] A. Adrover, M. Giona and M. Grassi, Controlled release of theophylline from water-swollen scleroglucan matrices, J. Membrane Sci., submitted. [26] P.M. Adler, C.G. Jacquin and J.A. Quiblier, Flow in simulated porous media, Int. J. Multiphase Flow, 16 (1990) 691.
v
ELSEVIER
Journal of Membrane Science 113 (1996) 21-30
Analysis of controlled release in disordered structures" a percolation model Alessandra Adrover a,,, Massimiliano Giona b Mario Grassi c '~Dipartimento di lngegneria Chimica, Universitil di Roma "La Sapien=a", Via Eudossiana 18, 00184 Roma, Itah' h Dipartimento di lngegneria Chimica, Universitd* di Cagliari, Piaz=a d'Armi, 09123 Cagliari, ltal3 • Dipartimento di lngegneria Chimiea, Universitd di Trieste, Pia:=a Europa I, 34127 Trieste, Italy Accepted 9 August 1995
Abstract
The analysis of controlled-release kinetics in disordered models (percolation clusters) of polymer matrices is developed in detail. It is shown that the "anomalous diffusion" experimentally found in the release kinetics from swollen gels, M , / M ~ ~ t", n > 1 / 2 , cannot possibly be related to hindered diffusive motion in disordered media, i.e. is not a consequence of the disordered structure of the matrix. A kinetic model is proposed to account for an exponent n greater than 1/2. This is based on a two-phase kinetics which takes into account both field effects (deriving from potential interaction between solute and polymeric matrix) and entrapping effects due to geometric constraints. Kevwords. Release kinetics: Percolation: Anomalous transport: Two-phase models
1. Introduction In recent years much work has been focused on the development of controlled release devices [ 1 - 4 ] . These drug-delivery systems provide release rates which are determined by the physical properties of the device itself. There exist different kinds of drug delivery systems, amongst which polymeric matrices are widely used. The rate-determining step to produce a specific dosage form is represented by solute diffusion inside the polymeric matrix, depending on the structural properties of the gel. Gel structures can be subdivided into two classes: chemical gels, which form a lattice that is stable also in the presence of a swelling agent; and physical gels, which may undergo structural deformation in time (dissolution). * Corresponding author. 0376-7388/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved S S D 1 0 3 7 6 - 7 3 8 8 ( 95 ) 0 0 2 2 0 - 0
Two different drug-release mechanisms occur in hydrophilic cross-linked polymer systems depending on the macromolecular state of the system, which can be in a glassy and rubbery state during the initial stage of the swelling agent diffusion or in a rubbery state after thermodynamic equilibrium is attained between the network and the surrounding dissolution medium ( usually water) [5]. In the study of swelling controlled systems, the overall solute release depends on the thickness of the swollen layer and on the velocity of the swelling interface. In this situation, the polymeric matrix structure changes in time and solute release develops simultaneously with the propagation of the swelling interface [6,7]. In the case of solute diffusion through swollen gels, the mechanism of transport depends on the relative size of the diffusing species, on the mesh formed by the macromolecular chains [8,9] and on the potential effects related to solute/polymeric
22
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
interactions. Morphological features associated with reduced macromolecular chain mobility and potential barriers to solute diffusion can be responsible of the decrease in the solute diffusivity. Assuming that field effects are negligible, the diffusion coefficient D of the solute molecules in liquid-filled porous matrices deviates from the corresponding coefficient in the dilute bulk solutions DE (Einstein diffusivity) by a factor G, called the accessibility factor, which depends exclusively on steric properties (size exclusion effect, SEE) [10]
D= DEG,
(1)
An approach to this kind of problem can be based on the simulation of transport in disordered structures representing the polymeric matrix: transport parameters are evaluated from reasonable models of the porenetwork [ 11,12]. Besides the hindered effects related to the porosity and to the crosslinking of the matrix, the key problem in the kinetics of controlled release from swollen gels is represented by the anomalous release behaviour from polymeric devices, which generally show, for intermediate time scales, an empiric power-law behaviour [6,13,14],
M,/M~~f', M,/M~ <0.5,
(2)
with an exponent n usually greater than 1/2 (pseudo non-Fickian diffusion). M,/M= is the fraction of solute released up to time t. A value of n > 1/2 is contrary to the classical release models of Fickian diffusion [ 3 ]. In order to understand the origin of the "anomalous" exponent n, we analyze the release kinetics from disordered models of polymeric matrices by taking percolation clusters as the paradigm of a disordered structure [ 15,16]. Percolation lattices are simple and general models of porous structures with spatially uncorrelated disorder [16]. We focus attention on swollen or non-swellable devices in the form of slabs and develop a percolation model of this kind of drugdelivery system in order to evaluate the influence of the complex and disordered nature of the polymeric matrices on the drug-diffusion coefficient and on the exponent n. The cases of lag and burst operating conditions [ 13 ] are analyzed. We show by means of numerical simulations that the transport mechanism in disordered structures is not responsible for the anomalies found in sorption kinet-
ics. In order to explain the apparent power-law release (Eq. (2) with an exponent n greater that 1/2), we propose a simple kinetic model which considers the entrapping effects of solute molecules inside the polymeric matrix, deriving from field interactions or geometric hindrances, in the form of a two-phase kinetics. The numerical results on percolation lattices prove to be in agreement with Eq. (2). The article is organized as follows. First, we discuss the simulation method for analyzing transport phenomena in disordered lattice models. Second, we analyze two typical release conditions (lag and burst effect) in percolation structures in order to evaluate G and n. We then present the two-phase kinetic model, showing that the entrapping kinetics may be responsible for the anomalous behaviour represented by Eq. (2). Finally, we analyze the differences between anomalous transport in fractal and disordered structures (in the sense discussed by Havlin and Ben-Avraham [ 17] ) and the notion of anomalous transport in controlled release systems in order to clarify this problem, which is not only a matter of semantics.
2. Percolation models and numerical simulations A percolation lattice is a spatially uncorrelated random distribution of empty and occupied sites characterized by a percolation probability p representing the frequency of empty sites (site percolation). In the application of percolation models to transport in polymeric matrices, the occupied sites are the elements of the matrix and the empty sites are the pores through which solute particles diffuse. From percolation theory [ 16], it is known that there exists a critical percolation probability, Pc, above which transport is allowed. For values of P>Pc there exists a connected cluster of empty sites (called the infinite cluster) all over the lattice. In two dimensional square lattices Pc = 0.593. For each value of p greater than Pc, in addition to the infinite cluster there exists a family of finite subclusters of various lattice sizes. In order to generate the infinite connected cluster alone, Leath's growth algorithm (spreading percolation) can be adopted [18]. As extensively discussed by many authors [ 15 ], percolation lattices constitute a "model for all seasons" to represent a disordered structure with uncorrelated spatial disorder.
A. Adrover et al. /Journal of Membrane Science 113 (1996) 21 30 The classical approach to the analysis of transport in disordered media and fractals is based on Monte Carlo simulation of random walks, which is grounded on the equivalence between Brownian motion and diffusion [ 19]. The random-walk approach is particularly suitable for analysis of the statistical properties associated with transport (e.g. the scaling of the mean square displacement with time) [ 17] but is of less utility if we want to solve specific transport boundary-value problems, as in the case of sorption kinetic simulations. For this reason, we propose a simple and somewhat elementary finite-difference numerical method for parabolic and elliptic differential equations in disordered lattice models. This satisfies the conservation principle [20] and furnishes promising results in the analysis of linear and nonlinear transport models on fractals and on disordered structures [21,22]. Let us consider the problem of diffusion in two dimensions,
Oc ~-[O~'c 02c~
(3)
where c is the concentration of the transported entity. A generic explicit algorithm for solving Eq. (3) takes the form (n+
%
I ~ --
.(n)
-~i
+fl'
~
jo,)
hk~j,
(4)
h,h } ~ I( i,]
where ¢3' = ~ A t / A x and JJ,'2Lo is the local flux at time in) n, Jhk~"~,j= (~hk'~"~--Gj)/Ax.. Let M be a disordered medium on a square lattice, (id) the lattice coordinates of a generic point and l ( i j ) the set of the first four nearest neighbours of (id). The medium M can be represented by means of its characteristic function
c o( n + l ) =c~j( n )
_~_ [3 ¢
~_,
23
(n} xIKk)J),k~,
(h,k)El(i,j)
'"'+[3
= ¢_!]
E
x(h,k)[d,i"-,,, I
~6)
{ h,k ~ E l ( i , . ] )
where [3 = DA t / A x z. It is important to observe that Eq. (6) satisfies the conservation principle, i.e.
~x(i,j) (i,j}
~
x(h,k)[Gk-c,] =0,
(7}
(It,k) ~ : l ( i , j )
which is a key condition lbr the applicability of a finitedifference scheme deriving from mass conservation. Of course, from the stability analysis of explicit schemes,/3 < 0.5. It should also be noted that the same approach makes it possible to develop more sophisticated numerical techniques (e.g. implicit alternating direction) in a straightforward way. Moreover it is easy to extend the same numerical technique to the case of multiphase transport problems (see section 4).
3. Lag and burst effects: evaluation of the accessibility factor Let us consider the controlled release of a solute from a slab of a swollen polymeric matrix of thickness L, total volume Vo, and void fraction e. The drug is released into a recovery volume V. Drug concentration in the slab is initially equal to zero and in x = 0 the slab is placed in contact with a solution in which drug concentration is kept constant and equal to co (infinite reservoir volume). In this case the polymeric matrix acts as a membrane. In a one-dimensional approach, the balance equation reads as
&,
32c
X( i,j )
X( i,j ) =
l, ( i , j ) ~ M 0, otherwise.
(5)
The initial and boundary conditions describing the lageffect experiment are [ 13 ] : c(x,0) = 0 , 0 < x < L
Eq. (4) is still valid in a disordered structure if care is taken that Eq. (4) is defined only in a point (id') ~ M , (n) i.e. for x ( i J ) = 1, and the local fluxes l~'hk~j are identically equal to zero if x ( i j ) = 0. According to these observations, the simplest explicit algorithm for the analysis of transport in a disordered lattice structure reads
[3c
(9)
hD Oc)
L G=,"
10)
where h = Voe/V. The latter boundary condition takes into account the fact that the recovery volume V is not infinite and consequently stationary conditions will be asymptotically reached leading to a finite value of M-~.
A. Adroveret al. / Journalof MembraneScience 113 (1996)21-30
24
Mt/Moo
From the solution of the differential scheme ( 8 ) - (10), the expression of the fraction of solute M,/M= released up to time t attains the form: y ' B,,exp(
Mr~M== 1
m=l
-
0.75
A2mDt)
~
a
(11) 0.50
~_Bm m=l
/ /
0.25
(A~L2 +h2)c°S(AmL) h,,COS(hmL)=h/L. Bin- L[(A2mL2+h2) + h ] '
a) p-o.95 b) p.o.a5 c) p-0.75 d) p-0.65
0.00 0
(12) For t tending to infinity the leading term in Eq. ( 11 ) corresponds to the smallest eigenvalue hj (fundamental mode) and consequently: (1 -mr~m=) ~ e x p ( - h 2 1 D t ) , t>> 1/(DA2).
(13)
The slope of the curve In (1 -M,/M=) vs. t therefore enables us to estimate the effective diffusion coefficient D. Eq. ( ! 3) is the typical exponential relaxation associated with diffusive motion. Simulation of the lag effect on percolation lattices can be simply performed by following the approach outlined in section 2 and imposing the boundary and initial conditions ( 9 ) - ( 1 0 ) on the lattice. Simulations were performed on a two-dimensional square lattice of size 5 0 × 150 for different values of the percolation probability p ( e = p ) , p > Pc and by choosing h = 0. lp. The size of the lattice is large enough for us to neglect finite-size effects. The choice of a larger size in the direction orthogonal to transport makes it possible to treat the phenomenon as one-dimensional by evaluating the concentration as the average over the direction orthogonal to transport. Fig. 1 shows the resulting sorption curves for different values of the percolation probability p. If s(p,h) is the slope of the curve ( 1 - M , / M = ) vs. t in a lognormal plot (see Fig. 2), then for large t we have
200
400
t
Fig. 1. M,/M~ vs. t for the lag-effect simulations [Eqs. (8)-( 10)] for different values of the percolation probabilityp: L = 1 a.u., D = 1 a.u..
taking into account its internal porosity. Conversely, we can define the accessibility factor GE(p) ( " E " stands for Euclidean in contrast to " p " , percolative) by taking the slab volume Vo as the reference volume and therefore by considering for the factor h the expression h = h E - Vo/V, since the effects of the porosity are encompassed in the accessibility factor. The relationship between G" and G E is therefore given by GE(p) = GP(p)p.
(15)
Fig. 3 shows the comparison of the values of GE(p) obtained by applying Eq. (14) to sorption kinetic simulations with the values of GE(p) obtained from Monte Carlo simulations [ 12]. The figure shows a fairly good level of agreement. It should also be observed that while s(p,h) and Al are functions o f p and h, the accessibility factor is only a function of p, i.e. it depends exclusively on the steric properties of percolation lattices and does not depend on the boundary condition. The burst con1-( Mr/M** ) 10 o ~
~
al °-O.g5 b) p-o.a5
E\\\
c,o.o, d,°-o°,
s(p,h) GP(p) =
AE(h,L ),
(14)
where GP(p) is the value of the accessibility factor for a fixed value p of the percolation probability evaluated by taking as reference volume the effective volume of the slab Vo'p in which the solute diffuses. The superscript " p " indicates simply that we assume the effective volume of the slab as the reference volume by
1o0[
\ 100\
200 t 300 Fig. 2. Log-normal plot of 1- M,IM=vs. t for the data of Fig. 1.
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30 GE
2oMC
10
~
burst effect
•
,oo,-,c,
05 i
i i
0.0
.
.
.
.
.
s(p) rfl/4L2.
GP(p) -
0 ql
.
.
.
00
.
.
.
.
.
.
.
0 5 (p-pc)/(1-Pc) 10
Fig 3 Comparison of the values of GE(p) obtained from Monte Carlo simulation [ 12] with the values of GE(p) for lag and bursteffect s i m u l a t i o n .
dition represents another physical situation often encountered in controlled release systems: the matrix is initially loaded up to a given solute concentration Co. This observation is confirmed by simulation results of controlled release in the burst condition [13], The transport equation describing this phenomenon is still given by Eq. (8), to be solved with the following initial and boundary conditions:
c(x,O)=co, O
(20)
In the case of burst-effect conditions, the accessibility factor is evaluated with respect to the effective slab volume allowed to the solute diffusion, i.e. with the above notations it is given by GP(p). The corresponding values of GE(p) are shown in Fig. 3, in agreement with the lag-effect simulation experiments. The second important question that can be addressed by means of percolative models is the evaluation of the exponent n appearing in Eq. (2) from sorption experiments on disordered lattice structures. By considering, for example, the burst effect (a similar result can be obtained from lag conditions), simulations on percolation lattices give a value of n equal to 0.5 for p at a distance from the critical threshold. This is evident from the slope of the curves ofM,/M~ ( forM,/Ms < 0.5 ) plotted 1-( Mr/Moo ) 10° %
\'\
-
\\
~\ '\
(16) (17)
The latter boundary condition comes from the assumption of perfect sink conditions in x = L. The solute fraction released up to time t attains the form
a
~
'
"\ e
0~xC, _ o = 0 , c(L,t) =0.
25
a) p-100 b) p-0.90 d
'c
c) p-080
e) p-065
d) 1:}-0.70
10-2 0 1 2 3 4 t 5 Fig 4 Log-normal plot of I -M,/M~ vs t for the bursteffect sim ulations for different values of the percolation probability p: L = I a u,D=l a u oMt/M=
M,/M=
1
8 77"
_
=
1 I
(2m
-
1 )2
exp
( 2m--1)27r2Dt] 4L=
]
(18)
'° I 1O- 1
and
In( 1 - M,/M~) ~ - ~ _ D t .
a) p-l.00 b) p-0.90 c) p-0,80 d) p-070 e) p-0.65
(19)
Again, if s(p) is the slope of the curve ( 1 - M , / M = ) vs t in a log-normal plot (see Fig. 4), for large values of t we have
/~/~//~
102
........ ................................. 10 -4 10 -a 10 "a 10 -~ 100 t 101 Fig 5 Log-log plot of M,/M= vs t for a burst-effect simulation for different values of the percolation probability p
26
A. Adroveret al. / Journal of MembraneScience 113 (1996) 21-30
in log-log scale, Fig. 5. Only forp = 0.65 does the slope change slightly and attain the apparent value of n = 0.38 _+0.03. This result can be explained by considering that for p tending toward criticality (Pc = 0.593) the highly heterogeneous nature of the porous matrix emerges, as appears from anomalous characteristics in diffusive motion (n4: 1/2). However, the resulting exponent n is always less than 1/2. An analysis of this topic is developed in section 5. To sum up, the apparent deviation from Fickian behaviour experimentally observed in many controlled-release devices cannot be explained by means of anomalous transport effects in disordered matrices. Consequently, the physical origins of the anomalies in M,/M=-kinetics should be sought in other controlling mechanisms.
4. Two-phase kinetic model The above results indicate that an exponent n greater than 1/2 is not related to anomalous diffusion effects in disordered media. Deviation from n = 1/2 can be due to different phenomena: concentration dependentdiffusion coefficient, erosion of the gel matrix, convective contribution to mass transport, etc. In the case of non-erodible matrices, under the assumption of a practical concentration-independent diffusion coefficient and negligible convective effects, a simple way of explaining the phenomenological Eq. (2) is to regard it as a consequence of interaction effects between solute and polymeric matrix. This can be modeled by assuming the presence, in the release matrix, of two phases: a gel phase ("entrapping phase") and a sol phase in the liquid solution. Two-phase models in the analysis of controlled release have been considered by other authors [ 13,23 ] in order to describe the effects of interfacial transfer kinetics on release rates for those systems which show significant free energy barriers for transport across liquid-liquid interfaces [ 13 ] or to model the release-condition above the solubility threshold [ 23 ]. In the proposed model we assume that a fraction es of solute molecules is initially entrapped in the gel phase and irreversibly released in the sol phase in which it diffuses with an effective diffusion coefficient D. The phenomenon is close to an adsorption and irreversible desorption first-order kinetics with a kinetic coefficient
k. The dynamic equations describing the phenomenon are: Oc 02c Oq o t = Dox 2 0 t
(21)
Oq --= -kq, Ot
(22)
where c(x,t) and q(x,t) are drug concentrations in the gel and sol phases respectively. Eqs. (21) and (22) are to be solved with the initial and boundary conditions (burst effect): c(x,O) -----(1--E~)coO
(23)
0c
(24)
Xx=
=0, c(L,t) =0. 0
In the simulations on percolation lattices, the stationary phase is represented by the first nearest neighbouring site of each site representing the polymeric matrix (see Fig. 6) and is therefore implicitly a function of position since it is defined at the matrix boundary sites. Fig. 7 shows the simulation results of the two-phase kinetics on the percolation cluster for different values of the kinetic constant k ( in arbitrary units), keeping constant the values of the initial partition coefficient es and of the percolation probability p (es = 1, p = 0.80). Fig. 8 shows the behaviour of the exponent n as a function of k for different values of p. For high values of k the internal release kinetics is so fast that the presence of
I
polymericmatrix~
sites
matrixboundary
sites
Fig. 6. Schematicrepresentationof two-phasesimulation in percolation lattices. In the simulationsthe concentrationq is definedonly on the nearest neighbors of a matrix site (matrix boundary site) while c is defined in all the sites not belonging to the polymeric matrix.
A. Adrover et al. /Journal of Membrane Science 113 (1996) 21-30
27 ]
100 Mt/M.__~
n /
1.2 !m ,i /
0.8-
10-; f bc
[°
~' ,//
.I k.5o0 b) k-2oo
e ~
c) k-lO0
d) k-50
e) k-lO
f) k-1 ,
¢
i m_~ 0.4 i
!
102L .... , ,,, , , , ,, , 10-2 10-1 10 o t 101 Fig. 7. Log-log plot of M,/M= vs. t for the two-phase percolation model for different values of k (a.u.) and fixed values o f p and e,
i
0.0
(p = 0.80, e+= 1).
~
0.00
0.25
.
.
.
0.50
.
0.75
1.00
6~
Fig. 10. Plot of n vs. ~+ for a fixed value of p and k (p = 0.80, k = 1 1.2
a.u.). a
\
n
1.0
c
0.8~ 0.6 0.4 0.2
a) p-0.90 b) p-0.80 c) p-0.70 d) p-0.65 10 -1
10 o
101
102
10 a k
104
Fig. 8. Plot o f n vs. k (a.u.) for a fixed value o f p and ~+ (p=0.80, am=l).
100 F
Mt/M°°
i
10 -1 ~s.O. 1 i i
10 -2 !
//
/
r-
/
[ 10-3~
/
/ /'
rC+.lO ....
i
F
'
10 -+ . . . . . . . . . . . . . . . . . . . . . . . 10-3 10"2 10-1 10° t
,ll
101
Fig. 9. Log-log plot of M,/M= vs. t for the two-phase percolation model for different values of e, and fixed values ofp and k (p = 0.80, k = l a.u.).
the gel phase is practically negligible and n tends to the Fickian value 0.5. The influence of the initial partition coefficient es of the solute molecules between the two phases is evident from Fig. 9, which shows the simulation results of the two-phase model for different values of a,, keeping constant the percolation probability and the kinetic constant (p = 0.80, k = 1 a.u. ). Evaluation of the slopes of the sorption curves of Fig. 9 for intermediate time scales reveals that n is a monotonically increasing function of es (see Fig. 10). Therefore, we may conclude that a simple kinetics model which takes into account the internal release kinetics between a sol and a gel phase is able to predict a release exponent n greater than 1/2. As the simulations show, the value of n depends on the kinetic parameters introduced in the model. It should be noted in passing that the predictions of the model proposed have been checked against experimental results on controlled release kinetics by developing a detailed data analysis. In [25] it is shown that the predictions of the model (21 ) - ( 2 2 ) are in agreement with experimental results on the release rate of theophylline from water-swollen scleroglucan matrices.
5. " A n o m a l o u s transport" in fractal m e d i a and in controlled release s y s t e m s
The empirical relation (2) is usually regarded as the deviation from Fickian behaviour in controlled release
28
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
systems. Since deviations from Fick's law are also observed in diffusion through disordered fractal structures, and since polymeric matrices can be regarded as disordered media, it is useful to analyze the differences between these two "anomalous transport phenomena" critically in the light of the results obtained in the previous sections. In the theory of dynamic phenomena on fractals, anomalous transport means deviation from Fick's law in the specific sense that the mean square displacement of a Brownian particle, (r2(t)), deviates from a linear behaviour with time and follows a scaling law,
(r2(t)) - t ~,
(25)
with an exponent/3 less than unity. Eq. (25) is valid asymptotically, i.e. for t ~ w. For some classes of deterministic fractals it is possible to obtain explicit expressions for the exponent/3 by means of renormalization group techniques (decimation methods) [ 17]. The propagation of a diffusion front in a fractal medium does not display Gaussian behaviour but has a stretched Gaussian form. For more details on all these topics see [ 17,24]. In Eq. (25) /3 is a universal exponent characteristic of the fractal structure of the medium. Eq. (25) implies that the heterogeneity at all the length scales of a fractal structure induces a characteristic slowing down (/3 < 1 ) in the diffusional propagation for all the time scales. This effect can be interpreted as the appearance of correlations (for uncorrelated processes/3 is strictly equal to unity) induced exclusively by geometric and topological disorder, although particle dynamics at microscopic scales is still Brownian. Anomalous transport in swollen polymers (the case of simultaneously swelling front propagation is not considered here since the phenomenology is slightly different) also shows an exponent n which deviates from the classical value n = 1/2 (which is strictly the counterpart of/3 = 1 for the scaling of the mean square displacement), but the physical origin of this deviation is completely different from the case of diffusion in disordered or fractal media. The scaling law (2) is not valid asymptotically, as in the case of Eq. (2), and is not a manifestation of universal physical properties, as in the case of the exponent/3. The simulations on percolation lattices developed in the previous sections have shown that for these disordered structures near criticality the exponent n in pure diffusion is definitely less than 1/2, as should be
expected since disorder tends to slow down diffusive motion. Consequently, an "anomalous" exponent n greater than 1/ 2 (which implies an enhancement of the release) cannot possibly be related to non-Fickian diffusion as in the case of fractal structures. The phenomenological origin of this exponent should be sought in the presence of correlations in the diffusive motion deriving from entrapping or field effects (which ultimately have the same physical origin). For these reasons, we have introduced in section 5 the simplest physical model consistent with these assumptions, and the simulation results confirm that such a model can predict an apparent scaling exponent n greater than 1 / 2. The simulation results confirm that the exponent n is not universal but depends on the kinetic constant k. To sum up, the appearance of a pseudo "non-Fickian" exponent n > 1/2 can be regarded as the superposition of regular Fickian transport and of kinetic effects described macroscopically with simple twophase model, but whose physical origin depends ultimately on the specific interactions between solute molecules and polymeric matrices.
6. Conclusions and discussion The expression "anomalous transport" usually refers to two different transport situations: the pseudo non-Fickian behaviour in controlled release kinetics (appearing as an exponent n > 1 / 2 in the sorption curve scaling with time) and the theory of transport in disordered and fractal structures, which are characterized by an exponent/3 less than unity in the scaling of the mean square displacement with time. It is intuitive to suppose that the pseudo non-Fickian behaviour observed in many controlled release experiments could be attributed to the disordered structure of the polymeric matrix. The simulations of the sorption kinetics in disordered models (represented by percolation clusters) in the case of both lag and burst conditions have clearly shown that an exponent n greater than 1/2 cannot be related to the disordered geometric and topological structure of the matrix. Indeed, in the case of percolation clusters near criticality, we observed an exponent n less than 1/2, as is to be expected from the scaling theory of transport in fractal structures [ 17,24].
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
This result is fairly important in the case of drug diffusion through swollen polymeric matrices (since in this case the anomalies cannot be attributed to the simultaneous propagation of the gel front) because it indicates that the origin of the apparent non-Fickian behaviour should be attributed to other kinetic mechanisms rather than to the microscopic disordered structure. We therefore propose a very simple kinetic model by including the effects associated with the solute/ polymeric matrix interactions. In macroscopic models, the interactions can be described by means of a twophase kinetic model, as developed in section 4. It is important to observe that the entrapping effects expressed macroscopically by means of two-phase models may be related to different physical mechanisms (see e.g. the different approaches in the analysis of usual adsorption phenomena): the presence of potential energy fields related to solute/matrix interactions, the entrapping effects related to hindered constraints and the physical adsorption of the solute molecules on the polymeric matrix. The analysis of these important aspects involves detailed microscopic modeling of the specific solute/matrix interactions, which goes beyond the scope of investigations of macroscopic models. However, in a macroscopic perspective it is important to observe that the simple kinetic model we propose is able to account for the fundamental features of the apparent non-Fickian release. Moreover, experimental results on the release of theophylline from water swollen scleroglucan matrices [25] confirm the results predicted by the two-phase model presented here. A final comment is related to the simulation method adopted in dealing with transport in disordered lattice models. The numerical method discussed, which has furnished good results in the simulation of Brownian motion on fractals [ 20], represents an improvement on Monte Carlo simulation and can be applied to a great variety of transport phenomena of interest in controlled-release modelling (e.g. the analysis of potential fields on release kinetics from polymeric materials). Of course, if the structural parameters of the polymeric matrix are known, this technique can be applied to a microscopic lattice model of the polymeric structure by constructing a lattice simulator of the matrix through the reconstruction method of porous solids proposed by Adler et al. [26].
29
7. List of symbols Bm e Co
D G g~ Gp h J k L M M,/M~ n
P p¢. q
(r~(t)) s(p,h ) t V
Vo X
see Eq. ( 11 ) drug concentration initial drug concentration effective diffusion coefficient see Eq. 3 Euclidean accessibility factor percolative accessibility factor
Voe/V local flux kinetic coefficient gel slab thickness disordered medium fractional solute release anomalous exponent [ see Eq. ( 2 ) ] percolation probability percolation threshold drug concentration in the gel phase mean square displacement see Eq. (14) time recovery volume gel slab total volume spatial coordinate
7.1. Greek letters
[3
universal exponent characteristics of a fractal medium DAt/Ax
A1
X( ij)
void fraction initial drug partition coefficient between the sol and the gel phase fundamental mode [see Eq. (13) ] characteristic function of the porous medium M
Acknowledgements We should like to thank A.R. Giona for constant encouragement and useful discussions. One of the authors (M.G.) would like to thank M. Liauw for useful discussions.
30
A. Adrover et al. / Journal of Membrane Science 113 (1996) 21-30
References [ 1] R.W. Baker and H.K. Lonsdale, Controlled release: mechanism and rates, in A.C. Tanquary and R.E. Lacey (Eds.), Controlled Release of Biologically Active Agents, Plenum Press, New York, NY, 1974, p 73. [2] R.S. Langer and N.A. Peppas, Chemical and physical structure of polymers as carriers for controlled release of bioactive agents: A review, J. Macrom. Sci., Rev. Macromol. Chem. Phys., C23 (1983) 61. [ 3 ] J. Crank, Mathematics of Diffusion, 2nd ed., Clarendon Press, Oxford, 1975, [4] L.S. Langer and D.L. Wise (Eds.), Medical Application of Controlled Release: Application and Evaluations, Vol. I-II, CRC Press, Boca Raton, FL, 1984. [5] P.I. Lee, Kinetics of drug release from hydrogel matrices, J. Controlled Release, 2 (1985) 277. [ 6 ] P.L. Ritger and N.A. Peppas, A simple equation for description of solute release II. Fickian and anomalous release from swellable devices, J. Controlled Release, 5 (1987) 37. [7] J.S. Vrentas, C.M. Jarzebski and J.L. Duda, A Deborah number for diffusion in polymer-solvent systems, AIChE J., 21 (1975) 894. [8] C.T. Reinhart and N.A. Peppas, Solute diffusion in swollen membrane. Part II. Influence of crosslinking on diffusive properties, J. Membrane. Sci., 18 (1984) 227. [9] R.A. MacDonald, Modelling macromolecular diffusion through a porous medium, J. Membrane Sci., 68 (1992) 93. [ 10] W.M. Deen, Hindered transport of large molecule in liquidfilled pores, AIChE J., 33 1409. [ 11 ] R.A. Siegel, J. Kost and R. Langer, Mechanistic studies of macromolecular drug release from macroporous polymers I. Experiments and preliminary theory concerning completeness of drug release, J. Controlled Release, 8 (1989) 223. [ 12] M. Giona, A. Adrover and D. Spera, Percolation analysis of hindered diffusion, Chem. Biochem. Eng. Q., 7 (1993) 199.
[ 13] J. Hadgraft, Calculation of drug release rates from controlled release devices, the slab, Int. J. Pharm., 2 (1979) 177. [ 14] N.R. Vyavahare, M.G. Kulkarni and R.A. Mashelkar, Zero order release from swollen hydrogels, J. Membrane Sci., 54 (1990) 221. [15] D. Deutscher, R. Zallen and J. Adler (Eds.), Percolation Structures and Processes, Adam Hilger, Bristol, 1993. [16] D. Stauffer and A. Aharony, Introduction to Percolation Theory, Taylor and Francis, London, 1992. [ 17 ] S. Havlin and D. Ben-Avraham, Diffusion in disordered media, Adv. Phys., 36 (1987) 695. [18] T. Vicksek, Fractal Growth Phenomena, World Science, Singapore, 1989. [19] N. Wax (Ed.), Selected Papers on Noise and Stochastic Processes, Dover, New York, 1954. [20] M. Giona, A. Adrover and A.R. Giona, Finite-difference methods in the analysis of transport phenomena in disordered structures, in Fractals in Natural and Applied Sciences, M.M. Novak (Ed.), North-Holland, 1994. [ 21 ] M. Giona and A. Viola, The pressure~tiffusion equation in the Darcy regime: scaling analysis and anomalies on fractals, in M.M. Novak (Ed.), Fractals in Natural and Applied Sciences, North-Holland, 1994. [ 22 ] M. Giona, A. Adrover and A.R. Giona, Convection~:liffusion transport in disordered structures: numerical analysis based on the exit-time equation, Chem. Eng. Sci., in press. [23] R.S. Harland, C. Dubernet, J.P. Benoit and N.A. Peppas, A model of dissolution-controlled, diffusional drug release from non-swellable polymeric microspheres, J. Controlled Release, 7 (1988) 207. [24] A. Bunde and S. Havlin, Fractals and Disordered Systems, Springer Verlag, Berlin, 1992. [25] A. Adrover, M. Giona and M. Grassi, Controlled release of theophylline from water-swollen scleroglucan matrices, J. Membrane Sci., submitted. [26] P.M. Adler, C.G. Jacquin and J.A. Quiblier, Flow in simulated porous media, Int. J. Multiphase Flow, 16 (1990) 691.