- Email: [email protected]
A general code to predict the drug release kinetics from different shaped matrices
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ejps
A general code to predict the drug release kinetics from different shaped matrices Anna Angela Barba a , Matteo d’Amore a , Serafina Chirico b , Gaetano Lamberti b,∗ , Giuseppe Titomanlio b a b
Dipartimento di Scienze Farmaceutiche, Università degli Studi di Salerno, Via Ponte don Melillo, I-84084 Fisciano (SA), Italy Dipartimento di Ingegneria Chimica e Alimentare, Università degli Studi di Salerno, Via Ponte don Melillo, I-84084 Fisciano (SA), Italy
a r t i c l e
i n f o
a b s t r a c t
Article history:
This work deals with the modeling of drug release from solid pharmaceutical systems
Received 15 May 2008
(matrices) for oral delivery.
Received in revised form
The attention was paid to the behavior of matrices made of hydrogels and drug, and the
19 August 2008
modeling was devoted to reproduce all the relevant phenomena (water up-take, gel swelling,
Accepted 22 October 2008
diffusivity increase, drug diffusion and polymer erosion). Thus, the transient mass balances
Published on line 30 October 2008
(for both drug and water), with the proper initial and boundary conditions were written, and a generalized numerical code was formulated; it is able to describe several geometries
Keywords:
(slab, sphere, infinite and finite cylinders; this latter was done by an approximation which
Drug release
reduces the 2D problem to an 1D scheme). The main phenomena observed in drug delivery
Modeling
from hydrogel-based matrix, i.e. polymer swelling and erosion, were taken into account.
Matrix geometry
The code was validated by comparison with analytical solutions, available for some simplified situation, and then it was tested with some experimental data taken from literature. © 2008 Elsevier B.V. All rights reserved.
1.
Introduction
The controlled release of drugs is a key topic in modern pharmacotherapeutics. Pharmaceutical systems able to release the dose of drug with a given kinetics are highly desired. In particular, even if not only, the so-called zero-order kinetics, i.e. the release of drug with a constant rate, to equate the drug consumption rate due to the metabolism is a valuable goal. Matrices made of hydrogels seem to be able to fulfill this requirement, i.e. properly designed matrix (in terms of composition, as well as dimensions and geometry) can give tailored drug release profiles. Indeed, solid pharmaceutical systems (matrices) for oral administration are the most widespread method to assume drugs, usually they were produced by compression or by granulation of a blend made of the drug (each in crystalline
∗
Corresponding author. Tel.: +39 089964077; fax: +39 089964057. E-mail address: [email protected] (G. Lamberti). 0928-0987/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ejps.2008.10.006
or in amorphous form) and of several excipients, generally polymers (mainly amorphous). The design of these systems requires a deep knowledge of the phenomena involved during matrix hydration/dissolution in the stomach and in the gastrointestinal tract. The phenomena can be described in term of mathematical equations (modeling) which can be solved by a properly designed software (numerical code). The mathematical modeling of drug release can indeed significantly facilitate the development of new and the optimization of existing pharmaceutical products. Certainly, the identification of the physical model and its mathematical description, in a way adequate to describe the real behavior of a pharmaceutical system can be used to simulate the effect of the device design parameters on the release kinetics. Availability of the right model allows to predict a priori the system formulation parameters to obtain the desired drug release kinetics.
360
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
A detailed overview of the modeling approaches followed in literature to describe the drug release from hydrogelbased matrices was presented by Siepmann and Peppas (2001). They started from the well known Higuchi treatment (Higuchi, 1961), developed to describe the drug release from an ointment containing a suspended drug and in contact with a perfect sink. The Higuchi equation is likely the most famous and the most often used tool to describe drug release processes. Basically, it is the solution to the transient mass transfer phenomenon (the Fick’s balance equation), obtained under some simplifications: perfect sink at the ointment/dissolution medium interface, constant diffusivity, no swelling and no erosion of the matrix. Thus, it is not directly applicable to more complex systems such as matrices made of hydrogels, which are subject to swelling and erosion, showing a water-content sensible diffusivity. The solution of Higuchi equation gives a drug release from a slab which is proportional to the square root of time. Since the experimental behavior showed drug release proportional to different powers of time, Peppas and Sahlin (1989) proposed a semi-empirical model in which the drug release was proportional to the sum of two different powers of time. The two terms accounted for contributions of the pure diffusivity on one side (the so-called Fickian transport) and of relaxation of the polymer molecules on the other side, the so-called non-Fickian or case-II “relaxational” contribution. The problem of the drug diffusion in more complex geometries was faced out analytically by Fu et al. (1976), who solved the diffusion problem in a finite cylinder in the absence of swelling and erosion, and under the hypothesis of constant diffusivities. More recently, Zhou and Wu (1997), Wu and Zhou (1998, 1999) and Frenning et al. (2005), applied the finite element method to solve the diffusion problem in matrix of various geometries, even not too simple (including convex matrices, hollow cylinders, doughnuts, inwards hemispheres) (Zhou and Wu, 1997), also in presence of moving boundaries (Wu and Zhou, 1999), or of slowly dissolving drugs (Frenning et al., 2005). However, this approach was, to our knowledge, never applied to swelling and eroding matrices. Kill and Dam-Johansen (2003) proposed a model to describe water up-take, polymer swelling and the drug diffusion in cylindrical matrices exposed to water only by their lateral surfaces. They validated the model by comparison with the experimental data presented by Bettini et al. (2001), related to matrices made of HPMC (hydroxypropylmethylcellulose) and BPP (buflomedilpyridoxal phosphate) confined between two Plexiglass discs and subjected to hydration. The model (in excellent agreement with the experimental data) the movements of swelling and erosion fronts, and the diffusion of drug and water in the swollen layer. The hypotheses: (i) the swelling takes place when water concentration overcomes a threshold value; (ii) the erosion phenomenon was negligible; (iii) the swollen gel was practically drug-free; and (iv) the un-swollen core was waterfree (dry). However the two last hypotheses require, more than the others, an experimental confirmation: in fact we obtained (Chirico et al., 2007a) results not in agreement with them.
Grassi et al. (2004, 2000) and Grassi and Grassi (2005) proposed several approaches, with increasing complexity, to model the drug release from solid pharmaceutical forms. A simple model based on the drug balance in the dissolution medium was developed taking into account the resistance to the release due to a layer of enteric coating (Grassi et al., 2004). A much more complex model was then developed to describe the release from matrices made of swellable hydrogels, with spherical drug particles, poly-dispersed in size and in different physical states (amorphous and crystalline). The water and drug fluxes due to diffusion and viscoelastic effects joint to the swelling phenomenon were described and the volume increase due to swelling was accounted for layer by layer (i.e. dividing the spheres in shells). The model was completed by the drug balance in the dissolution medium (Grassi et al., 2000; Grassi and Grassi, 2005). Along more than a decade, Peppas and co-workers (Korsmeyer et al., 1986a,b; Siepmann et al., 1999a,b, 2000; Narasimhan and Peppas, 1996a,b; Siepmann and Peppas, 2000) developed a comprehensive mathematical theory able to describe the full set of observed experimental phenomena during the drug release from matrices made of swellable hydrogels. Firstly, Korsmeyer et al. (1986a,b) modeled the drug release and the water up-take in thin films of hydrogels, accounting for any kind of water concentration-dependent diffusivity (the model requires numerical solution); then Siepmann, Peppas and co-worker (Siepmann et al., 1999a,b, 2000; Narasimhan and Peppas, 1996a,b; Siepmann and Peppas, 2000) extended the approach to cylindrical matrices, accounting also for swelling and erosion phenomena. The full model was named “sequential layer”, since it describes all the phenomena (water diffusion, swelling, drug diffusion, polymer erosion), layer-by-layer, outer to inner part of the matrix. Even if all the phenomena have received attention by the researchers, under both the experimental and the modeling point of view, it would be nice:
(i) To have a full model (and the related numerical code), useful to describe the drug release from matrices of different shape (slabs, disks or short cylinders, needles or long cylinders, spheres). (ii) To confirm the nature and the extension of the phenomena postulated in describing the drug release from hydrogel-based matrices. Traditional experiments (measurements of drug release by USP apparatuses) are not sufficient to observe and quantify swelling, erosion and diffusion phenomena which were hypothesized to take place during the drug release. These (experimental) issues are currently under investigation in our laboratory (Chirico et al., 2007a), and will be the subjects of future communications.
Aim of this work is thus to formulate a physical model, and the related numerical code, able to describe all the phenomena observed during the drug release from different shaped matrices made of swellable/erodible polymer matrices.
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
2.
Modeling
2.1.
Balance equations
Barba and Lamberti (2003) proposed a generalization of the differential balance equation, and they built a code able to solve transient balance of mass and/or energy across slabs, cylinders of infinite length and spheres. The geometrical generalization consists in the introduction of a geometrical factor f which is 0 for slabs, 1 for infinite cylinders and 2 for spheres. Accounting for variable density and diffusivity, the generalized balance equation reads as follows:
∂ω f ∂2 ω ∂(D) = D 2 + D + ∂t ∂ ∂
∂ω ∂
(1)
where is the density, D is the diffusivity (or a pseudodiffusion coefficient if the balance is written for a system different from two gases or vapours), ω is the mass fraction (dependent variable), t is the time (independent variable), is the position across the body (independent variable, which reads as the thickness direction within a slab, and the radial direction within a cylinder or a sphere) and f is the geometrical factor mentioned above. Here the code was improved by adding the ability to simulate also the behavior of a finite cylinder. Following the reasoning suggested by Coviello et al. (2003), which is reported in Appendix A, the release from a finite cylinder (with equal to the semi-height/radius ratio, = z/r) was described in analogy with the release from a sphere (f = 2), and using a further geometrical factor (), given by Eq. (A4) as (22 + 1)/32 . The generalized balance equation becomes ∂ω i = ∂t
∂2 ω f ∂(Di ) Di 2 i + Di + ∂ ∂
∂ωi ∂
()
(2)
where i = 1 is for the water and i = 2 is for the drug (thus, i = 3 is for the polymer). The balance equations require initial and boundary conditions. They are summarized as follows:
⎧ ⎪ ⎪ (a) I.C. ⎪ ⎪ ⎪ ⎨ (b) B.C.1
@t = 0, ∀ @ = 0, ∀t
⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (c) B.C.2 @ = RE (t), ∀t
ω1 = ω1, in ω2 = ω2, in
∂ω ∂ω2 1 = 0, =0 ∂ ∂ ω1 = ω1, eq ω2 = 0
outer medium, ω1,eq , and the drug concentration in the dissolution medium is negligible, so that its equilibrium value is ω2,eq = 0 (perfect sink condition). The numerical nature of the solution scheme adopted allows for the use of different initial/boundary conditions (however, along this work they were the ones listed in Eq. (3)). It means that, in principle, the initial contents of the water and of the drug in the matrix can be position dependent (ω1,in = f1 () and/or ω2,in = f2 ()), as well as the boundary conditions can account for time dependent concentration (ω1 (t, = RE ) = g1 (t) and/or ω2 (t, = RE ) = g2 (t)). Furthermore, the code allows different kind of boundary conditions (e.g., the mass transfer at interface instead of a given value of concentration, see Barba and Lamberti (2003) for further details).
2.2.
Constitutive equations
To solve Eq. (2), the pseudo-diffusion coefficients, Di (for i = 1, 2), have to be evaluated. In polymeric systems subjected to swelling, the diffusion coefficients are not constant, being low in the dry polymer and increasing as the water content increases (into the gel). They can be modeled according to Siepmann et al. (1999b):
Di (ω1 ) =
The initial condition (3a) states that the matrix initially contains water at mass fraction ω1,in (if the matrices were initially dry, ω1,in = 0) and drug at mass fraction ω2,in ; the axial boundary condition (3b) is a standard symmetry (no flux) condition for both the components (water and drug); this condition holds at the slab mid-plane, as well as at the cylinder axis and at the sphere centrum, all these conditions are in fact described by = 0; the surface condition (3c) takes into account that, at the external boundary, identified by means of an “equivalent” erosion radius (RE , which holds exactly for a large height/radius cylinder and for a sphere, and which is the semi-thickness for a slab) the water mass fraction is in equilibrium with the
D∗i
exp −ˇi
ω1 1− ω1, eq
(4)
where D∗i /exp(ˇi ) is the value of the pseudo-diffusion coefficient in the dry matrix (ω1 = 0), and D∗i is the value of the pseudo-diffusion coefficient in the fully swollen matrix (ω1 = ω1,eq ).
2.3.
Solution procedure
Since the pseudo-diffusion coefficients are concentrationdependent, and since the outer boundary is moving, RE = RE (t) and the Eq. (2) has to be solved numerically. The erosion front position RE (t) was calculated on the basis of the matrix’s mass and density. This step thus requires the knowledge of water and polymer mass in the matrix at any given time instant, t.
2.3.1. (3)
361
Water mass calculation
At each time step, the resolution of the diffusion equation (2) gives the water mass fraction profiles, ωi (t, r) (for i = 1, 2). A proper integration of this profile gives the actual value of the water and drug masses (mi (t), by Eq. (5) below, as well as the value of the water and drug mass fractions averaged on the matrix volume, ωi (t), by Eq. (6)):
mi (t) =
(ω1 , ω2 )ωi (t, r) dV
(5)
V
ωi (t) =
2.3.2.
1 V
ωi (t, r) dV
(6)
V
Polymer mass calculation
The polymer mass, m3 (t), is given by the resolution of a balance equation which accounts for the erosion phenomenon, postulated to be proportional to the surface area, AE (t) (to
362
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
our knowledge, this is the simplest and most widely used approach (Siepmann et al., 1999a), in which the erosion is considered as a first-order surface phenomena):
dm3 = −ke AE (t) dt m3 (t = 0) = m3, in
(7)
At each time step the erosion radius, RE , and thus the erosion surface, AE , are known, allowing for the calculation of the polymer mass, m3 (t), by solving Eq. (7), once a proper value for erosion constant, kE , is given.
2.3.3.
1 ω1 ω2 1 − ω1 − ω2 = + + tot 1 2 3
(8)
where 1 , 2 and 3 are the water, the drug and the polymer densities, respectively. Once they are known, the drug and the polymer masses are known, the current volume of the matrix can be calculated: mtot (t) tot
(9)
Finally, the required value for the erosion radius was obtained by the volume value.
2.3.4.
Geometrical issues
In Eq. (5), (6) and (9) the matrix volume, V, is involved; similarly, on the other hand, in Eq. (7) the external (erosion) surface, AE , has to be known. Recalling the meaning of the geometrical factor f, which is 0 for slabs, 1 for infinite cylinders, 2 for spheres and for finite cylinder, one can write: AE (t) =
f CA RE (t)
(10)
Results and discussions
3.1.
Preliminary validation
The code’s core was validated elsewhere (Barba and Lamberti, 2003) by comparison with known analytical solutions. In this work but here the code was largely modified, thus a new step of preliminary validation was required. Therefore, the code was used to simulate three cases for which the analytical solutions are available: the release from (i) a slab of constant thickness, (ii) an “infinite” cylinder (a cylinder with a large ratio height/radius) of constant radius, (iii) a sphere of constant radius (no moving boundaries means constant density, thus the matrix does not “swell” because of water up-take and no erosion, and the polymer does not dissolve in the medium). For these three cases the analytical solutions are known, under the further restriction of constant diffusivities, and are reported in Table 2. The solutions were obtained by solving Eq. (2) with initial and boundary conditions given by Eq. (3). For each case, Table 2 reports the normalized mass fraction profile of the ith component, ϕi (t, ), as a function of time, t, and of the generic position variable, (which is the axial direction in a slab, and the radial direction in an “infinite” cylinder and in a sphere). Table 2 also reports the cumulative drug release of the ith component, DRi (t). The normalized profile is defined as ϕi (t, ) =
ωi,in − ω(t, ) ωi,in − ωi,eq
(12)
And the cumulative release is defined as the ratio of the ith component mass released at the time t over its initial value:
and f +1
V(t) = CV RE
(t)
(11)
The use of coefficients CA and CV allows to reproduce the desired shape of the matrix, leading to the code generalization. Their values are summarized in Table 1. Thus, during the solution procedure, at each time step the code evaluates: the drug and water masses (Eq. (5)) (and by this way it accounts for the swelling), the polymer mass (Eq. (7)) (and by this way it accounts for the erosion), the overall density (Eq. (8)), the matrix volume (Eq. (9)), and at last the “erosion radius”, i.e. the actual size of the matrix (Eq. (11)). The erosion radius change (increase due to swelling and/or decrease due to the erosion) requires of course an update of the spatial grid, which is performed for each time step.
2.3.5.
3.
Erosion radius calculation
The average density of the partially hydrated matrix can be calculated from the simplest mixing rule which can be written for the specific volume:
V(t) =
equation (2), and the Euler finite difference scheme for the erosion equation (7), the integral equations ((5) and (6)) being approximated by summations. The code was developed using the Mathcad© software (Mathsoft Engineering & Education, Inc.). To ensure numerical solution stability, each simulation was carried out several times by doubling time and space step numbers, until the obtained solution become independent from numbers of these steps.
Numerical issues
The model equations were solved numerically, adopting the Crank & Nicolson finite difference scheme for the diffusion
DRi (t) = 1 −
V(t)
ωi (t, ) dV
ωi (0, ) dV V0
(13)
The code was validated by comparison of its calculations with analytical solutions, obtained assigning to the model parameters reasonable values, based on literature or physical data. Here the following values were used: ω1,in = 0.00 (initial water content in the matrix), ω2,in = 0.60 (initial drug content in the matrix), ω1,eq = 0.80 (equilibrium water mass fraction), ω2,eq = 0.00 (equilibrium drug mass fraction), D1 = 5.6 × 10−9 m2 s−1 (water diffusivity within the matrix), D2 = 12 × 10−9 m2 s−1 (drug diffusivity within the matrix) and RE = 6 mm (slab thickness or cylinder radius or sphere radius). Results of the calculations are reported in Fig. 1, in which the analytical solutions and the code predictions perfectly agree. From these comparisons the code can be considered as validated.
363
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Table 1 – Geometrical factors for code generalization. Geometry
f
CA
CV
Slab (square of side L) Infinite cylinder (height H RE ) Sphere Finite cylinder with constant
0 1 2 2
Not of use Not of use Not of use H/2RE
1 1 1 (22 + 1)/32
2L2 2H 4 2(1 + 2)
2L2 H 4/3 2
3.2.
Modeling for a finite cylinder
der height) and for > 3 the cylinder can be taken as “infinite” (the height being considerably larger than the radius). The values of the correction factor s which causes the numerical solutions to be in agreement with the analytical ones are reported as symbols in Fig. 3 (left). These data were also fitted by a suitable relationship, i.e.
Following the reasonings proposed by Coviello et al. (2003), as summarized in Appendix A, which lead to Eq. (2), the code should be able to describe the behavior of a “finite” cylinder, i.e. a cylinder with comparable height and radius. The diffusion equation in a cylinder (which would arise from Eq. (A2)) is different from the diffusion equation in a sphere (see Eq. (A4), the basis for Eq. (2)). This means that, to make the two solutions coincident (the solution of the 2D diffusion problem into a cylinder of radius RE and the solution of the 1D diffusion problem into a sphere of radius ()RE ), the code needs for a correction coefficient. Actually, the simulations of a finite cylinder with = 2 (no swelling–no erosion – constant diffusivities), obtained by the analytical solution and by the numerical code (Eq. (2)), do not give the same results (see Fig. 2, left). To make the code able to reproduce the analytical result, the right hand side of Eq. (2) has to be multiplied by a factor which is roughly s = 0.17. By using this value, the numerical predictions match the analytical ones (see Fig. 2, right). Then, the same procedure can be followed for several values of the ratio between height and diameter, = height/diameter = semi-height/radius. There is no need to investigate very low and very large values of , since for < 0.2 the cylinder is practically a slab (the thickness being the cylin-
2 A 2 e−(ln(/c )) /2ω s() = s0 + √ 2 w
(14)
in which s0 = 4.32 × 10−3 , lc = 1.7745, w = 0.8769 and A = 0.7652. It is worth to be noted that Eq. (14) does not have any physical meanings, it only assures that the full model, by using f = 2 (i.e. for a sphere) and given by Eq. (A4), can be used in the description of a finite cylinder. In fact, the correction factor given by Eq. (14) allows the code to reproduce the behavior of a finite cylinder (with 0.2 < < 3). Cylinders with values of lower than 0.2 can be simulated as slabs, and cylinders with values of higher than 3 can be simulated as infinite cylinders. Fig. 3 (the right) also reports the simulation obtained by using the numerical code, with the correction factor s given by Eq. (14), compared with the analytical solution, in term of drug release with time. Summarizing, the code given by Eqs. (2)–(11), including the correction coefficient given by Eq. (14) (which holds only for “finite” cylinders), is able to reproduce the diffusion of both
Table 2 – Summary of the equations which hold for constant diffusivity – no moving boundary, and for several geometries. Slab (square of side L and thickness 2RE ) (§ 4.3.2, p. 47, Crank, 1975)
4
ϕi,slab (t, ) = 1 −
∞
∞
ϕi,cylinder (t, ) = 1 −
2 RE
(2n+1)2 2 Di t 4 R2 E
8 (2n+1)2 2
n=0
∞
exp −
J0 (˛n ) ˛n J1 (RE ˛n )
∞
4 R2 ˛2 E n
cos
(2n+1)2 2 Di t 4 R2 E
(2n+1) 2
RE
exp(−Di ˛2n t) where ˛n ’s are the positive roots of J0 (RE ˛n ) = 0
⎧ n=1 ⎪ ⎪ ⎪ ⎨ =0 ϕi,sphere (t, ) =
exp(−Di ˛2n t)
1+2
∞
⎪ ⎪ ⎪ ⎩>0
DRi,sphere (t) = 1 −
6 2
n
2 2
(−1) exp −n
2RE 1+
∞ n (−1)
n n=1
1 n2
exp
−n2 2
n=1
ϕi,finite
cylinder (t, )
DRi,finite
cylinder (t)
Di t R2 E
= ϕi,slab(t, )ϕi,cylinder (t, ) RE
=1−
n=1
∞
Finite cylinder (height H and radius RE , = H/2RE )
n=1
DRi,cylinder (t) = 1 −
Sphere (of radius RE ) (§ 6.3.1, p. 90, Crank, 1975)
exp −
n=0
DRi,slab (t) = 1 −
Infinite cylinder (height H RE ) (§ 5.3.1, p. 73, Crank, 1975)
(−1)n 2n+1
0RE 0
f dRE E f ωi,in ()CV R dRE E ωi (t,)CV R
Di t RE2 2 2
exp −n
Di t RE2
sin
n
RE
364
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Fig. 1 – (Left) Drug and water mass fraction profiles within the matrix (after a given time from immersion). (Right) Evolution with time of the drug release. Symbols: results of numerical code; lines: results of analytical solutions as reported in Table 2. From the top: results for a slab, an “infinite” cylinder and a sphere.
water and drug within the solid pharmaceutical form (“the matrix”), for several shapes: slabs, infinite cylinders, sphere and finite cylinders. The code’s ability was tested by comparison with the analytical solutions, which are available under some restrictions (no swelling–no erosion – constant diffusivities).
3.3.
Case history from literature
The last, most important, challenge for the code is to test its ability in predicting a real case. Here we take a real case from Siepmann et al. (1999b). Cylindrical matrices made of HPMC (a swellable hydrogel commonly used in controlled release
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
365
Fig. 2 – Analytical and numerical drug and water profiles for the diffusion within a finite cylinder with = 2. The analytical solution was obtained as reported in Table 2; the numerical solution was obtained by solving Eq. (2). (Left) The numerical solutions came from Eq. (2) without any correction. (Right) The numerical solutions came from Eq. (2) using a correction factor s().
systems) and of propanolol HCl (a non-selective beta blocker, mainly used in the treatment of hypertension) were prepared by direct compression of powders mixture (5% (w/w) drug, 28 mg of HPMC and 1.5 mg of propanolol HCl). The matrices were 2.5 mm in radius and 0.7 mm in semi-height. Thus the ratio was 0.28. The matrices were immersed in phosphate buffer (pH 7.4) following the USP rotating paddle method. The drug release was followed by periodical assaying of the release medium (spectroscopically, at wavelength of 290 nm). Experimental results were reported as symbols in Fig. 4 (left). In the same graph, the prediction by a 2D model developed in Siepmann et al. (1999b) was also reported (dashed line). The model, proposed by Siepmann, Peppas and co-worker (Siepmann et al., 1999b), is based on the same assumptions
used in the present work, the main difference being the ability of the code here presented to describe the behavior of different shaped matrices (it has to be noted that the code is mathematically simpler, being 1D). To obtain the simulation in Fig. 4, the parameters to be used are: ω1,eq = 0.802 (in Eqs. (3) and (4)), ˇ1 = 2.5, ˇ2 = 9.5, D∗1 = 56 × 10−11 , D∗2 = 6.3 × 10−11 m2 /s (in Eq. (4)) and ke = 5.5 × 10−6 kg/(m2 s) (in Eq. (7)). The parameters have been taken from Siepmann et al. (1999b), no further optimization were performed in this work. In Fig. 4 the predictions of the model presented in this work (with the same values of the parameters) were reported as a solid line. The predictions of the two codes are practically coincident, confirming the ability of the proposed code
Fig. 3 – (Left) The values of the correction factor s(), as obtained by comparison between numerical and analytical solutions, together with a fitting curve. (Right) the cumulative drug release, in the case of = 2, as predicted by analytical solution (Eq. (13)) and by numerical solution (Eq. (2)), using the correction factor s().
366
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Fig. 4 – (Left) The drug release in phosphate buffer (pH 7.4) from a cylindrical matrix of hydroxypropylmethylcellulose (28 mg)/propanolol HCl (1.5 mg), with = 0.28. Experimental data and 2D model predictions came from Siepmann et al. (1999b); the 1D model predictions was obtained in this work. (Right) The time evolutions of water, drug and polymer masses, as predicted by the code.
to reproduce both the experimental data and the predictions of a more complex and less versatile model, only applicable to cylindrical matrices. Fig. 4 reports also the predictions in term of mass evolutions for drug, water and polymer. These evolutions were simulated by the code, and the expected behavior can be easily recognized: the drug and the polymer masses show a monotonically decreasing behavior (the drug diffuses toward the dissolution medium, the polymer is subjected to erosion). The water content firstly increases (due to the up-take from the external medium), and then decreases (due to the erosion which causes the external layers of gel to be destroyed). To be sure that the predicted phenomena really takes place, and when they take place, a comparison with experimental data is required, i.e. the availability of experimental data to be compared to above predictions would be a very powerful piece of information. Currently, our research group is working to define the experimental procedures necessary to obtain these data.
Fig. A1 – The scheme for differential volume identification in finite cylinder analysis.
4.
Conclusions
Controlled release of drugs can be obtained by using matrices made of swellable hydrogels. The availability of a mathematical tool, able to predict the drug release kinetics from matrices of several shapes, dimensions and made of several materials, is of great aid in design and manufacture of new pharmaceutical systems. In this work, a physical model was formulated, to describe all the phenomena commonly observed during hydration of hydrogels made matrices. Water up-take, which causes the hydrogels swelling and the (water’s and drug’s) diffusivities increase, polymer erosion and drug release were mathematically described. The equations were implemented in a numerical code, based on a scheme previously proposed. The code was built on the basis of a generalization (Barba and Lamberti, 2003) which made it able to predict the behavior of matrix with several simple shapes: an infinite slab (thickness negligible with respect to the other dimensions), an infinite cylinder (radius negligible with respect to the height) and a sphere. The first step, was to validate the code by comparing the results with some analytical solutions known for simple cases (constant diffusivity, no swelling–no erosion). Then, on the basis of reasonings suggested in Coviello et al. (2003), the code was made able to describe also the behavior of a finite cylinder, with constant height/radius ratio. The comparison of the modified code with analytical solutions for this case allowed to add a correction factor to the code itself. Finally, the full code was used to simulate the behavior of a real system, reported in Siepmann et al. (1999b). The simplified (1D+time) code proposed here was able to reproduce the experimental data with the same accuracy of the more complex model (2D+time) available in literature. Summarizing, in this work a model suitable to describe the drug release from different shaped matrices made of swellable hydrogels, was formulated and implemented into a numerical
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
code. The code was validated by comparison with analytical solutions, and it was found able to reproduce a real case from literature. Properly designed experiments are highly desirable to further test the code ability to describe real systems. These experiments have to be designed to gather detailed information on the phenomena involved in the drug release (water up-take, matrix swelling, polymer erosion, drug diffusion). Such experiments are currently under development in our laboratory (see for example Chirico et al., 2007b). The code proposed here could be of use both as an aid in quantification of the phenomena observed in R&D studies, and useful in design and manufacture of new pharmaceutical systems.
Appendix A The basic idea (Coviello et al., 2003) to get a simplified description (2D → 1D) of a finite cylinder is to consider a cylindrical shell with a semi-height/radius ratio = z/r equal to its initial value, = Z0 /R0 , and to take this shell as the computational element (Fig. A1). The differential volume of the computational element is 2
V = [(r + r) (z + z) − r2 z] = 6 r2 r
(A1)
The mass balance within the computational volume can be written on the basis of constant concentration in the shell as
V
∂ω = (2 r2znr ) − (2 r2znr ) r r+ r ∂t
+ (2 r2 nz ) − (2 r2 nz ) z
z+ z
(A2)
The axial flux can be related to the radial flux as follows: nz = −D
∂ω D ∂ω 1 = − ∂r = nr ∂z z r
(A3)
Insertion of (A3) in (A2), followed by dividing by V, recalling that = z/r and taking the limit 0, gives ∂ω ∂t = 2 for r1 → 4 ∂ 2 2 ∂ 2 1 ∂ 2 − 6 r2 ∂r (r nr ) − 6 2 r2 ∂r (r nr ) = − 3 + 32 r2 ∂r (r nr )or
∂ω =− ∂t
22 + 1 32
1 ∂ 2 1 ∂ (r nr ) = −() 2 (r2 nr ) r2 ∂r r ∂r
(A4)
Eq. (A4) is similar to the (1D) mass balance within a sphere, but for the correction coefficient (), for which Eq. (A4) is the definition.
references
Barba, A.A., Lamberti, G., 2003. Preliminary validation of a numerical code for heat transfer simulations. Heat Mass Transfer 39, 429–433. Bettini, R., Catellani, P.L., Santi, P., Massimo, G., Peppas, N.A., Colombo, P., 2001. Translocation of drug particles in HPMC matrix gel layer: effect of drug solubility and influence on release rate. J. Control. Rel. 70, 383–391.
367
Chirico, S., Lamberti, G., Nunziata, V., Titomanlio, G., 2007a. Analysis and modeling of radial water up-take in pure HPMC matrices. In: Proceedings of the European Congress of Chemical Engineering (ECCE-6), Copenhagen, September. Chirico, S., Dalmoro, A., Lamberti, G., Russo, G., Titomanlio, G., 2007b. Analysis and modeling of swelling and erosion behavior for pure HPMC matrix. J. Control. Rel. 122, 181–188. Coviello, T., Grassi, M., Lapasin, R., Marino, A., Alhaique, F., 2003. Scleroglucan/borax: characterization of a novel hydrogel system suitable for drug delivery. Biomaterials 24, 2789–2798. Crank, J., 1975. The Mathematic of Diffusion. Clarendon Press, Oxford. Frenning, G., Brohede, U., Strømme, M., 2005. Finite element analysis of the release of slowly dissolving drugs from cylindrical matrix systems. J. Control. Rel. 107, 320–329. Fu, J.C., Hagemeir, C., Moyer, D.L., Ng, E.W., 1976. A unified mathematical model for diffusion from drug–polymer composite matrices. J. Biomed. Mater. Res. 10, 743–758. Grassi, M., Grassi, G., 2005. Mathematical modelling and controlled drug delivery: matrix systems. Curr. Drug Deliv. 2, 97–116. Grassi, M., Colombo, I., Lapasin, R., 2000. Drug release from an ensemble of swellable crosslinked polymer particles. J. Control. Rel. 68, 97–113. Grassi, M., Zema, L., Sangalli, M.E., Maroni, A., Giordano, F., Gazzaniga, A., 2004. Modeling of drug release from partially coated matrices made of a high viscosity HPMC. Int. J. Pharm. 276, 107–114. Higuchi, T., 1961. Rate of release of medicaments from ointment bases containing drugs in suspensions. J. Pharm. Sci. 50, 874–875. Kill, S., Dam-Johansen, D., 2003. Controlled drug delivery from swellable hydroxypropylmethylcellulose matrices: model-based analysis of observed radial front movements. J. Control. Rel. 90, 1–21. Korsmeyer, R.W., Lustig, S.R., Peppas, N.A., 1986a. Solute and penetrant diffusion in swellable polymers. I. Mathematical modeling. J. Polym. Sci. B: Polym. Phys. 24 (2), 395–408. Korsmeyer, R.W., von Meerwall, E., Peppas, N.A., 1986b. Solute and penetrant diffusion in swellable polymers. II. Verification of theoretical models. J. Polym. Sci. Polym. Phys. Ed. 24, 409–434. Narasimhan, B., Peppas, N.A., 1996a. Disentanglement and reptation during dissolution of rubbery polymers. J. Polym. Sci. Polym. Phys. 34, 947–961. Narasimhan, B., Peppas, N.A., 1996b. On the importance of chain reptation during dissolution of rubbery polymers. Macromolecules 29, 3283–3291. Peppas, N.A., Sahlin, J.J., 1989. A simple equation for the description of solute release. III. Coupling of diffusion and relaxation. Int. J. Pharm. 57, 169–172. Siepmann, J., Peppas, N.A., 2000. Hydrophilic matrices for controlled drug delivery: an improved mathematical model to predict the resulting drug release kinetics (the “Sequential Layer” model). Pharm. Res. 17, 1290–1298. Siepmann, J., Peppas, N.A., 2001. Modeling of drug release from delivery systems based on hydroxypropyl methylcellulose (HPMC). Adv. Drug Deliv. Rev. 48, 139–157. Siepmann, J., Poudal, K., Sriwongjanya, M., Peppas, N.A., Bodmeier, R., 1999a. A new model describing the swelling and drug release kinetics from hydroxypropyl methylcellulose matrices. J. Pharm. Sci. 88, 65–72. Siepmann, J., Kranz, H., Bodmeier, R., Peppas, N.A., 1999b. HPMC-matrices for controlled drug delivery: a new model combining diffusion, swelling, and dissolution mechanism and predicting the release kinetics. Pharm. Res. 16, 1748–1756.
368
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Siepmann, J., Kranz, H., Peppas, N.A., Bodmeier, R., 2000. Calculation of the required size and shape of hydroxypropyl methylcellulose matrices to achieve desired drug release profiles. Int. J. Sci. 201, 151–164. Wu, X.Y., Zhou, Y., 1998. Finite element analysis of diffusional drug release from complex matrix systems. II. Factors influencing release kinetics. J. Control. Rel. 51, 57–71.
Wu, X.Y., Zhou, Y., 1999. Studies of diffusional release of a dispersed solute from polymeric matrixes by finite element method. J. Pharm. Sci. 88, 1050–1057. Zhou, Y., Wu, X.Y., 1997. Finite element analysis of diffusional drug release from complex matrix systems. I. Complex geometries and composite structures. J. Control. Rel. 49, 277–288.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ejps
A general code to predict the drug release kinetics from different shaped matrices Anna Angela Barba a , Matteo d’Amore a , Serafina Chirico b , Gaetano Lamberti b,∗ , Giuseppe Titomanlio b a b
Dipartimento di Scienze Farmaceutiche, Università degli Studi di Salerno, Via Ponte don Melillo, I-84084 Fisciano (SA), Italy Dipartimento di Ingegneria Chimica e Alimentare, Università degli Studi di Salerno, Via Ponte don Melillo, I-84084 Fisciano (SA), Italy
a r t i c l e
i n f o
a b s t r a c t
Article history:
This work deals with the modeling of drug release from solid pharmaceutical systems
Received 15 May 2008
(matrices) for oral delivery.
Received in revised form
The attention was paid to the behavior of matrices made of hydrogels and drug, and the
19 August 2008
modeling was devoted to reproduce all the relevant phenomena (water up-take, gel swelling,
Accepted 22 October 2008
diffusivity increase, drug diffusion and polymer erosion). Thus, the transient mass balances
Published on line 30 October 2008
(for both drug and water), with the proper initial and boundary conditions were written, and a generalized numerical code was formulated; it is able to describe several geometries
Keywords:
(slab, sphere, infinite and finite cylinders; this latter was done by an approximation which
Drug release
reduces the 2D problem to an 1D scheme). The main phenomena observed in drug delivery
Modeling
from hydrogel-based matrix, i.e. polymer swelling and erosion, were taken into account.
Matrix geometry
The code was validated by comparison with analytical solutions, available for some simplified situation, and then it was tested with some experimental data taken from literature. © 2008 Elsevier B.V. All rights reserved.
1.
Introduction
The controlled release of drugs is a key topic in modern pharmacotherapeutics. Pharmaceutical systems able to release the dose of drug with a given kinetics are highly desired. In particular, even if not only, the so-called zero-order kinetics, i.e. the release of drug with a constant rate, to equate the drug consumption rate due to the metabolism is a valuable goal. Matrices made of hydrogels seem to be able to fulfill this requirement, i.e. properly designed matrix (in terms of composition, as well as dimensions and geometry) can give tailored drug release profiles. Indeed, solid pharmaceutical systems (matrices) for oral administration are the most widespread method to assume drugs, usually they were produced by compression or by granulation of a blend made of the drug (each in crystalline
∗
Corresponding author. Tel.: +39 089964077; fax: +39 089964057. E-mail address: [email protected] (G. Lamberti). 0928-0987/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ejps.2008.10.006
or in amorphous form) and of several excipients, generally polymers (mainly amorphous). The design of these systems requires a deep knowledge of the phenomena involved during matrix hydration/dissolution in the stomach and in the gastrointestinal tract. The phenomena can be described in term of mathematical equations (modeling) which can be solved by a properly designed software (numerical code). The mathematical modeling of drug release can indeed significantly facilitate the development of new and the optimization of existing pharmaceutical products. Certainly, the identification of the physical model and its mathematical description, in a way adequate to describe the real behavior of a pharmaceutical system can be used to simulate the effect of the device design parameters on the release kinetics. Availability of the right model allows to predict a priori the system formulation parameters to obtain the desired drug release kinetics.
360
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
A detailed overview of the modeling approaches followed in literature to describe the drug release from hydrogelbased matrices was presented by Siepmann and Peppas (2001). They started from the well known Higuchi treatment (Higuchi, 1961), developed to describe the drug release from an ointment containing a suspended drug and in contact with a perfect sink. The Higuchi equation is likely the most famous and the most often used tool to describe drug release processes. Basically, it is the solution to the transient mass transfer phenomenon (the Fick’s balance equation), obtained under some simplifications: perfect sink at the ointment/dissolution medium interface, constant diffusivity, no swelling and no erosion of the matrix. Thus, it is not directly applicable to more complex systems such as matrices made of hydrogels, which are subject to swelling and erosion, showing a water-content sensible diffusivity. The solution of Higuchi equation gives a drug release from a slab which is proportional to the square root of time. Since the experimental behavior showed drug release proportional to different powers of time, Peppas and Sahlin (1989) proposed a semi-empirical model in which the drug release was proportional to the sum of two different powers of time. The two terms accounted for contributions of the pure diffusivity on one side (the so-called Fickian transport) and of relaxation of the polymer molecules on the other side, the so-called non-Fickian or case-II “relaxational” contribution. The problem of the drug diffusion in more complex geometries was faced out analytically by Fu et al. (1976), who solved the diffusion problem in a finite cylinder in the absence of swelling and erosion, and under the hypothesis of constant diffusivities. More recently, Zhou and Wu (1997), Wu and Zhou (1998, 1999) and Frenning et al. (2005), applied the finite element method to solve the diffusion problem in matrix of various geometries, even not too simple (including convex matrices, hollow cylinders, doughnuts, inwards hemispheres) (Zhou and Wu, 1997), also in presence of moving boundaries (Wu and Zhou, 1999), or of slowly dissolving drugs (Frenning et al., 2005). However, this approach was, to our knowledge, never applied to swelling and eroding matrices. Kill and Dam-Johansen (2003) proposed a model to describe water up-take, polymer swelling and the drug diffusion in cylindrical matrices exposed to water only by their lateral surfaces. They validated the model by comparison with the experimental data presented by Bettini et al. (2001), related to matrices made of HPMC (hydroxypropylmethylcellulose) and BPP (buflomedilpyridoxal phosphate) confined between two Plexiglass discs and subjected to hydration. The model (in excellent agreement with the experimental data) the movements of swelling and erosion fronts, and the diffusion of drug and water in the swollen layer. The hypotheses: (i) the swelling takes place when water concentration overcomes a threshold value; (ii) the erosion phenomenon was negligible; (iii) the swollen gel was practically drug-free; and (iv) the un-swollen core was waterfree (dry). However the two last hypotheses require, more than the others, an experimental confirmation: in fact we obtained (Chirico et al., 2007a) results not in agreement with them.
Grassi et al. (2004, 2000) and Grassi and Grassi (2005) proposed several approaches, with increasing complexity, to model the drug release from solid pharmaceutical forms. A simple model based on the drug balance in the dissolution medium was developed taking into account the resistance to the release due to a layer of enteric coating (Grassi et al., 2004). A much more complex model was then developed to describe the release from matrices made of swellable hydrogels, with spherical drug particles, poly-dispersed in size and in different physical states (amorphous and crystalline). The water and drug fluxes due to diffusion and viscoelastic effects joint to the swelling phenomenon were described and the volume increase due to swelling was accounted for layer by layer (i.e. dividing the spheres in shells). The model was completed by the drug balance in the dissolution medium (Grassi et al., 2000; Grassi and Grassi, 2005). Along more than a decade, Peppas and co-workers (Korsmeyer et al., 1986a,b; Siepmann et al., 1999a,b, 2000; Narasimhan and Peppas, 1996a,b; Siepmann and Peppas, 2000) developed a comprehensive mathematical theory able to describe the full set of observed experimental phenomena during the drug release from matrices made of swellable hydrogels. Firstly, Korsmeyer et al. (1986a,b) modeled the drug release and the water up-take in thin films of hydrogels, accounting for any kind of water concentration-dependent diffusivity (the model requires numerical solution); then Siepmann, Peppas and co-worker (Siepmann et al., 1999a,b, 2000; Narasimhan and Peppas, 1996a,b; Siepmann and Peppas, 2000) extended the approach to cylindrical matrices, accounting also for swelling and erosion phenomena. The full model was named “sequential layer”, since it describes all the phenomena (water diffusion, swelling, drug diffusion, polymer erosion), layer-by-layer, outer to inner part of the matrix. Even if all the phenomena have received attention by the researchers, under both the experimental and the modeling point of view, it would be nice:
(i) To have a full model (and the related numerical code), useful to describe the drug release from matrices of different shape (slabs, disks or short cylinders, needles or long cylinders, spheres). (ii) To confirm the nature and the extension of the phenomena postulated in describing the drug release from hydrogel-based matrices. Traditional experiments (measurements of drug release by USP apparatuses) are not sufficient to observe and quantify swelling, erosion and diffusion phenomena which were hypothesized to take place during the drug release. These (experimental) issues are currently under investigation in our laboratory (Chirico et al., 2007a), and will be the subjects of future communications.
Aim of this work is thus to formulate a physical model, and the related numerical code, able to describe all the phenomena observed during the drug release from different shaped matrices made of swellable/erodible polymer matrices.
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
2.
Modeling
2.1.
Balance equations
Barba and Lamberti (2003) proposed a generalization of the differential balance equation, and they built a code able to solve transient balance of mass and/or energy across slabs, cylinders of infinite length and spheres. The geometrical generalization consists in the introduction of a geometrical factor f which is 0 for slabs, 1 for infinite cylinders and 2 for spheres. Accounting for variable density and diffusivity, the generalized balance equation reads as follows:
∂ω f ∂2 ω ∂(D) = D 2 + D + ∂t ∂ ∂
∂ω ∂
(1)
where is the density, D is the diffusivity (or a pseudodiffusion coefficient if the balance is written for a system different from two gases or vapours), ω is the mass fraction (dependent variable), t is the time (independent variable), is the position across the body (independent variable, which reads as the thickness direction within a slab, and the radial direction within a cylinder or a sphere) and f is the geometrical factor mentioned above. Here the code was improved by adding the ability to simulate also the behavior of a finite cylinder. Following the reasoning suggested by Coviello et al. (2003), which is reported in Appendix A, the release from a finite cylinder (with equal to the semi-height/radius ratio, = z/r) was described in analogy with the release from a sphere (f = 2), and using a further geometrical factor (), given by Eq. (A4) as (22 + 1)/32 . The generalized balance equation becomes ∂ω i = ∂t
∂2 ω f ∂(Di ) Di 2 i + Di + ∂ ∂
∂ωi ∂
()
(2)
where i = 1 is for the water and i = 2 is for the drug (thus, i = 3 is for the polymer). The balance equations require initial and boundary conditions. They are summarized as follows:
⎧ ⎪ ⎪ (a) I.C. ⎪ ⎪ ⎪ ⎨ (b) B.C.1
@t = 0, ∀ @ = 0, ∀t
⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (c) B.C.2 @ = RE (t), ∀t
ω1 = ω1, in ω2 = ω2, in
∂ω ∂ω2 1 = 0, =0 ∂ ∂ ω1 = ω1, eq ω2 = 0
outer medium, ω1,eq , and the drug concentration in the dissolution medium is negligible, so that its equilibrium value is ω2,eq = 0 (perfect sink condition). The numerical nature of the solution scheme adopted allows for the use of different initial/boundary conditions (however, along this work they were the ones listed in Eq. (3)). It means that, in principle, the initial contents of the water and of the drug in the matrix can be position dependent (ω1,in = f1 () and/or ω2,in = f2 ()), as well as the boundary conditions can account for time dependent concentration (ω1 (t, = RE ) = g1 (t) and/or ω2 (t, = RE ) = g2 (t)). Furthermore, the code allows different kind of boundary conditions (e.g., the mass transfer at interface instead of a given value of concentration, see Barba and Lamberti (2003) for further details).
2.2.
Constitutive equations
To solve Eq. (2), the pseudo-diffusion coefficients, Di (for i = 1, 2), have to be evaluated. In polymeric systems subjected to swelling, the diffusion coefficients are not constant, being low in the dry polymer and increasing as the water content increases (into the gel). They can be modeled according to Siepmann et al. (1999b):
Di (ω1 ) =
The initial condition (3a) states that the matrix initially contains water at mass fraction ω1,in (if the matrices were initially dry, ω1,in = 0) and drug at mass fraction ω2,in ; the axial boundary condition (3b) is a standard symmetry (no flux) condition for both the components (water and drug); this condition holds at the slab mid-plane, as well as at the cylinder axis and at the sphere centrum, all these conditions are in fact described by = 0; the surface condition (3c) takes into account that, at the external boundary, identified by means of an “equivalent” erosion radius (RE , which holds exactly for a large height/radius cylinder and for a sphere, and which is the semi-thickness for a slab) the water mass fraction is in equilibrium with the
D∗i
exp −ˇi
ω1 1− ω1, eq
(4)
where D∗i /exp(ˇi ) is the value of the pseudo-diffusion coefficient in the dry matrix (ω1 = 0), and D∗i is the value of the pseudo-diffusion coefficient in the fully swollen matrix (ω1 = ω1,eq ).
2.3.
Solution procedure
Since the pseudo-diffusion coefficients are concentrationdependent, and since the outer boundary is moving, RE = RE (t) and the Eq. (2) has to be solved numerically. The erosion front position RE (t) was calculated on the basis of the matrix’s mass and density. This step thus requires the knowledge of water and polymer mass in the matrix at any given time instant, t.
2.3.1. (3)
361
Water mass calculation
At each time step, the resolution of the diffusion equation (2) gives the water mass fraction profiles, ωi (t, r) (for i = 1, 2). A proper integration of this profile gives the actual value of the water and drug masses (mi (t), by Eq. (5) below, as well as the value of the water and drug mass fractions averaged on the matrix volume, ωi (t), by Eq. (6)):
mi (t) =
(ω1 , ω2 )ωi (t, r) dV
(5)
V
ωi (t) =
2.3.2.
1 V
ωi (t, r) dV
(6)
V
Polymer mass calculation
The polymer mass, m3 (t), is given by the resolution of a balance equation which accounts for the erosion phenomenon, postulated to be proportional to the surface area, AE (t) (to
362
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
our knowledge, this is the simplest and most widely used approach (Siepmann et al., 1999a), in which the erosion is considered as a first-order surface phenomena):
dm3 = −ke AE (t) dt m3 (t = 0) = m3, in
(7)
At each time step the erosion radius, RE , and thus the erosion surface, AE , are known, allowing for the calculation of the polymer mass, m3 (t), by solving Eq. (7), once a proper value for erosion constant, kE , is given.
2.3.3.
1 ω1 ω2 1 − ω1 − ω2 = + + tot 1 2 3
(8)
where 1 , 2 and 3 are the water, the drug and the polymer densities, respectively. Once they are known, the drug and the polymer masses are known, the current volume of the matrix can be calculated: mtot (t) tot
(9)
Finally, the required value for the erosion radius was obtained by the volume value.
2.3.4.
Geometrical issues
In Eq. (5), (6) and (9) the matrix volume, V, is involved; similarly, on the other hand, in Eq. (7) the external (erosion) surface, AE , has to be known. Recalling the meaning of the geometrical factor f, which is 0 for slabs, 1 for infinite cylinders, 2 for spheres and for finite cylinder, one can write: AE (t) =
f CA RE (t)
(10)
Results and discussions
3.1.
Preliminary validation
The code’s core was validated elsewhere (Barba and Lamberti, 2003) by comparison with known analytical solutions. In this work but here the code was largely modified, thus a new step of preliminary validation was required. Therefore, the code was used to simulate three cases for which the analytical solutions are available: the release from (i) a slab of constant thickness, (ii) an “infinite” cylinder (a cylinder with a large ratio height/radius) of constant radius, (iii) a sphere of constant radius (no moving boundaries means constant density, thus the matrix does not “swell” because of water up-take and no erosion, and the polymer does not dissolve in the medium). For these three cases the analytical solutions are known, under the further restriction of constant diffusivities, and are reported in Table 2. The solutions were obtained by solving Eq. (2) with initial and boundary conditions given by Eq. (3). For each case, Table 2 reports the normalized mass fraction profile of the ith component, ϕi (t, ), as a function of time, t, and of the generic position variable, (which is the axial direction in a slab, and the radial direction in an “infinite” cylinder and in a sphere). Table 2 also reports the cumulative drug release of the ith component, DRi (t). The normalized profile is defined as ϕi (t, ) =
ωi,in − ω(t, ) ωi,in − ωi,eq
(12)
And the cumulative release is defined as the ratio of the ith component mass released at the time t over its initial value:
and f +1
V(t) = CV RE
(t)
(11)
The use of coefficients CA and CV allows to reproduce the desired shape of the matrix, leading to the code generalization. Their values are summarized in Table 1. Thus, during the solution procedure, at each time step the code evaluates: the drug and water masses (Eq. (5)) (and by this way it accounts for the swelling), the polymer mass (Eq. (7)) (and by this way it accounts for the erosion), the overall density (Eq. (8)), the matrix volume (Eq. (9)), and at last the “erosion radius”, i.e. the actual size of the matrix (Eq. (11)). The erosion radius change (increase due to swelling and/or decrease due to the erosion) requires of course an update of the spatial grid, which is performed for each time step.
2.3.5.
3.
Erosion radius calculation
The average density of the partially hydrated matrix can be calculated from the simplest mixing rule which can be written for the specific volume:
V(t) =
equation (2), and the Euler finite difference scheme for the erosion equation (7), the integral equations ((5) and (6)) being approximated by summations. The code was developed using the Mathcad© software (Mathsoft Engineering & Education, Inc.). To ensure numerical solution stability, each simulation was carried out several times by doubling time and space step numbers, until the obtained solution become independent from numbers of these steps.
Numerical issues
The model equations were solved numerically, adopting the Crank & Nicolson finite difference scheme for the diffusion
DRi (t) = 1 −
V(t)
ωi (t, ) dV
ωi (0, ) dV V0
(13)
The code was validated by comparison of its calculations with analytical solutions, obtained assigning to the model parameters reasonable values, based on literature or physical data. Here the following values were used: ω1,in = 0.00 (initial water content in the matrix), ω2,in = 0.60 (initial drug content in the matrix), ω1,eq = 0.80 (equilibrium water mass fraction), ω2,eq = 0.00 (equilibrium drug mass fraction), D1 = 5.6 × 10−9 m2 s−1 (water diffusivity within the matrix), D2 = 12 × 10−9 m2 s−1 (drug diffusivity within the matrix) and RE = 6 mm (slab thickness or cylinder radius or sphere radius). Results of the calculations are reported in Fig. 1, in which the analytical solutions and the code predictions perfectly agree. From these comparisons the code can be considered as validated.
363
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Table 1 – Geometrical factors for code generalization. Geometry
f
CA
CV
Slab (square of side L) Infinite cylinder (height H RE ) Sphere Finite cylinder with constant
0 1 2 2
Not of use Not of use Not of use H/2RE
1 1 1 (22 + 1)/32
2L2 2H 4 2(1 + 2)
2L2 H 4/3 2
3.2.
Modeling for a finite cylinder
der height) and for > 3 the cylinder can be taken as “infinite” (the height being considerably larger than the radius). The values of the correction factor s which causes the numerical solutions to be in agreement with the analytical ones are reported as symbols in Fig. 3 (left). These data were also fitted by a suitable relationship, i.e.
Following the reasonings proposed by Coviello et al. (2003), as summarized in Appendix A, which lead to Eq. (2), the code should be able to describe the behavior of a “finite” cylinder, i.e. a cylinder with comparable height and radius. The diffusion equation in a cylinder (which would arise from Eq. (A2)) is different from the diffusion equation in a sphere (see Eq. (A4), the basis for Eq. (2)). This means that, to make the two solutions coincident (the solution of the 2D diffusion problem into a cylinder of radius RE and the solution of the 1D diffusion problem into a sphere of radius ()RE ), the code needs for a correction coefficient. Actually, the simulations of a finite cylinder with = 2 (no swelling–no erosion – constant diffusivities), obtained by the analytical solution and by the numerical code (Eq. (2)), do not give the same results (see Fig. 2, left). To make the code able to reproduce the analytical result, the right hand side of Eq. (2) has to be multiplied by a factor which is roughly s = 0.17. By using this value, the numerical predictions match the analytical ones (see Fig. 2, right). Then, the same procedure can be followed for several values of the ratio between height and diameter, = height/diameter = semi-height/radius. There is no need to investigate very low and very large values of , since for < 0.2 the cylinder is practically a slab (the thickness being the cylin-
2 A 2 e−(ln(/c )) /2ω s() = s0 + √ 2 w
(14)
in which s0 = 4.32 × 10−3 , lc = 1.7745, w = 0.8769 and A = 0.7652. It is worth to be noted that Eq. (14) does not have any physical meanings, it only assures that the full model, by using f = 2 (i.e. for a sphere) and given by Eq. (A4), can be used in the description of a finite cylinder. In fact, the correction factor given by Eq. (14) allows the code to reproduce the behavior of a finite cylinder (with 0.2 < < 3). Cylinders with values of lower than 0.2 can be simulated as slabs, and cylinders with values of higher than 3 can be simulated as infinite cylinders. Fig. 3 (the right) also reports the simulation obtained by using the numerical code, with the correction factor s given by Eq. (14), compared with the analytical solution, in term of drug release with time. Summarizing, the code given by Eqs. (2)–(11), including the correction coefficient given by Eq. (14) (which holds only for “finite” cylinders), is able to reproduce the diffusion of both
Table 2 – Summary of the equations which hold for constant diffusivity – no moving boundary, and for several geometries. Slab (square of side L and thickness 2RE ) (§ 4.3.2, p. 47, Crank, 1975)
4
ϕi,slab (t, ) = 1 −
∞
∞
ϕi,cylinder (t, ) = 1 −
2 RE
(2n+1)2 2 Di t 4 R2 E
8 (2n+1)2 2
n=0
∞
exp −
J0 (˛n ) ˛n J1 (RE ˛n )
∞
4 R2 ˛2 E n
cos
(2n+1)2 2 Di t 4 R2 E
(2n+1) 2
RE
exp(−Di ˛2n t) where ˛n ’s are the positive roots of J0 (RE ˛n ) = 0
⎧ n=1 ⎪ ⎪ ⎪ ⎨ =0 ϕi,sphere (t, ) =
exp(−Di ˛2n t)
1+2
∞
⎪ ⎪ ⎪ ⎩>0
DRi,sphere (t) = 1 −
6 2
n
2 2
(−1) exp −n
2RE 1+
∞ n (−1)
n n=1
1 n2
exp
−n2 2
n=1
ϕi,finite
cylinder (t, )
DRi,finite
cylinder (t)
Di t R2 E
= ϕi,slab(t, )ϕi,cylinder (t, ) RE
=1−
n=1
∞
Finite cylinder (height H and radius RE , = H/2RE )
n=1
DRi,cylinder (t) = 1 −
Sphere (of radius RE ) (§ 6.3.1, p. 90, Crank, 1975)
exp −
n=0
DRi,slab (t) = 1 −
Infinite cylinder (height H RE ) (§ 5.3.1, p. 73, Crank, 1975)
(−1)n 2n+1
0RE 0
f dRE E f ωi,in ()CV R dRE E ωi (t,)CV R
Di t RE2 2 2
exp −n
Di t RE2
sin
n
RE
364
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Fig. 1 – (Left) Drug and water mass fraction profiles within the matrix (after a given time from immersion). (Right) Evolution with time of the drug release. Symbols: results of numerical code; lines: results of analytical solutions as reported in Table 2. From the top: results for a slab, an “infinite” cylinder and a sphere.
water and drug within the solid pharmaceutical form (“the matrix”), for several shapes: slabs, infinite cylinders, sphere and finite cylinders. The code’s ability was tested by comparison with the analytical solutions, which are available under some restrictions (no swelling–no erosion – constant diffusivities).
3.3.
Case history from literature
The last, most important, challenge for the code is to test its ability in predicting a real case. Here we take a real case from Siepmann et al. (1999b). Cylindrical matrices made of HPMC (a swellable hydrogel commonly used in controlled release
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
365
Fig. 2 – Analytical and numerical drug and water profiles for the diffusion within a finite cylinder with = 2. The analytical solution was obtained as reported in Table 2; the numerical solution was obtained by solving Eq. (2). (Left) The numerical solutions came from Eq. (2) without any correction. (Right) The numerical solutions came from Eq. (2) using a correction factor s().
systems) and of propanolol HCl (a non-selective beta blocker, mainly used in the treatment of hypertension) were prepared by direct compression of powders mixture (5% (w/w) drug, 28 mg of HPMC and 1.5 mg of propanolol HCl). The matrices were 2.5 mm in radius and 0.7 mm in semi-height. Thus the ratio was 0.28. The matrices were immersed in phosphate buffer (pH 7.4) following the USP rotating paddle method. The drug release was followed by periodical assaying of the release medium (spectroscopically, at wavelength of 290 nm). Experimental results were reported as symbols in Fig. 4 (left). In the same graph, the prediction by a 2D model developed in Siepmann et al. (1999b) was also reported (dashed line). The model, proposed by Siepmann, Peppas and co-worker (Siepmann et al., 1999b), is based on the same assumptions
used in the present work, the main difference being the ability of the code here presented to describe the behavior of different shaped matrices (it has to be noted that the code is mathematically simpler, being 1D). To obtain the simulation in Fig. 4, the parameters to be used are: ω1,eq = 0.802 (in Eqs. (3) and (4)), ˇ1 = 2.5, ˇ2 = 9.5, D∗1 = 56 × 10−11 , D∗2 = 6.3 × 10−11 m2 /s (in Eq. (4)) and ke = 5.5 × 10−6 kg/(m2 s) (in Eq. (7)). The parameters have been taken from Siepmann et al. (1999b), no further optimization were performed in this work. In Fig. 4 the predictions of the model presented in this work (with the same values of the parameters) were reported as a solid line. The predictions of the two codes are practically coincident, confirming the ability of the proposed code
Fig. 3 – (Left) The values of the correction factor s(), as obtained by comparison between numerical and analytical solutions, together with a fitting curve. (Right) the cumulative drug release, in the case of = 2, as predicted by analytical solution (Eq. (13)) and by numerical solution (Eq. (2)), using the correction factor s().
366
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Fig. 4 – (Left) The drug release in phosphate buffer (pH 7.4) from a cylindrical matrix of hydroxypropylmethylcellulose (28 mg)/propanolol HCl (1.5 mg), with = 0.28. Experimental data and 2D model predictions came from Siepmann et al. (1999b); the 1D model predictions was obtained in this work. (Right) The time evolutions of water, drug and polymer masses, as predicted by the code.
to reproduce both the experimental data and the predictions of a more complex and less versatile model, only applicable to cylindrical matrices. Fig. 4 reports also the predictions in term of mass evolutions for drug, water and polymer. These evolutions were simulated by the code, and the expected behavior can be easily recognized: the drug and the polymer masses show a monotonically decreasing behavior (the drug diffuses toward the dissolution medium, the polymer is subjected to erosion). The water content firstly increases (due to the up-take from the external medium), and then decreases (due to the erosion which causes the external layers of gel to be destroyed). To be sure that the predicted phenomena really takes place, and when they take place, a comparison with experimental data is required, i.e. the availability of experimental data to be compared to above predictions would be a very powerful piece of information. Currently, our research group is working to define the experimental procedures necessary to obtain these data.
Fig. A1 – The scheme for differential volume identification in finite cylinder analysis.
4.
Conclusions
Controlled release of drugs can be obtained by using matrices made of swellable hydrogels. The availability of a mathematical tool, able to predict the drug release kinetics from matrices of several shapes, dimensions and made of several materials, is of great aid in design and manufacture of new pharmaceutical systems. In this work, a physical model was formulated, to describe all the phenomena commonly observed during hydration of hydrogels made matrices. Water up-take, which causes the hydrogels swelling and the (water’s and drug’s) diffusivities increase, polymer erosion and drug release were mathematically described. The equations were implemented in a numerical code, based on a scheme previously proposed. The code was built on the basis of a generalization (Barba and Lamberti, 2003) which made it able to predict the behavior of matrix with several simple shapes: an infinite slab (thickness negligible with respect to the other dimensions), an infinite cylinder (radius negligible with respect to the height) and a sphere. The first step, was to validate the code by comparing the results with some analytical solutions known for simple cases (constant diffusivity, no swelling–no erosion). Then, on the basis of reasonings suggested in Coviello et al. (2003), the code was made able to describe also the behavior of a finite cylinder, with constant height/radius ratio. The comparison of the modified code with analytical solutions for this case allowed to add a correction factor to the code itself. Finally, the full code was used to simulate the behavior of a real system, reported in Siepmann et al. (1999b). The simplified (1D+time) code proposed here was able to reproduce the experimental data with the same accuracy of the more complex model (2D+time) available in literature. Summarizing, in this work a model suitable to describe the drug release from different shaped matrices made of swellable hydrogels, was formulated and implemented into a numerical
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
code. The code was validated by comparison with analytical solutions, and it was found able to reproduce a real case from literature. Properly designed experiments are highly desirable to further test the code ability to describe real systems. These experiments have to be designed to gather detailed information on the phenomena involved in the drug release (water up-take, matrix swelling, polymer erosion, drug diffusion). Such experiments are currently under development in our laboratory (see for example Chirico et al., 2007b). The code proposed here could be of use both as an aid in quantification of the phenomena observed in R&D studies, and useful in design and manufacture of new pharmaceutical systems.
Appendix A The basic idea (Coviello et al., 2003) to get a simplified description (2D → 1D) of a finite cylinder is to consider a cylindrical shell with a semi-height/radius ratio = z/r equal to its initial value, = Z0 /R0 , and to take this shell as the computational element (Fig. A1). The differential volume of the computational element is 2
V = [(r + r) (z + z) − r2 z] = 6 r2 r
(A1)
The mass balance within the computational volume can be written on the basis of constant concentration in the shell as
V
∂ω = (2 r2znr ) − (2 r2znr ) r r+ r ∂t
+ (2 r2 nz ) − (2 r2 nz ) z
z+ z
(A2)
The axial flux can be related to the radial flux as follows: nz = −D
∂ω D ∂ω 1 = − ∂r = nr ∂z z r
(A3)
Insertion of (A3) in (A2), followed by dividing by V, recalling that = z/r and taking the limit 0, gives ∂ω ∂t = 2 for r1 → 4 ∂ 2 2 ∂ 2 1 ∂ 2 − 6 r2 ∂r (r nr ) − 6 2 r2 ∂r (r nr ) = − 3 + 32 r2 ∂r (r nr )or
∂ω =− ∂t
22 + 1 32
1 ∂ 2 1 ∂ (r nr ) = −() 2 (r2 nr ) r2 ∂r r ∂r
(A4)
Eq. (A4) is similar to the (1D) mass balance within a sphere, but for the correction coefficient (), for which Eq. (A4) is the definition.
references
Barba, A.A., Lamberti, G., 2003. Preliminary validation of a numerical code for heat transfer simulations. Heat Mass Transfer 39, 429–433. Bettini, R., Catellani, P.L., Santi, P., Massimo, G., Peppas, N.A., Colombo, P., 2001. Translocation of drug particles in HPMC matrix gel layer: effect of drug solubility and influence on release rate. J. Control. Rel. 70, 383–391.
367
Chirico, S., Lamberti, G., Nunziata, V., Titomanlio, G., 2007a. Analysis and modeling of radial water up-take in pure HPMC matrices. In: Proceedings of the European Congress of Chemical Engineering (ECCE-6), Copenhagen, September. Chirico, S., Dalmoro, A., Lamberti, G., Russo, G., Titomanlio, G., 2007b. Analysis and modeling of swelling and erosion behavior for pure HPMC matrix. J. Control. Rel. 122, 181–188. Coviello, T., Grassi, M., Lapasin, R., Marino, A., Alhaique, F., 2003. Scleroglucan/borax: characterization of a novel hydrogel system suitable for drug delivery. Biomaterials 24, 2789–2798. Crank, J., 1975. The Mathematic of Diffusion. Clarendon Press, Oxford. Frenning, G., Brohede, U., Strømme, M., 2005. Finite element analysis of the release of slowly dissolving drugs from cylindrical matrix systems. J. Control. Rel. 107, 320–329. Fu, J.C., Hagemeir, C., Moyer, D.L., Ng, E.W., 1976. A unified mathematical model for diffusion from drug–polymer composite matrices. J. Biomed. Mater. Res. 10, 743–758. Grassi, M., Grassi, G., 2005. Mathematical modelling and controlled drug delivery: matrix systems. Curr. Drug Deliv. 2, 97–116. Grassi, M., Colombo, I., Lapasin, R., 2000. Drug release from an ensemble of swellable crosslinked polymer particles. J. Control. Rel. 68, 97–113. Grassi, M., Zema, L., Sangalli, M.E., Maroni, A., Giordano, F., Gazzaniga, A., 2004. Modeling of drug release from partially coated matrices made of a high viscosity HPMC. Int. J. Pharm. 276, 107–114. Higuchi, T., 1961. Rate of release of medicaments from ointment bases containing drugs in suspensions. J. Pharm. Sci. 50, 874–875. Kill, S., Dam-Johansen, D., 2003. Controlled drug delivery from swellable hydroxypropylmethylcellulose matrices: model-based analysis of observed radial front movements. J. Control. Rel. 90, 1–21. Korsmeyer, R.W., Lustig, S.R., Peppas, N.A., 1986a. Solute and penetrant diffusion in swellable polymers. I. Mathematical modeling. J. Polym. Sci. B: Polym. Phys. 24 (2), 395–408. Korsmeyer, R.W., von Meerwall, E., Peppas, N.A., 1986b. Solute and penetrant diffusion in swellable polymers. II. Verification of theoretical models. J. Polym. Sci. Polym. Phys. Ed. 24, 409–434. Narasimhan, B., Peppas, N.A., 1996a. Disentanglement and reptation during dissolution of rubbery polymers. J. Polym. Sci. Polym. Phys. 34, 947–961. Narasimhan, B., Peppas, N.A., 1996b. On the importance of chain reptation during dissolution of rubbery polymers. Macromolecules 29, 3283–3291. Peppas, N.A., Sahlin, J.J., 1989. A simple equation for the description of solute release. III. Coupling of diffusion and relaxation. Int. J. Pharm. 57, 169–172. Siepmann, J., Peppas, N.A., 2000. Hydrophilic matrices for controlled drug delivery: an improved mathematical model to predict the resulting drug release kinetics (the “Sequential Layer” model). Pharm. Res. 17, 1290–1298. Siepmann, J., Peppas, N.A., 2001. Modeling of drug release from delivery systems based on hydroxypropyl methylcellulose (HPMC). Adv. Drug Deliv. Rev. 48, 139–157. Siepmann, J., Poudal, K., Sriwongjanya, M., Peppas, N.A., Bodmeier, R., 1999a. A new model describing the swelling and drug release kinetics from hydroxypropyl methylcellulose matrices. J. Pharm. Sci. 88, 65–72. Siepmann, J., Kranz, H., Bodmeier, R., Peppas, N.A., 1999b. HPMC-matrices for controlled drug delivery: a new model combining diffusion, swelling, and dissolution mechanism and predicting the release kinetics. Pharm. Res. 16, 1748–1756.
368
e u r o p e a n j o u r n a l o f p h a r m a c e u t i c a l s c i e n c e s 3 6 ( 2 0 0 9 ) 359–368
Siepmann, J., Kranz, H., Peppas, N.A., Bodmeier, R., 2000. Calculation of the required size and shape of hydroxypropyl methylcellulose matrices to achieve desired drug release profiles. Int. J. Sci. 201, 151–164. Wu, X.Y., Zhou, Y., 1998. Finite element analysis of diffusional drug release from complex matrix systems. II. Factors influencing release kinetics. J. Control. Rel. 51, 57–71.
Wu, X.Y., Zhou, Y., 1999. Studies of diffusional release of a dispersed solute from polymeric matrixes by finite element method. J. Pharm. Sci. 88, 1050–1057. Zhou, Y., Wu, X.Y., 1997. Finite element analysis of diffusional drug release from complex matrix systems. I. Complex geometries and composite structures. J. Control. Rel. 49, 277–288.