Predicting the Drug Release Kinetics of Matrix Tablets

In this paper we develop two mathematical models to predict the release kinetics of a water soluble drug from a polymer/excipient matrix tablet. The first of our models consists of a random walk on a weighted graph, where the vertices of the graph represent particles of drug, excipient and polymer, respectively. The graph itself is the contact graph of a multidisperse random sphere packing. The second model describes the dissolution and the subsequent diffusion of the active drug out of a porous matrix using a system of partial differential equations. The predictions of both models show good qualitative agreement with experimental release curves. The models will provide tools for designing better controlled release devices.


1.
Introduction. The first Workshop on the Application of Mathematics to Problems in Biomedicine took place from December 17-19, 2007 at the University of Otago in Dunedin, New Zealand. During that workshop, our group of nine, comprising pharmaceutical scientists and mathematicians worked on the problem of predicting the drug release kinetics from a sustained release matrix tablet. The present paper grew out of this event.
Sustained release (also called extended release) tablets are a common dosage form. A sustained release (SR) tablet is typically designed to release drug over such as magnesium stearate might be used at a 1% level, whereas the level of a rate-controlling polymer would be much higher (10-50%) but preferably the level would be kept as low as possible for cost reasons.
When such tablets are placed in water, or in the intestinal fluid of the patient, the coherent polymer matrix sustains release of the drug by keeping the tablet largely intact and by providing a tortuous network through which water penetrates and dissolved drug and soluble excipient molecules diffuse out. That is, the postulated release mechanism involves penetration of fluid, dissolution of the drug and soluble excipient in the fluid, and outward diffusion of molecules of dissolved drug due to the concentration difference between the solution in the tablet and the intestinal fluid. Once released from the tablet, the drug is rapidly absorbed through the patient's intestine into the blood stream. Ideally, the delivery rate of the drug from the tablet should be such as to maintain the preferred blood level for an extended time. The left panel of figure 1 shows the structure of the polymer matrix of a tablet after dissolution of the solvable excipient and the drug.
In this paper we propose two mathematical models that predict the release kinetics of matrix tablets, taking into account parameters such as the composition of the powder mixture, powder particle sizes and applied curing temperature. The first model models the diffusion process as a random walk through the channels in the tablet that are formed in the compaction and heating process. In this model, the particle distribution within the tablet is represented by a graph embedded in R 3 whose vertices are generated from a random sphere packing and in which two vertices are joined by an edge if the corresponding spheres are in close proximity. On this contact graph we perform random walks that start at the location of each drug particle. These random walks mimic the diffusion of the drug molecules from their initial position to the edge of the tablet. By adjusting the probabilities of moving to vertices representing different types of particles, the composition of the tablet affects the length distribution of these exit paths in a realistic way. The right panel of figure 1 shows a sample path through a random packing of disks. Field emission scanning electron microscope (SEM) image of a matrix tablet after dissolution of drug and excipient at 400× magnification. Right panel: A sample path through a schematic two-dimensional tablet from an inner particle to the edge. Small solid circles represent polymer particles. Notice that there are multiple paths from each initial particle location.
Low levels of drug in some polymer systems can lead to trapping and incomplete release. Percolation theory has been applied to understand these processes [3-6, 10, 17]. Monte Carlo simulations of drug release from a binary powder mixtures with inert (water insoluble) excipient have been performed by Villalobos et al. [30,31]. The drug and excipient particles were placed on the sites of a cubic lattice and the drug particles performed random walks on the accessible empty sites under a volume exclusion constraint. The authors reported different exponents for the release profiles as function of time, with √ t-like behavior for large drug fractions. In our work we take into account drug, excipient and polymer, where the polymer plays the role of the inert matrix. More importantly, we also allow for differences in particle sizes and irregular placements of the particles in space.
The second model investigates a different approach to predicting release kinetics of tablets using partial differential equations, where time and space are treated as continuous variables. Previous work using this approach includes [16,24]. Suitable versions of the diffusion equation (often in cylindrical coordinates, resembling the natural shape of a tablet) are formulated for the diffusion of dissolved drug and (e.g. in [24]) water penetrating the tablet. These models can take into account the dissolution kinetics of the drug in different types of buffer. Siepmann and Peppas [24] consider hydrophilic polymer matrices that swell and then dissolve, unlike the polymers used in the tablets that we are studying. While these results may not be directly comparable to our results for inert matrices, they do show that an increasing drug load results in faster release kinetics [24, Fig. 4].
2. Model 1: The diffusion of drug molecules is a random walk on a weighted graph. This first model is constructed in three distinct steps. In the first step, the positions of the particles in the tablet are taken from a random sphere packing, where each particle is a sphere of a fixed radius. In Step 2, this random sphere packing is represented as a graph embedded in R 3 , whose vertices are the particles, with edges between particles that are close to each other in space. The compression of the tablet is modeled as a deformation of this graph, and the heating of the tablet is modeled by removing edges. Finally, in Step 3, the diffusion of the drug particles to the exterior of the tablet is modeled as a random walk on the graph generated in Step 2. Drug release profiles are then generated as the distribution of exit times of these particles, using the following argument.
Since we model the diffusion process as a random walk, the distance that a particular drug molecule must traverse from its position inside the tablet to the exterior is a random variable, L. Let f L (s) denote the density of these path lengths. Assume that a dissolved drug molecule performs Brownian motion with diffusion constant D along these paths. If the particle has to travel along a path of length s to get to the exterior of the tablet, then the travel time T is a random variable with density function (see [11,VI.2 In other words, f T |s is the density of T conditioned on the path-length being equal to s. Thus, we can integrate this conditional density against the density of path lengths to get the density of escape times of the particles from the tablet It's cumulative distribution function gives the release profile of the matrix tablet. In Step 3, we estimate F T (t) by performing many random walks on the graph created in Step 2.
We now describe each of these three steps in more detail. Results from preliminary implementations of this model are given in Section 4. In this work we are not interested in supplying the optimal algorithm for completing each step in the modeling process. Rather, we are trying to validate the idea and, if the results appear sound, we can then optimize the algorithms. We assume that drug, excipient and polymer particles can be modeled as perfect hard (i.e. non-deformable) spheres, where all particles of each type have the same radius, but the three radii of the three different types of particles may differ from one another. The assumption that all particles of one species have the same radius is a simplification. In practice, their sizes are distributed around a mean with a given variance. Relaxing this assumption, as well as including particles of other, non-spherical shapes are topics of future research.
Step 1: We assume that prior to compaction the powder mixture is a dense multidisperse random sphere packing. Dense packings of disks or spheres have attracted the interest of mathematicians, physicists and engineers for a long time, see e.g. [9,14,18,27,29] and the references therein. It turns out to be extremely difficult to give a precise mathematical meaning to the concept of a "dense random sphere packing", as the desired properties of high density and high degree of randomness conflict with each other [29]. The more dense a packing is, the more ordered it tends to be, culminating in the highly ordered lattice packings such as the tridiagonal packing of disks in R 2 and the face-centered cubic packing of balls in R 3 . (The optimality of the latter packing was conjectured by Johannes Kepler in 1611 and proved by Thomas Hales in 1998.) While this question is certainly interesting from a mathematical point of view, we shall not pursue it further. Rather, we will work with the outcome of an experimental protocol that produces jammed configurations of hard spheres.
The protocol used in this paper was first suggested by Lubachevsky and Stillinger [18] and later expanded by Knott and coworkers [14]. The paper by Knott et al. [14] discusses the problem of packing spheres of oxidizer and fuel in a solid rocket propellant. Briefly, the Lubachevsky-Stillinger (LS) protocol proceeds as follows. A random initial configuration of sphere centers (x 1 (0), x 2 (0), . . . x N (0)) and a random set of initial velocities (v 1 (0), v 2 (0), . . . v N (0)) are chosen. The spheres are partitioned into M radius classes, each class contains N i spheres (i = 1, . . . , M ) and N = M i=1 N i . The spheres move within a container and collide with each other and with the walls of the container. As time evolves, the radius of class i spheres grows according to r i (t) = a i t. Thus the ratios r i (t)/r j (t) are preserved throughout the process. An event is a collision of two spheres or the collision of a sphere with the wall. At each collision of two spheres the velocities of the colliding spheres are updated in such a way that the spheres move away from each other after the collision. This requires an increase in mechanical energy since the radii are growing at a constant rate. Eventually, as the packing becomes more and more dense, the time between two collisions approaches zero and the process is stopped. The resulting configuration of spheres is (nearly) rigid, in the sense that no sphere can be moved while keeping the centers of all other spheres fixed. For further details we refer to [14,18]. It should be noted that there exist alternative methods to construct random dense sphere packings [13,25,28].  table 2). The packing density is approximately 0.63.
Step 2: Having obtained a multidisperse dense random packing of spheres, this sphere packing is compressed as the powder is compressed in the die [25]. We model the die as a right circular cylinder whose axis is the z-axis and whose lower surface lies and remains in the xy-plane. The top surface of the cylinder which was originally in the plane {z = 1} is displaced by α ≤ 1. A particle originally located at the level z is displaced to the level αz. Using the new coordinates, a graph G is constructed, the proximity graph of the compacted sphere packing. This is a Euclidean graph G = {V, E} whose vertices are the compressed sphere centers. Each vertex carries a label L : V → {D, P, X} that indicates whether the sphere is a drug, polymer or excipient particle, respectively. Two vertices are connected by an edge if the corresponding spheres are sufficiently close to each other. We also introduce a surface indicator s : V → {0, 1} where s(v) = 1 if the sphere centered at v is in contact with the exterior of the tablet.
Next, we form the heated contact graphĜ. If the centers of two spheres are connected by an edge, this edge is removed with a probability p that is an increasing function of the heating temperature T h and the duration of the heating process t h . This reflects the fact that polymer particles will amalgamate and block the access of solvent to certain particles. In the extreme case we observe drug or excipient particles that are completely surrounded by a coat of molten polymer particles. This results in a certain amount of drug that is not released within a physiologically relevant time span. Such behavior has been observed in experimental release curves, see figure 7.
Step 3: The estimation of the distribution of path lengths, f L (t), and the modeling of the diffusion process on these graphs is streamlined in this implementation. Theoretically, each molecule undergoes one dimensional Brownian motion along a path from its initial location until it reaches the edge of the tablet.The onedimensional diffusion can be written as a limit of one-dimensional random walks, where the step size, ∆x, and the time between steps, ∆t, go to zero in such a what that, in the limit, the variances do not tend to either zero or infinity. This requirement is met by requiring that where d is the diffusion constant, see [11, XIV.6, p. 325]. As a coarse approximation to the one-dimensional diffusion, therefore, we allow the molecules to perform random walks along the vertices of the heated contact graphĜ, starting at every vertex of type D. These random walks approximate the diffusion process that each drug molecule undergoes in order to exit the tablet, and therefore give us an approximation of the distribution of exit times, f T . In practice, the physical length of each step in the random walk is related to the diameters of the spheres representing the original particles, hence the actual exit time is proportional to the sum of the lengths. In the simulations used to generate the distribution of exit times shown in figures 4 and 5, the radii of all the original particles were approximately the same, hence, in this case, the number of steps in the random walk was used as a surrogate for the time to exit the tablet. Further work must be done to investigate the scaling properties of these random walks, as well as the effect of including particles of very different sizes. However, our preliminary investigations indicate that the distribution of the number of steps in the walks is at least qualitatively the same as the distribution of exit times. The effect of the polymer, in particular the effect of the fusing of polymer particles in the heating process, is modeled by adjusting the probability of traversing a particular edge in the contact graph. Different rules can be implemented and tested. The simplest is to assign a conductivity c ij to the edge that connects vertices v i and v j . The conductivity is large for edges that connect vertices of type D or X, and small for an edge to or from a particle of type P . Let N (i) be the set of vertices neighboring vertex i, v i . Note that i ∈ N (i), so a molecule can remain at its initial position. If the walker, representing a drug molecule, is situated at vertex v i then the probability of moving to the adjacent vertex v j is given by The walk is stopped if either an exterior vertex or a prescribed maximum number of steps N max is reached. Notice that the metric onĜ is the Euclidean distance between its vertices inherited from its embedding in R 3 , so that the distance a molecule travels on the graph corresponds to the actual distance it must travel in space. If a random walk does not reach one of the exterior vertices within N max steps, then the molecule is considered trapped inside the tablet. This will be used to explain partial release of drug at higher polymer concentrations.
3. Model 2: The continuous model. In the second model we regard time and space as continuous quantities and we set up a system of partial differential equations for the contents of dissolved and undissolved excipient and drug in the tablet. Our starting point is the assumption that fluid quickly permeates the tablet [16] and that dissolution is limited by the saturation of the fluid with excipient and drug. This causes a delay in the initial release of drug molecules until porosity increases through diffusion at the boundary. This delay can be observed as a change in concavity near t = 0 in the experimental release profiles, see figure 7.
Let Ω be the spatial domain of the tablet. We introduce cylindrical coordinates (r, θ, z) such that and assume that our tablet dissolves symmetrically, i.e. concentrations do not depend on the angular variable θ [24]. Let u 1 be the concentration of dissolved excipient in the solvent. Let u 2 be the concentration of undissolved excipient in the solid remainder of the tablet. Likewise, we denote by v 1 and v 2 the concentration of solved drug in the solvent and the content of undissolved drug in the tablet, respectively. All these quantities have dimension M L −3 . By κ we denote the porosity of the tablet, κ is a dimensionless value between 0 and 1 and will increase as more and more excipient and drug are dissolved in the solvent. Then κu 1 is the concentration of solved excipient in the tablet. We model the diffusion classically by assuming that the flux of solved excipient across an interface is proportional to the concentration gradient times the area of the interface. The area of an interface is also proportional to κ, i.e. the flux is given by where D u is the diffusion constant of the dissolved excipient (and D v is the diffusion constant of the dissolved drug). The conservation of mass equation then yields where g(u 1 , u 2 ) is the rate of concentration increase from dissolving excipient. The higher the porosity, the more contact there is between the solid excipient and the fluid in the pore space and hence the rate of dissolution increases with porosity. We assume this to be linear for simplicity but plan to explore more complex relationships in future work. A similar modeling assumption has been made by Lemaire et al. [16], see in particular equation (1) in that paper. The higher the concentration of dissolved excipient, the slower the rate of dissolution with no dissolution taking place at the saturation point C u max of the fluid. Hence we assume that PREDICTING DRUG RELEASE KINETICS 9 and we obtain the following system of evolution equations The rates of dissolution of excipient and drug are denoted by α u and α v , respectively. The porosity κ(u 2 , v 2 ) depends on the concentration of undissolved excipient and drug in the tablet. We assume that where κ 0 is the initial porosity and κ end is the porosity of the tablet once all the excipient and drug are dissolved. The initial concentrations of undissolved excipient and drug are denoted by u 0 2 and v 0 2 . We assume that initially no excipient and drug are dissolved and that the solid excipient and drug are uniformly distributed, , for positive constants u 0 2 and v 0 2 . If drug and excipient are completely undissolved then the initial porosity κ 0 will be about 2 %, that is the porosity after compaction. If drug and excipient are completely dissolved then the porosity will be about κ end = 60 %. Equations (3) are completed by homogeneous Dirichlet boundary conditions for u 1 and v 1 u 1 = 0, v 1 = 0 on ∂Ω, as we assume that any dissolved excipient or drug outside the tablet is immediately carried away.
We add the first two equations of system (3), integrate over Ω and apply the divergence theorem. We obtain, after exchanging the order of integration and differentiation where n denotes the outward normal and σ the surface measure. This is the rate of change of the excipient load inside the tablet and defines the flux of excipient −J u (t) across the boundary ∂Ω of the tablet. (Alternatively we could integrate the normal component of the flux in (1) over the boundary). A similar calculation gives the flux of drug The cumulative amount of drug released is then given by The solutions u 1 , u 2 , v 1 and v 2 remain positive and lim t→∞ ||u 1 (t)|| ∞ + ||u 2 (t)|| ∞ + ||v 1 (t)|| ∞ + ||v 2 (t) ∞ = 0.

It follows then that lim
that is, the initial drug load will be completely released. We will return to this point in section 5.

4.
Results. In this section we describe the implementation of our models outlined in sections 2 and 3. All procedures have been implemented in matlab and programs will be available upon request from the corresponding author.
4.1. The Lubachevsky-Stillinger protocol. We work with a version of the Lubachevsky-Stillinger protocol [14,18] that allows for classes of spheres with different radii. In the first step we convert the mass fractions of drug f D , polymer f P and excipient f X into particle fractions, using the known particle radii r D , r P and r X and the densities D , P and X of the pure powders. For i ∈ {D, P, X} let m i = r 3 i i and setñ these are the particle fractions of drug, polymer and excipient in the powder mixture. The particle sizes, densities and mass fractions for different mixtures of the drug indomethacin (a commonly used anti-inflammatory drug), the polymer Eudragit RLPO and the excipient mannitol are given in table 1. The corresponding particle fractions are given in table 2.
A difficulty in the implementation of the Lubachevsky-Stillinger protocol is the decision when to stop. In theory, the time between two collisions (either spheresphere or sphere-wall collisions) approaches zero and hence the increase in the packing fraction approaches zero. We stop the execution if the times between successive collisions t k+1 − t k , t k+2 − t k+1 , . . . , t k+n − t k+n−1 < , for a number n and a constant > 0 that are set by the user. The role of n is to avoid premature termination if accidentally a time between two collisions is very short. In our implementation we choose = 10 −6 and n = 50. This results in packing fractions of ≈ 0.54 for packings of spheres with roughly equal radii (figure 2). This value is somewhat below 0.64, the "generally accepted" packing fraction for a random dense monodisperse packing [29]. However, we consider it sufficient in this preliminary work.
At the termination of the Lubachevsky-Stillinger protocol at time t * we obtain a set of positions (x i (t * )) N i=1 and radii (r i (t * )) N i=1 . After possible compression of this sphere packing, we define the graph G by joining vertices x i and x j if where λ ≥ 1 is a constant. Graphs of this kind of graph are also known as proximity graphs on random point sets [19,20]. The heating process removes edges with a probability p and results in the heated contact graphĜ. Figure 3 shows a contact graph of a packing before and after heating. Choosing λ slightly greater than 1 takes into account that the molecules can travel through interstitial spaces between particles that are close, but not touching. The spheres that have at least one surface point close to a boundary of the domain are given the "exterior" label. For example, if for a constant µ ≥ 1, then that sphere is a possible point of exit in the x-direction. Figure 3. The proximity graphs on the sphere centers of a random dense sphere packing. The constants in equations (6) and (7) are chosen λ = µ = 1.03. The left panel shows the contact graph G before the heating process. The edges are removed with a probability p = 0.3 during the heating process, this results in the heated contact graphĜ on the right.

Release curves obtained from random walks on the contact graph.
On the heated contact graphĜ we perform random walks that start from each drug particle and end when an exterior vertex is reached, or the maximum number of steps, N max , is achieved. As explained above, the probability of traversing edge E ij is determined by the type of the terminal vertex, we use here c + = 100 and c − = 1. Notice that the current position is also a possible terminal vertex, that is, the random walker can remain in a fixed place for (perhaps long) periods of time. For every path we record the number of steps in the path and the length of the edges traversed. If a particle remains at a vertex, we still add "one" to the number of steps, and we add the diameter of the current vertex' particle type to the path length. Thus, even though a molecule may be close to where it started after one time step, the time step is still counted as "time spent before exiting". We present simulated release profiles in Figures 4 and 5. In these figures we vary the polymer fraction and keep all other parameters constant. We observe that as the polymer fraction increases, the drug release is slowed down. Further, we test different edge removal probabilities (p = 0.3 and p = 0.5) that can be interpreted as different curing temperatures. An increase in p decreases the fraction of drug released.

4.3.
Release curves predicted by the continuous model. As a first simplification we disregard the height of the tablet and collapse it to a disk, so that functions are now only dependent on the variable r ∈ [0, 1]. We split the numerical solution procedure into a diffusion step for u 1 and v 1 and a reaction step for all four concentrations. It is assumed that during the diffusion step the porosity κ(u 2 + v 2 ) does not change. We discretize the radial Laplace operator on the uniform grid r k = (k − 1)∆r, k = 1, . . . , N using standard finite differences and we use the Crank-Nicolson scheme at every diffusion step [23,Section 17.3]. The reaction step is calculated using the Euler forward method. Release curves obtained from numerical solution of equation ( (6) and (7) and a probability for edge removal p = 0.1. The edge propensities are given in equation (8). The maximum number of steps is N max = 10 3 . and their release profiles were determined, see [8] for a detailed description. Briefly, powder mixtures were compressed at different pressures and tablets were heated at selected temperatures. The tablets were then placed in a phosphate buffer medium and samples were collected at different time points over a period of 8 h (see figure  7). We observe a good qualitative agreement with our simulated release profiles in figures 4 and 5.
5. Discussion and Conclusion. Several studies have investigated percolation on regular lattices and graphs obtained from random sphere packings. Villalobos and collaborators [30,31] used Monte Carlo simulations on cubic lattices to predict drug release profiles from binary drug/excipient mixtures. In these works the excipient played the role of the inert matrix, which in our case is the polymer. It was reported that if the particle fraction of the inert matrix (the inaccessible sites in the lattice) is below a threshold of ≈ 69%, then the release of the drug is slowed down with increasing matrix fraction, but will still be complete. Only above the critical matrix concentration an entrapment of drug will be observed. The value of 69% is complementary to the estimated site percolation threshold of 0.3116 for the cubic lattice [26]. A similar value was reported by Powell [21,22] for percolation on the contact graph of a monodisperse random sphere packing. In this work we include the heating of the tablet, which melts the polymer. The blocked polymer particles and the removal of certain edges provide a case of what is known as site-bond percolation [7]. In their paper [7], Chang and Odagaki studied site and bond percolation processes where sites and bonds are removed independently from each other and eventually open sites are connected if they share an open bond. They found that for the cubic lattice, the percolation threshold as a function of the probabilities that sites respectively bonds are open, is well described by a hyperbola. While this was not the primary goal of our study, we ended up addressing the problem of determining percolation thresholds for contact graphs of dense random sphere packings. While we have seen a good qualitative agreement between our simulated release profiles and experimental data, further research is needed to provide more reliable and quantitative predictions. In our simulations we find many particles with very short escape paths. However, their number is going to diminish as we consider sphere packings with more spheres. Intuitively, the fraction of spheres "close to the exterior" is decreasing as the number of spheres in the packing increases, a problem that was discussed already in [21]. Larger simulations, perhaps using parallel computing methods, are required. Further work will also include careful calibration of the parameters λ, µ and p to experimental release curves.
Our continuous model captures nicely the change of the release curves from convex to concave. However, as figure 6 and equation (4) suggest, the release will always be complete. No percolation behavior is exhibited by the system (3). In Figure 7. Release of indomethacin (mass fraction 10%) from Eudragit RLPO matrix tablets containing mannitol (90 − 125 µm particle diameter), a plastic excipient, using the USP basket apparatus at 100 rpm in 900 ml phosphate buffer pH = 7.2 (0.2 M ) at 37 • C. The top figure shows a tablet that was compressed at 74 M P a and cured 24 h at 40 • C, the bottom figure shows a tablet that was cured 24 h at 70 • C. Data are representative for three repetitions of the experiment (from [8]). the future we need to take into account the permeability of the porous matrix, a notion that arises in hydrology [2]. In the derivation of the dissolution kinetics (2), we have assumed that the dissolution of drug and excipient is helped by an increase in porosity. This is contrary to the commonly made assumption of a receding boundary in the pharmaceutical literature [12]. There, the drug dissolves in a way that decreases the area of contact between drug and solvent. It poses an interesting inverse problem to distinguish these two concepts, given the experimentally determined release profiles.
Another interesting suggestion of Villalobos et al. is that the release profiles have a Weibull cumulative distribution function, [30,Equation 2].
where the parameters are influenced by the size of the matrix tablet. More precisely, the ratio of the number of exit sites versus the total number of sites, N exit /N should influence the release kinetics, with faster release for smaller values of N exit /N . If we assume that N exit ∝ N