- Email: [email protected]
Development and temporal evolution of erosion fronts in bioerodible controlled release devices
0009-2509/90 9;3 .00 + 0.00
Chemieul Engineering Science, Vol . 45, No . 8, pp . 2 359 -2366, 1990 .
T'. 1990 Pergamon Press plc
Printed in Great Britain .
DEVELOPMENT AND TEMPORAL EVOLUTION OF EROSION FRONTS IN BIOERODIBLE CONTROLLED RELEASE DEVICES Kyriacos Zygourakis Department of Chemical Engineering Rice University Houston, Texas 77251-1892
ABSTRACT Cellular automata and discrete iterations are used to model and simulate the release of bioactive agents (drugs) from bioerodible pellets used as controlled release devices . The presented discrete simulation approach is a computationally efficient tool for designing delivery systems with optimal release properties, since it can accurately treat the transient erosion processes in multicomponent release systems of arbitrary geometry and with different dissolution rates for each component . Results from two-dimensional simulations are presented to show how the intrinsic dissolution rates, drug loadings, matrix porosities etc . affect the overall release rates and to establish the relationship between the ratio of dissolution rates and the fractal dimension of the erosion front. KEYWORDS Cellular Automata, Dissolution, Controlled Release of Drugs, Erosion Fronts, Fractal Dimension . INTRODUCTION Drug treatment of medical problems is most often accomplished by administering a drug at fixed time intervals, a method that causes wide oscillations of the drug concentration in the bloodstream . After a dosage is administered, drug levels rise, reach a peak and decrease continuously thereafter . The drug loses its therapeutic effect when its concentration falls to very low levels . If we attempt to prevent this by increasing the drug dosage, however, the peak concentration may exceed the toxic level causing undesirable side effects . An ideal controlled or sustained release device should deliver the drug at constant rate, thus maintaining the drug level in the blood within a narrow therapeutic range . Several studies (see reviews by Langer, 1981 ; Rosen et al ., 1988) have developed polymeric and other systems that can be used to fabricate controlled release devices for specific drugs . Currently, however, the design of devices with specific release characteristics requires extensive experimentation and the associated large development costs lead to high prices for the final pharmaceutical products . This communication will present a modeling approach that can systematize and facilitate the design of controlled release devices with desired properties . We will consider chemically controlled release systems where particles of a bioactive agent (drug) are dispersed in a porous or nonporous bioerodible matrix consisting of a polymer or some other inert filler material . When the release matrix is exposed to an environmental fluid (solvent), the filler erodes exposing previously inaccessible drug particles . These particles will then dissolve releasing the bioactive agent into the solvent phase (Langer, 1981) . Let mdo and md(t) be respectively the total amount of drug initially present in the device and the amount released from the device at time t . An ideal controlled release device must be designed so that the fractional release X defined by x = ( md(t) / mdo) is a linear function of time with slope equal to some desired value. When the drug and matrix material have low solubilities, surface detachment is the rate determining step for the erosion processes occurring at the solvent/drug and solvent/filler interfaces (Cooney, 1972) . If Vd and V p are the volumes of dissolved drug and polymer (filler) respectively, the two dissolution rates can be expressed as Rd
dVd = ka Sd (cd dt Pd
Rp _ dVp _ tit
kp Sp (
cd) = kd Sd ACd = vd Sd Pa
4 - 4)
pp
_ kp Sp ilcp _
pp
SP VP
(1)
(2)
where kd, kp are the two dissolution rate constants, Sd, S p are the drug/solvent and filler/solvent interfacial areas, Ach, oc p are the differences between the saturation and bulk concentrations of the drug and filler and Pd, pp are the drug and filler densities respectively . Finally, vd = (kd Acj ) / pd and v p = (kPAcp ) / p p and the ratio of dissolution rate constants is defined as S = ( vd / V p) . If we assume that there are no diffusional resistances to the transport of the dissolved drug to the bulk of the fluid, the drug release rate equals the drug dissolution rate Rd given by Eq . (1) and the specific drug release rate R* d is given by
2359
KYRIACOS ZYGOURAKIS
2360
dvd = Rd = _L Vo dt V0
Cl
Vd Sd
(3)
where Vo is the initial total solid volume (drug plus filler(s)) in the release device . Previous studies (Cooney, 1972 ; Langer, 1981) considered a simplified dissolution problem . If Vp (2) may be combined to obtain
=
vd=-v, Eqs. (1) and
dVt = v St tit
(4)
where Vt = Vd+Vp and S t = Sd+S p. If Sd = aS t where a is a constant, the release rates can be easily obtained by following the temporal evolution of the overall interfacial area S t for the drug/polymer matrix . This simple problem has been solved for release matrices with slab, cylindrical or spherical geometries (Cooney, 1972) . In practice, however, Vd z vi, and Eqs . (1) and (2) must be solved simultaneously, a problem that is significantly complicated when interfacial areas S d and S p vary with time . When 8 ;& 1 (i .e . the drug and the polymer matrix dissolve at different rates), erosion fronts with tortuous and continuously changing morphology will develop (Zygourakis, 1988) . The temporal evolution of the interfacial areas and, thus, the release rates will be strongly affected by several chemical and structural properties of the release system such as dissolution rates, the matrix porosity, the connectivity of the pore network and the dispersion of the drug . ('FT I ULAR AUTOMATA MODELS We will solve the dissolution problem described by Eqs . (1) and (2) using models based on cellular automata and discrete iterations . The evolving structure of multicomponent release systems is simulated on rectangular arrays of square computational cells . Each cell can exist in one of four possible states : it can represent a small volume element of drug, polymer, solvent or pore void . By assigning the same state to groups of adjacent cells, domains of arbitrary shape or size can be created on the initial cellular arrays to model drug particles, pores, the polymer matrix or solvent-filled regions . Fig . 1 is a pictorial representation of a 1024x1024 cellular array modeling a cross-section taken through a cubical "tablet" prepared by compressing a mixture of drug particles (black clusters in Fig . 1) and "filler" material particles (light grey areas) . The pores created by the pelletization process are indicated by the white clusters in the interior of the tablet . All sides of the tablet except the top are coated with an insoluble film (dark grey strip in Fig . 1 ; see also sides of Figs . 2,3 and 4). Finally, the top side of the tablet is exposed to a solvent (white phase) .
Fig . 1 :
Image showing the initial configuration of a cellular array modelling the cross-section of a controlled release tablet.
The initial configuration of the cellular arrays may be chosen to match the structural properties of a given controlled release device . For a tablet consisting of a drug and a filler with volume fractions f d and fp respectively, the total areas of drug Ad and polymer A that will be measured on an arbitrary planar cross-section can be estimated from the following formulae (Weibel, 1980 fd = A ~fp =
At
AP At
(5)
where A t is the total area of the cross section . If we also use the easily measured values of tablet porosity and particle (domain) size for each phase, we can construct cellular arrays having any desired values of e, A d and A p by randomly
Development and temporal evolution of erosion fronts
Cl
2361
distributing on a square lattice circles with a known size distribution (Weibel, 1980) . We followed this procedure to construct the array of Fig . 1 that models a cube with side equal to 10 mm where 150 - 212 lun drug particles are imbedded in a porous (e = 061) matrix. Image analysis techniques were used to check the grid properties . We used 1024x1024 and 2048x2048 cellular arrays for the two-dimensional simulations presented in this study and each cell corresponded to a square with sides equal to 10 or 5 µm respectively . The high spatial resolution of our simulation approach allows us to accurately describe the transient dissolution phenomena . The initial cellular array is updated at equally spaced time instants tr, t 2,.. ., tr, tr+t, ... as follows : If at time tr a drug or polymer cell is adjacent to one or more solvent cells, then this drug or polymer cell may be changed to a solvent cell (or dissolved) at the next time level tr+ 1 = tr + at. The dissolution of drug or polymer cells proceeds at different rates according to either a probabilistic or a periodic deterministic algorithm that changes the state of different cells at unequal time intervals that are proportional to the corresponding dissolution rates . Similarly, pore void cells are changed to solvent cells to simulate the penetration of pore space by the solvent . Figs. 2, 3 and 4 present the configurations of the cellular array of Fig . 1 at various levels of fractional release . Note that both the drug loading (40% wt .) and the ratio of dissolution rate constants (8 = 10) are high for this simulation . Due to the relatively large drug clusters (black areas) and the large difference between the dissolution rates of the drug and the filler, a highly tortuous erosion front develops leading to significant enhancements of the dissolution rates as we will see later.
Fig. 2:
Image showing the configuration of the cellular array of Fig . 1 after 25% of the drug has been released ( x = 0.25 ) . Note the highly tortuous erosion front .
Fig. 3:
Image showing the configuration of the cellular array of Fig . 1 after 50% of the drug has been released (x = 0.50 ).
2362
C1
KYRIACOS ZYGOURAKIS
Fig. 4:
Image showing the configuration of the cellular array of Fig . 1 after 75% of the drug has been released (x = 0 .75 )-
If the dissolution rate is uniform over the entire drug/solvent interfacial area, a layer of drug with uniform thickness dz will be dissolved over a time period dt . Thus, dVd = Sd dz and Eq . (1) becomes -dZ kd Acd = = Vd dt Pd
(6)
If every simulation step dissolves a single layer of cells at the drug/solvent interface, Eq . (6) implies that every step corresponds to the time interval At necessary to dissolve a solid drug layer with thickness equal to the cell size Az, or AZ ^ Vd
or
At
At = Az Vd
(7)
The drug release rate can be easily approximated after each grid update as (tr) = 1 (AN )` Rd
No
(8) At
where (AN)r is the number of drug cells dissolved in the r-th grid update . A dimensionless release rate rd` and the fractional release x can finally be computed as (AN)' (AN)r rd (tr) _- R d (tr )At = No
x(tr ) =
k-t
(9) Ndo
where Ndo is the total number of drug cells contained in the initial cellular array . RESULTS AND DISCUSSION Fig . 5 presents the fractional release vs time for a series of porous release devices (e = 0 .061) with different drug loadings and different ratios a of dissolution rate constants . Since the solvent will quickly penetrate the pore voids, the connectivity of the pore network is important . For these runs, we assumed that the internal pore voids are not initially accessible to the solvent . They open up, however, as the erosion front moves through the slab . Almost zero-order release (straight lines of x vs. time) is predicted for all values of 8 . The slight deviations from linearity observed on the release curves are attributed to the presence of large pores . As these pores are rapidly penetrated by the solvent creating characteristic "fingers"and increasing the drug/solvent and filler/solvent interfacial areas, the release rates may jump to higher values . The drug release rates decrease significantly with increasing S due to a "shielding" effect . Low filler dissolution rates slow down the movement of the erosion front and decrease the rate at which drug particles are exposed to solvent. This
C1
2363
Development and temporal evolution of erosion fronts
shielding effect diminishes, however, as the drug loading increases . High drug loadings lead to the formation of large drug clusters that are rapidly dissolved when exposed to the solvent . We would expect that better dispersion of the drug (in the form of smaller particles) should partially compensate for this effect .
C9 N d O
m
W C O W O i CL
1500
1000
500
Dimensionless Time
I
Fig . 5 :
r
.
r
,
I
,
r
,
r
500 1000 Dimensionless Time
1500
500 1000 Dimensionless Time
1500
Porous Matrices : The effects of drug loading and ratio 8 of dissolution rate constants on the predicted fractional release curves .
Fig. 6 presents the fractional release curves for non-porous matrices with the same drug loadings and ratios of dissolution rate constants as those considered in Fig . 5 . The absence of internal pores leads to lower release rates especially for the higher values of 8 . Also, the release curves are now smoother as we would expect from the previous discussion .
Cl
KYRIACOS ZYGOURAKIS
2364
10
(A) 20% drug loading 500 1000 Dimensionless Time
I
I
1500
.
500 1000 Dimensionless Time
0
i
0
Fig . 6:
I
500 1000 Dimensionless Time
1500
,
1500
Non-porous Matrices : The effects of drug loading and ratio S of dissolution rate constants on the predicted fractional release curves .
The release curves presented in Figs . 5 and 6 were fitted with a straight line (for x between 0 .05 and 0 .95) to compute the average release rates
2365
Development and temporal evolution of erosion fronts
Cl
10000 a> r ca OC: a)
1000
N
01
01 Lr N 0, t5 L d
100
Q
10 2
4
6
10
8
12
Ratio of Dissolution Rate Constants The effects of drug loading, matrix porosity and ratio 8 of dissolution rate constants on the average release rates (Closed symbols : porous samples ; open symbols: nonporous samples) .
Fig . 7
Fig. 8 presents instantaneous dissolution (release) rates computed according to Eq . (9) for some typical simulations . Note that the dissolution rates oscillate widely (even when 8 = 1 ) as the moving erosion front progressively exposes drug particles to the solvent .
0.007 AD
0 0 .006 N
0 .005
to
m 0 .004 Qf
¢ 0.003 v
w 0.002 C
Q 0 .001 to 0
200 Fig . 8 :
400 600 Dimensionless Time
800
1000
Temporal evolution of dissolution rates for some typical simulations with a porous matrix .
In order to systematize the analysis of simulation data and to answer some of the fundamental questions pertaining to the temporal evolution of dissolution rates, we studied the fractal dimension of the erosion fronts . The fractal dimension provides a quantitative measure of the roughness of the erosion front and should correlate well with the total dissolution rate (i.e. the sum of drug and filler dissolution rates here) . We identified the erosion front (i .e the union of drug/solvent and filler/solvent interfaces) using dilation and subtraction operations on binary images of the cellular array and computed its fractal dimension at several release levels for each run using Flook's method (Flook, 1978) . This method is based on dilation operations and is highly suitable for image processing systems . Fig . 9 shows how the fractal dimension of the erosion front evolves with the release level for three different values of 8 and a sample with 40% drug loading. The average values of the fractal dimension are also shown as straight lines on the same graph . Although the morphology of the erosion front predicted for a given cellular array and set of parameters changes drastically with time, its roughness (as measured by the front's fractal dimension) remains at fairly constant levels . But, the fractal dimension depends strongly on the problem parameters (drug loading, particle size and ratio of dissolution constants 8) as shown in Fig . 10 . ACKNOWLEDGMENT This work was supported by a grant provided by the Texas Advanced Technology Program .
CES 4$f$-ae
2366
KYRIACOS ZYGOURAKIS
Sample 2 :
CI
40% drug loading
1 .4 C 0 1 .3 _A 0 C ar E 0 1 .2
A A
S
8= 10
a
A 5
0
n m
∎
1 .1
∎
8 = 2
∎
∎ ∎
U-
.∎∎
∎ ∎
1 0
Fig . 9 :
0.2
0 .4 0.6 Fractional Release
0 .8
1
Temporal evolution of the fractal dimension of the erosion fronts obtained for different values of the ratio S of dissolution rate constants .
2
4 6 8 10 12 Ratio of Dissolution Rate Constants
14
Fig . 10 : The effects of drug loading and ratio S of dissolution rate constants on the fractal dimension of the erosion fronts .
REFERENCES Cooney, D.O. (1972). Effect of Geometry on the Dissolution of Pharmaceutical Tablets and Other Solids : Surface Detachment Kinetics Controlling_ AIChE J., Za, 446-449 . Flook, A .G. (1978) . The Use of Dilation Logic on the Quantimet to Achieve Fractal Dimension Characterisation of Textured and Structured Profiles . Powder Technol ., 21, 295-301 . Langer, R . (1980) .Polymeric Delivery Systems for Controlled Drug Release . Chem. Eng. Commun ., ¢, 1-48 . Rosen, H .B ., J. Kohn, K . Leong and R . Langer (1988) . Bioerodible Polymers for Controlled Release Systems . In : Controlled Release Systems : Fabrication Technology (D.S .T. Hsieh, Ed.). vol . II, pp . 83-110, CRC Press, Boca Raton, FL. Weibel, E .R . (1980) . Stereological Methods . Vol . 2. Academic Press, New York . Zygourakis, K . (1988) . On Modeling the Controlled Release of Drugs from Bicerodible Polymer S ystems. Proceed.
Intern . Symp . Control. Rel. Bioact . Mater., 11, 306-307 .