Release mechanism and parameter estimation in drug-eluting stent systems: analytical solutions of drug release and tissue transport

Drug-eluting stents have significantly improved the treatment of coronary artery disease. They offer reduced rates of restenosis compared with their bare-metal predecessors and are the current gold standard in percutaneous coronary interventions. Drug-eluting stents have been approved for use in humans since 2002 and yet, despite the intensive research activity over the past decade, the drug release mechanism(s) and the uptake into the arterial wall are still poorly understood. While stent manufacturers have focussed primarily on empirical methods, several mathematical models have appeared in the literature considering the release problem, the uptake problem and also the coupled problem. However, two significant challenges that remain are in understanding the drug release mechanism(s) and also the determination of the various parameters characterizing the system. These include drug diffusion coefficients and dissolution constants in the stent polymer coating as well as drug diffusion coefficients, binding/uptake rates and the magnitude of the transmural convection in the arterial wall. In this paper we attempt to address these problems. We provide analytical solutions which, when compared with appropriate experiments, may allow the various parameters of the system to be estimated via the inverse problem. The analytical solutions which we provide here for drug release in vitro may thus be used as a tool for providing insights into the mechanism(s) of release.


Introduction
An arterial stent is a small scaffold-like medical device used to increase the size of the lumen when it has become narrowed due to the formation of atherosclerotic plaque. These stents are now routinely coated with a drug to counteract the inflammatory response following insertion into the diseased artery. Compared with their bare metal predecessors, the so-called drug-eluting stents (DESs) offer reduced rates of restenosis and thus represent the current gold standard in percutaneous coronary interventions. While DESs have been around for over a decade, the drug release mechanism(s) and the uptake into the arterial wall are still poorly understood.
The problem of modelling drug release from arterial stents was presented at the fourth UK Mathematics in Medicine Study Group (see Green et al., 2005) and many publications on this topic have appeared in the literature to date. One of the most important aspects in the performance of DESs is the drug release profile. Owing to the expense and time associated with in vivo experiments, stent manufacturers routinely test the release of drug from their stents in an in vitro environment. While this is unlikely to replicate the in vivo situation, where flowing blood, pulsatility, wound healing, proliferation, migration of cells and complex uptake/binding no doubt all play some part, it nonetheless provides the manufacturer with an idea of the shape of the release profile and allows for comparison between different stent designs. Furthermore, it allows the manufacturer to test the quality of the product and the repeatability of the release profile. While stent manufacturers primarily use empirically based methods, when coupled with appropriate modelling, the in vitro experiments may provide insights into the mechanism(s) of release: this is a topic of some controversy in the literature with diffusion, dissolution, polymer swelling and degradation all of which are cited as possible release mechanisms. Several authors have attempted to model the drug release from these devices. For example, referring to polymer-coated DESs which is the focus of this article, Zhao et al. (2012) presented an analytic solution of a cylindrical diffusion model to describe the experimental drug release of everolimus from a Dynalink-E stent while Hossainy and Prabhu (2008) presented a mathematical model for predicting the drug release from a DES coating, based on two discrete modes of transport. Most models assume that Fickian diffusion plays an important role in the release process. Siepmann and Siepmann (2012) in their review on modelling of diffusion-controlled drug delivery provide a series of analytical solutions for drug release from reservoir and monolithic drug delivery systems, some of which may be applied to DESs under certain assumptions. Many of their solutions are early or late time approximations, and in some cases steadystate solutions. A number of authors have focussed on drug release from biodegradable polymers: we do not consider such polymers in this article, but the reader is referred to Formaggia et al. (2010), Rossi et al. (2012) and Fredenberg et al. (2011).
In reality, of course, the stent and the arterial wall are a coupled system and this has also received much attention in the literature, with most of the models necessarily requiring to be solved numerically due to the complexity of the problem. One exception is the model of Pontrelli and de Monte (2010) who consider diffusion-based drug release from the stent coupled with a convection-diffusion equation in the arterial wall which also accounts for drug consumption via a linear reaction. They obtain an analytical solution through separation of variables and use this to show the effect of the transmural velocity, drug metabolism and the amount of drug in the tissue. More recently, Pontrelli et al. (2013a) presented a semi-analytical expression for the drug concentration and mass in each layer of the arterial wall for the case where the drug must dissolve in the polymer before it can diffuse. Models which incorporate more sophisticated binding/reaction terms includes those of McGinty et al. (2010), McGinty et al. (2013), Horner et al. (2010 and Abraham et al. (2013), who all assume an equilibrium reaction. While McGinty et al. (2013) made some analytical progress, the other authors employed numerical techniques to solve their equations. Perhaps the most sophisticated model is that of Tzafriri et al. (2012), who provide a second-order saturable reversible binding model: this model not only requires to be solved numerically, but also includes a number of parameters which are difficult to measure experimentally.
In most of the models described above, many parameters characterizing the system are required to be measured, thus posing a great challenge to experimentalists. These include drug diffusion coefficients and dissolution constants in the stent polymer coating as well as drug diffusion coefficients, binding/uptake rates and the magnitude of the transmural convection in the arterial wall. Values for these parameters are typically taken from the available data in the literature which encompasses a great number of studies using different experimental species (for example, rabbit or pig). Furthermore the experimental techniques are often quite different. For example, many experiments are conducted in vitro when the transmural velocity has necessarily been neglected and in many cases 'lumped diffusion coefficients', which inherently include such effects as convection and binding, are measured. Thus great care must be taken to ensure the accurate estimation of the model parameters. In this paper we attempt to address this problem by providing analytical solutions which, when compared with appropriate experiments, may allow the various parameters of the system to be estimated via the inverse problem. Of course, the inverse approach we propose here to estimate model parameters may also be applied to other drug delivery systems, for example, delivery by drug-filled balloons and transdermal patches (see, for example, Stark et al., 2013;Pontrelli et al., 2013b). This paper is organized as follows. Firstly we present a collection of models to describe diffusionand dissolution-based release of drug from polymer-coated arterial stents in an in vitro environment. The boundary conditions have been chosen to reflect the conditions maintained in the experiments. All of the models are solved analytically, providing user-friendly solutions which can be quickly compared with observed release profiles. It is intended that when used in conjunction with appropriate experimental data, these solutions should help us to clarify the release mechanism(s) and thus aid in the development of DESs. By consideration of the inverse problem, the various parameters characterizing the release may be determined. We then proceed to consider a model of drug transport through the arterial wall. This model too is solved analytically, providing three interesting characterizations of the solution. We demonstrate the equivalence of these solutions and, by comparison with experimental data via an inverse problem, attempt to uncouple and estimate both the tissue drug diffusion coefficient and convection parameter.

Diffusion-based models of drug release
Most of the modelling of drug release from arterial stents in the literature has thus far been concerned with first generation DESs. These stents typically comprise a stainless steel platform with a drug containing polymer coating attached to the stent struts (Tzafriri et al., 2012;Stefanini & Holmes, 2013). The philosophy behind this design was not only to allow the drug to be released gradually so as to avoid toxic levels of drug initially, but also to permit sustained delivery over many weeks. The drugs used (sirolimus and paclitaxel) are lipophilic and are able to inhibit SMC proliferation and migration. Drug release from these stents has been modelled as a diffusion-dominated process (see, for example, McGinty et al., 2010;Pontrelli and de Monte, 2010), with the drug concentration in the polymer C p satisfying a diffusion equation with drug diffusion coefficient D p . Several simplifying assumptions are made. Firstly, it is assumed that the device geometry is that of a thin film of thickness L p with no edge effects so that the modelling may be restricted to one dimension. The diffusion of the drug in the polymer is thus considered to be isotropic and it may be assumed that the diffusion coefficient is independent of time, space and concentration. Furthermore, the initial drug concentration is usually taken to be uniform. For the case of in vitro drug release a zero flux boundary condition is normally assumed at the impermeable stent boundary and either an infinite sink or Robin-type boundary condition at the interface with the release medium. If the release medium volume is sufficiently large and is replenished between sampling times then the infinite sink condition is most appropriate, whereas if the concentration of drug in the release medium is allowed to build up, the Robin-type condition is preferred. The model is as follows: where γ is some constant of proportionality with dimensions ms −1 . Of course (2.4a) can be obtained from (2.4b) by letting γ → ∞ which physically means that the flux across the interface is extremely rapid. For the case of the infinite sink boundary condition, the system (2.1-2.4a) are readily solved using the method of Laplace transforms (or otherwise) to provide (see, e.g. Crank, 1975) (2.5) We can obtain the mass of drug per unit area in the polymer, M (t), by integrating C p (x, t) over L p . From this it is straightforward to show that the cumulative fraction of drug released, where M 0 = L p C 0 has been assumed. For the case where the concentration is allowed to rise in the release medium (i.e. employing (2.4b) rather than (2.4a)), we can again utilize Laplace transforms along with the Residue Theorem to provide the solution where α n are the roots of In this case the cumulative fraction of drug released turns out to be . (2.8) In Fig. 1 we compare the release profiles obtained from Equations (2.6) and (2.8). We assume the typical values of diffusion coefficient and polymer thickness of D p = 10 −16 m 2 s −1 and L p = 10 −5 m, respectively (McGinty et al., 2010). The solution (2.8) is plotted for three values of γ , 10 −10 ms −1 , 10 −11 ms −1 and 10 −12 ms −1 . It is evident that for the case where the drug is allowed to rise in the release medium, the release of the drug can be considerably slowed down. Comparison of release profiles for the cases of release medium being replenished (infinite sink) and drug rising in the release medium (Robin boundary condition). The top profile displays the solution (2.6), while as we move down, the profiles display the solution of (2.8) with γ = 10 −10 ms −1 , γ = 10 −11 ms −1 and γ = 10 −12 ms −1 .

Diffusion-dissolution-based models of drug release
In the preceding section, we considered models of drug release whereby the drug was assumed to be free to diffuse. However, this need not be the case. In fact it is common for the drug to be loaded on the polymer in such a way that the initial concentration exists in two forms: crystalline (bound) and free (unbound). In this case the drug needs to dissolve before it can freely diffuse. This can be achieved by the ingress of fluid into the polymer after placement in the release medium. In this section we describe three models which incorporate aspects of dissolution as well as diffusion. In all of these models we assume an infinite sink boundary condition.

Instantaneous dissolution model
In this case we assume that the initial concentration of drug in the polymer exceeds the solubility limit, say C s (assumed to be constant). Thus drug exists in a crystalline form and is permitted only to diffuse after it has dissolved. Where the initial drug concentration is below the solubility limit, C 0 < C s , the dissolved drug is free to diffuse. In this case the solutions given by (2.5-2.6) describe the concentration profile and cumulative fraction of drug released. However, if C 0 > C s then there is an excess of drug which is immobile and must dissolve before diffusion can commence. This problem was first considered by Higuchi (1961) who derived an expression for the amount of drug absorbed at time t per unit area. He obtained his solution from geometrical considerations and Fick's law. A more mathematical approach has been considered by, among others, Paul & McSpadden (1976). However, these approaches derive a solution that is only valid until the time at which the concentration of drug throughout the polymer falls below the solubility, say t = t . Here we provide the complete solution which is valid for all time.
Suppose that there is a small region immediately adjacent to the release medium boundary where the concentration of drug is less than the solubility, i.e. it is permitted to diffuse. Suppose further that to the left of the region, the drug concentration is above solubility. Let there be a moving boundary between these two regions whose time varying location is given by x = s(t). As drug is released from the polymer the moving boundary tracks back from the (initially fixed) boundary with the release medium to the impermeable stent boundary. We assume that the dissolution is instantaneous so that the concentration on the moving boundary is equal to the solubility. Then the model is coupled with the Stefan-type boundary condition The problem is readily seen to have a similarity structure.
The solution for t < t is (see Vo, 2012or Paul & McSpadden, 1976: Note t is known and is given by The expression for the cumulative fraction of drug released is found to be Now, when all of the drug in the polymer is below solubility, the form of solution changes. At t = t , the concentration of drug in the polymer is clearly given by (3.14) The solution for the duration of drug release is then obtained by solving (2.1), (2.3) and (2.4a) with (3.14) replacing (2.2). Proceeding to solve by separation of variables, the concentration of the drug in the polymer can be shown to be (see Appendix A) where {z} denotes the real part of z. The corresponding expression for the cumulative fraction of drug Recall that t = (L p /θ ) 2 where θ is given by (3.12).

Dissolution kinetics involving a first-order reaction
In the previous case we considered a situation where the dissolution of drug from crystalline into dissolved form is instantaneous. This may be reasonable where the volume of fluid available to dissolve the drug is large. Consider now the case where the dissolution rate is finite, say K 1 . We will now denote by C b the bound drug which is in a crystalline form and C f the free (dissolved) drug which is initially zero. The model in this case is The solution of the system (3.17-3.22) is detailed in Appendix B and is given by The solution for the cumulative mass of drug released is . (3.25)

Dissolution kinetics involving a first-order reversible reaction
In this case we allow for a reversible reaction, i.e. the drug can dissolve from the crystalline form into free form and also from free form into crystalline form. We assume that the forward reaction rate is K 1 and the backward reaction is K 2 . The model in this case is coupled with (3.19-3.22). The solution as detailed in Appendix B is where (3.30) The cumulative fraction of drug released is then given by . (3.31)

Graphical results
In this section we compare the release profiles obtained from the solutions to the diffusion-dissolution models described in Sections 3.1-3.3. For each of the plots we assume typical values of the diffusion coefficient and the polymer thickness, D p = 10 −16 m 2 s −1 and L p = 10 −5 m, respectively, as in Section 2. Figure 2 plots the cumulative percentage of drug released (Equations (3.13) and (3.16)) for the instantaneous dissolution model described in Section 3.1. We consider four different ratios of initial drug concentration to solubility, C 0 /C s = 10, 5, 2, 1. Figure 2 demonstrates that by increasing the fraction C 0 /C s , we can significantly prolong the duration of release. We note that when C 0 /C s = 1, i.e. when the initial drug concentration is not above solubility, we obtain a release profile identical to that provided by (2.6) as expected. Figure 3 plots the cumulative percentage of drug released (Equation (3.25)) for the first-order reaction model described in Section 3.2. We consider three different values of K 1 , namely 10 −4 s −1 , 10 −5 s −1 and 10 −6 s −1 , which correspond to Damkohler numbers (Da = L 2 p K 1 /D p ) of 100, 10 and 1. We demonstrate that, as we reduce the Damkholer number (by reducing the value of K 1 ), the release profile can be significantly slowed down as diffusion becomes faster than reaction: in this case the system is thus reaction-limited. In Fig. 4, we display the cumulative percentage of drug released (Equation (3.31)) for the first-order reversible reaction model described in Section 3.3. We choose K 1 = 10 −6 s −1 and consider four different ratios of K 1 /K 2 , namely 100, 10, 1, 0.1 and demonstrate that in the case of a reversible reaction, the release can also be significantly slowed down. Furthermore, we note that for large K 1 /K 2 , the profile approaches that of the non-reversible reaction model. We note, too, that there is a qualitative difference in the shape of the profile of the three different models.

Mechanism of drug release: the inverse problem
In Sections 2 and 3 we provided a series of analytical solutions to models which encompass diffusionand dissolution-based drug release. Here we test our models on in vitro experimental data that we have obtained in our laboratory to try to ascertain the mechanism of release in a commercially available DES, as well as estimating the model parameter(s) via the inverse problem. The experiments involved measuring drug release from the Cypher DES. We have previously reported on the drug release from the Cypher DES (McGinty et al., 2013). Briefly, the experiments consisted of placing four Cypher DESs in separate sealed glass vials containing physiological release medium (phosphate buffered saline:ethanol (90:10)). The experiments were carried out at 37 • C. At several time points up to 60 days, each stent was removed and placed in a separate vial containing fresh release medium (to maintain perfect sink conditions), with the mass of drug in the original solution subsequently quantified using UV-spectroscopy. Solving an inverse problem consists of adjusting the parameters of a model function so as to best fit the data set. In this case, the data set should be the experimentally measured mass of drug released at various times. Adopting a least squares approach (and solving the resulting nonlinear equations via Newton's method), we find that drug release from the Cypher stent is best described by the diffusion only model given by equations (2.1-2.4a). The best fitting value of the diffusion coefficient is  Fig. 4. Release profiles for the case of reversible dissolution. Here K 1 = 10 −6 s −1 and we vary the value of K 2 . For K 1 /K 2 1, the effect of the reversible dissolution is relatively small, as expected. However, when K 1 /K 2 < 1 so that the backward reaction is faster than the forward reaction, the release may be significantly slowed down. D p = 6.3 × 10 −17 m 2 s −1 . Figure 5 displays a comparison between the experimentally measured cumulative percentage of drug released, and the solution of (2.6). The good agreement between the model and the experiments serves to demonstrate that diffusion is the dominant mechanism of release, at least for this particular stent. For the case of the instantaneous diffusion model, the best fit was obtained when C 0 /C s =1, i.e. when the initial drug concentration is not above solubility. In this case the solution is equivalent to (2.6). For the first-order reaction model, the best fit was obtained when the reaction rate, K 1 , was so large that the dissolution could be considered instantaneous. Similarly, for the first-order reversible reaction model, the best fit was obtained when K 1 was very large and K 2 very small, again suggesting dissolution can be regarded as instantaneous, while also indicating a negligible backward reaction.

Estimating physiological transport parameters
In many physiological processes, it is extremely difficult in practice to measure diffusion coefficients. This is in part, at least, because there is often a small but significant convection flow; this is particularly true near to arterial walls where there is a transmural pressure gradient and a consequent transmural flow. In addition, the species diffusing may be consumed or lost from the system through such vehicles as the vasa vasorum. In this section we indicate how a simple mathematical approach may allow the diffusion coefficient to be uncoupled from these other effects.

Drug transport through the arterial wall
Consider the initial value problem where C T is the drug concentration in the tissue, C 0 some constantly applied concentration, v the magnitude of the transmural velocity, D T the diffusion coefficient in the tissue and drug is removed from the system in proportion to α. For the case of the drug-eluting stent, the tissue is considered to be the region (0 < x < ∞) and while a constant applied concentration of drug at the polymer/media interface is somewhat unrealistic, it may, nonetheless provide a good approximation for early times. A more comprehensive model has been considered by Pontrelli and de Monte (2009), where the equation for drug transport in the tissue (5.1) has been coupled with a diffusion equation for the release of drug from the coating. Here, to effect comparisons with ex vivo experimental data, we restrict our attention to the equation in the tissue.
Taking Laplace transforms of (5.1) and making use of the boundary conditions leads tō Now let x/ √ D T = √ a and (v 2 /4D T ) + α = b, so that (5.2) can be written in the more concise form Applying Lemmas 1 and 2 (as detailed in Appendix C) directly provides two forms of the solution: where erfc is the complementary error function. Clearly, letting v, α → 0 in (5.4) and (5.5), returns the well known solution for the diffusion equation, that is The two forms of solution, (5.4) and (5.5), can be shown to be equivalent (see Appendix D) and further, a third equivalent form of solution is given by (5.7)

Inverse problem
The motivation for developing a solution to (5.1) was to allow the various model parameters to be inferred via the inverse problem. The data set, in this case, should be a series of experimental tissue concentration profiles where the drug transport is governed by convection as well as diffusion and drug is lost from the system. The data we would require to fit to our analytic solution do not appear to be available at present in the literature. However, Creel et al. (2000) provide concentration profiles based on experiments which, while admittedly do not account for drug loss (and thus render our parameter α redundant), do still allow us to estimate D T and v. Figure 6 shows experimentally measured tissue paclitaxel concentration, normalized with respect to the applied endovascular concentration. The data are obtained from experiments where arterial samples were perfused ex vivo for 15 min, 1 h and 4 h with a physiological transmural pressure gradient. Paclitaxel was applied to the endovascular aspect of the artery in buffer solution and drug distribution determined through en-face cryosectioning. In applying a least squares approach, starting guesses of v = 10 −  the parameter v converges to a value an order of magnitude higher while the parameter D T converges to a value an order of magnitude lower. The initial guesses along with an assumed arterial thickness of the order 10 −4 m (McGinty et al., 2010), results in a Peclet number O(1). However, the resulting Peclet number from the least squares analysis is O(100). This only serves to demonstrate the importance of accurate parameter estimation. While we have only been able to estimate v and D T from the experimental data available, the parameter α could also in principle be inferred if more appropriate experimental data were available.

Conclusions
This paper has been concerned with release mechanisms and parameter estimation in drug-eluting stent systems. Analytical solutions were obtained for mathematical models that were diffusion-and diffusiondissolution based. In addition, analytical solutions were also provided for dissolution kinetics involving first-order reversible and non-reversible reactions. Experiments were performed in vitro and the subsequent data employed, via an inverse problem, to determine which of the models best fitted the in vitro data. It was found that the simplest model, involving only diffusion, provided the best description of in vitro drug release from the Cypher stent. However, it is anticipated that the other models developed in this paper may well provide a better description of the release from different stent platforms. Finally, a simple model of the transport of drug through the arterial wall was formulated; its analytical solution was then employed using the existing data from Creel et al. (2000) to determine the underlying diffusion coefficient and the transmural velocity through the arterial wall.

Appendix A
In this appendix we consider the model defined in Section 3.1 and derive the expressions for the concentration of drug in the polymer and the cumulative fraction of drug released when the concentration of drug in the polymer has dropped below the solubility. The model is stated as

Separation of variables results in the solution
where B n (n = 1, 2, . . .) are constants to be determined. Imposing (A.4), multiplying through by cos ((m − 1/2)π x/L p ) and integrating from 0 to L p provides: Consider firstly the integral on the left-hand side of (A.6): on interchanging the integrals and performing the internal integration. Writing sine in terms of cosine in (A.7) and performing the integration results in Thus, returning to (A.6), and making use of orthogonality: we may write (A.9) as Substitution of (A.10) into (A.5) provides the solution as (3.15). Now, the cumulative fraction of drug released is given by where M (t) = L p 0 C p (x, t) dx is the mass of drug per unit area and M 0 = L p C 0 . Evaluating the integration, the solution is readily seen to be (3.16).

Appendix B
In this appendix we provide an outline of the solution of the model presented in Section 3.3. The model presented in Section 3.2 is then shown to be a special case and we deduce the solution by taking the limit as K 2 → 0. Consider the model of drug release by dissolution and diffusion as presented in Section 3.3: Substitution of (B.7) into (B.1) and taking Laplace transforms we find that We note that By writing the hyperbolic cosine in terms of exponentials it is readily seen that the roots of (B.13) must satisfy Γ (s)L p = iπ(n − 1 2 ), n = ±(1, 2, 3, . . .). (B.14) Thus we have a countably infinite number of singularities. We proceed to invert (B.11) by use of the Residue Theorem. Substituting the expression for Γ from (B.9) into (B.14) we see that the roots must satisfy the quadratic: The roots s 1 n and s 2 n are given by (3.30) from which we note that they all lie to the left of s = 0. The residues are thus given by (−1) n D p (2n − 1)π(s j n + K 1 ) 2 cos ((2n − 1)π x/2L p ) exp {s j n t} L 2 p s j n (s 2 j n + 2K 1 s j n + K 1 (K 2 + K 1 ))(K 2 + K 1 + s j n ) .