1. Introduction

Coronary artery disease (CAD) is a common cause of heart disease and heart attacks. It is caused by the buildup of atheroma, also known as plaque, on the inner walls of the coronary arteries. Drug-eluting stents (DESs) are currently one of the preferred treatment options of CAD. A DES generally consists of a metallic scaffold to hold the artery open and a polymer coating containing a drug that diffuses into its surroundings subsequent to deployment. The drug, which is usually an anti-proliferant, helps to prevent re-blockage of the artery due to restenosis. The drug is contained within the polymer coating and then slowly released into the arterial wall from the polymer source over a period on the order of 60-120 days [5,17]. The first DESs were designed with nondegradable polymer coatings; however, some of the newer DESs are manufactured with biodegradable polymer coatings [7,14]. In Figure 1, the deployment of a DES in a diseased coronary artery is schematically represented.

Mathematical modelling has a huge potential to inform both the design of drug-eluting stents and the choice of the appropriate DES for a specific lesion. For that potential to be realised mathematical models of drug-eluting stents must be capable of describing the phenomenon of restenosis and its modulation by the drugs released by the stent. In other words mathematical models of drug-eluting stents must model not just the chemistry of drug release and transport but also the biology of restenosis and the interaction between the two. In this series of two papers we present a modelling framework that can be used to address these questions. Our aim is firstly to show that such models can be created and can give results consistent with available experimental data, and secondly to show where additional experimental measurements would be most useful to enable the refinement of these models. This is the first of the two papers in which we use nondimensionalisation and asymptotic analysis to investigate models of the drug release into the artery wall and transport of drugs through the artery wall. In doing so we develop reduced models of the drug concentration within the artery wall which are used as inputs for the models of the biology of restenosis developed in the paper [19], which is currently in preparation.

Accurately modelling the kinetics of drug in-vivo subsequent to its release from a stent coating is very challenging because there are numerous factors that can affect the drug behaviour. For example, drug redistribution will depend on its diffusive and convective character in the artery wall, and the artery wall is known to contain three distinct substructures through its thickness [17,13]. Furthermore, compressed atherosclerotic plaque is likely to form part of the drug's tissue environment near to the inner wall of the artery because of the stent implantation procedure [9]. Also, the drug may bind specifically and non-specifically with receptors in the tissue, and specific binding in particular can have a very strong effect on drug kinetics [26,27,25]. Other complicating factors include the details of the construction and drug loading of the polymer coating, and drug washout through the inner and outer walls of the artery. Added to this, and for obvious reasons, there is a scarcity of experimental data for drug release from stents implanted in-vivo.

Mathematical models of varying complexity and sophistication have been proposed to model drug release from a DES. In a recent study, McGinty et al. [17] have developed a hierarchy of mathematical models to describe elution from stents that incorporate many of the phenomena referred to above; earlier references for stent modelling can also be found in this study. Bozsak et al. [6] developed a computational model for drug eluted from a drug-eluting stent into the arterial wall. The model took into account the multilayered structure of the arterial wall and incorporated a reversible binding model to describe drug interactions with the constituents of the arterial wall. They assumed drug transport to be purely diffusive in the polymer coating. The model considered here is closely related to the models described in Borghi et al. [5], Sakharov et al. [21], and Tzafriri et al. [24].

In [5,21], the transport of drug within the polymer is assumed to be dominated by diffusion while in the artery wall, the effects of both diffusion and reversible binding of drug to receptors were incorporated in the modelling. Only numerical solutions were calculated. The model in [24] included drug convection, diffusion and accounted for saturable binding of drug to both specific and non-specifics binding sites. However they neglected the actual geometry of the drug-eluting stent and approximated it using an equivalent phantom surface that elutes a defined drug load to the arterial lumen and wall. They provided both numerical and analytical results and also in-vivo experimential data.

McGinty et al. [15,16] presented a collection of models to describe drug release from the polymer including diffusion-based models and diffusion-dissolution-based models. When the initial concentration of drug in the polymer exceeds the solubility limit and the drug can only diffuse after it has dissolved, the models that incorporate both dissolution and diffusion are often used [8,10,4]. However, when the drug solubility is high or very low initial concentration of drug, a pure diffusion model is more appropriate [18]. Sirianni et al. [22] evaluated models that accounted for drug release from polymer coatings by different mechanisms. They observed that Fickian diffusion, dissolution and osmotic gradient models were capable of fitting the data equally well, and concluded that when the mechanism of drug release is not known, the simplest model with good predictive value is desired.

In this paper, we analyse a somewhat simple model to describe the redistribution of drug in a coronary artery wall subsequent to its elution from an implanted DES. We apply Fick's laws to model the evolution of the drug concentration in the polymer coating and in the tissue, and also include for the tissue the effects of drug convection, diffusion and a reversible binding of drug to specific receptors. Both numerical and analytical solutions are calculated. Although the model contains many simplifying assumptions, it is in our view capable of providing some useful order of magnitude estimates for the key quantities of interest. It will be seen that three independent small dimensionless parameters usually arise in the system, which complicates an asymptotic analysis, but does allow for useful qualitative information to be extracted. The advantage of a simple model is that we can obtain analytical results which can provide useful qualitative insights into the mechanisms governing release behaviour. From the analysis, we shall obtain the time scales over which the drug traverses the artery wall, empties from the polymer coating and resides in the arterial tissue. Also a formula for the total amount of drug in the artery tissue as a function of time will be derived, as well as a formula for the release profile from the polymer.

---

2. Model of drug release from DES

For simplicity, we consider a one-dimensional problem, and suppose that the polymer coating is located at $-L_p<x<0$, with $x=0$ giving the interface between the polymer and the artery wall, and $L_p, L_a$ denoting the thickness of the polymer coating and the artery wall, respectively. It is supposed that there is a stent strut located at $x=-L_p$ through which the drug cannot penetrate. We denote by $c_p(x,t)$ the concentration of drug in the polymer at penetration $x$ and time $t$, and suppose that this concentration is governed by Fick's law, so that (in dimensional variables):

$\begin{aligned} &\frac{{\partial c_p}}{{\partial t}} = D_p \frac{\partial^2 c_p}{\partial x^2} \text{ in } -L_p<x<0, t>0,\\ & \frac{{\partial c_p}}{{\partial x}} (-L_p, t) = 0 \text{ for } t \geq 0, \\ &c_p(x, 0) =c^* \text{ for } -L_p < x < 0, \label{eq:governing_polymer} \end{aligned}$

(1)

where $D_p$ is the constant diffusivity of the drug in the polymer coating, and $c^*$ is the uniform initial drug concentration in the polymer. It should be emphasised that all of the assumptions here concerning the drug and the polymer are certainly not true of all stents. For example, the polymer coating may not be uniformly loaded, and the drug concentration need not be below solubility throughout the coating [24,22]. Furthermore, some modern stent coatings are manufactured using biodegradable materials. However, it should also be remembered here that the manufacturer has control over the design of the polymer coating and the drug loading, and there is a case to be made for designing the system so that its behaviour may be adequately described by a simple mathematical model.

In the arterial tissue, it is supposed that the drug can associate with and dissociate from its specific binding sites, and that it can diffuse in its free form and be convected by the outward movement of plasma through the arterial wall. We also suppose that the effect of non-specific binding is negligible. We denote by $C, B, A,$ a free molecule, a specific binding site, and a drug-specific binding site complex (bound drug), respectively. The reversible binding reactions can then be represented by:

$C + B\underset{{k_{\rm {off}}}}{\overset{{k_{\rm {on}} }}{\rightleftharpoons}}A,$

where $k_{\rm {on}}, k_{\rm {off}}$ are rate constants. We denote by $a(x,t), b(x,t), c(x,t)$ the concentrations of $A,B,C$, respectively, at location $x$ and time $t$. Following [17,24], we suppose that the convection velocity of the plasma is constant and denote it by $V_a$. The governing equations for $a,b,c$ are now:

$\begin{align}&\frac{{\partial a}}{{\partial t}} = k_{\rm {on}} bc - k_{\rm {off}} a, 0<x<L_a, t > 0, \notag \\ &\frac{{\partial b}}{{\partial t}} = - k_{\rm {on}} bc + k_{\rm {off}} a, \ 0<x<L_a, t >0, \notag\\ &\frac{{\partial c}}{{\partial t}} + V_a \frac{{\partial c}}{{\partial x}}= D_a\frac{\partial^2 c}{\partial x^2} - k_{\rm {on}} bc + k_{\rm {off}} a, 0<x<L_a, t >0, \label{eq:governing_artery} \\ &c(L_a,t)= 0 \text{ for } t > 0, \notag \\ &a(x,0)=0, b(x,0)=b^*, c(x,0)=0 \text{ for } 0< x < L_a, \notag \end{align}$

(2)

where $D_a$ is the diffusivity of the drug in the arterial tissue, and $b^*$ is the equilibrium concentration of specific binding site for the drug in the arterial tissue. It is noteworthy in this model that the drug can only diffuse in its free form, and that both the binding sites and the drug-binding site complexes are immobile. Therefore binding will clearly have an effect on the rate at which the drug can move through the tissue. The first three equations in (2) may be written in the equivalent form

$\begin{align}&\frac{{\partial }}{{\partial t}} (a+c) = D_a\frac{\partial^2 c}{\partial x^2} - V_a \frac{{\partial c}}{{\partial x}}, 0<x<L_a, t > 0, \notag \\ &\frac{{\partial }}{{\partial t}} (a+b) = 0, 0<x<L_a, t >0, \label{eq:governing_artery_sum} \\ &\frac{{\partial a}}{{\partial t}} = k_{\rm {on}} bc - k_{\rm {off}} a, 0<x<L_a, t >0, \notag \end{align}$

(3)

where, for example, equation (3)₁ is obtained by forming (2)₁+(2)₃.

The problem is completed by imposing continuity in the drug concentration and drug flux at the polymer-artery wall interface, so that

$c_p(0^-,t)=c(0^+,t), \left( - D_p \frac{\partial c_p}{\partial x} \right)_{x=0^{-}} = \left(- D_a \frac{\partial c}{\partial x}\right)_{x=0^{+}} + \left( V_a c \right)_{x=0^{+}} \text{ for } t \geq 0. \label{eq:interface_boundary}$

(4)

---

2.1. Non-dimensionalization

We introduce non-dimensional variables as follows

$\bar{t}=\dfrac{t}{(L_p^{2}/D_p)}, \bar{x}=\frac{x}{L_p}, \bar{a}=\dfrac{a}{b^*}, \bar{b}=\dfrac{b}{b^*}, \bar{c}=\frac{c}{c^*}, \bar{c}_p=\frac{c_p}{c^*},$

to obtain the following dimensionless equations (dropping the over-bars for convenience)

Polymer coating:

$\begin{aligned} & \frac{{\partial c_p}}{{\partial t}} = \frac{\partial^2 c_p}{\partial x^2}, -1<x<0, t>0, \\ & \frac{{\partial c_p}}{{\partial x}}(-1,t) =0 \;{\text{ for }} t\geq 0, \\ & c_p(x,0) =1\;{\text{ for }} -1<x<0; \label{eq:nondimension_polymer} \end{aligned}$

(5)

Arterial tissue:

$\begin{aligned} & \varepsilon \frac{{\partial}}{{\partial t}} (\eta a +c) = \frac{\partial^2 c}{\partial x^2} - \frac{Pe}{L} \frac{{\partial c}}{{\partial x}}, 0<x<L, t > 0, \\ & a+b=1,\; \delta \frac{{\partial a}}{{\partial t}} = \frac{K_b}{\eta}bc - a, 0<x<L, t > 0,\\ & c (L, t ) = 0\;{\text{ for }} t \geq 0, \\ & c(x,0)=0\;{\text{ for }} 0<x< L; \label{eq:nondimension_artery} \end{aligned}$

(6)

Polymer/artery wall interface:

$c_p(0^-,t)=c(0^+,t), \left( -\varepsilon \frac{\partial c_p}{\partial x} \right)_{x=0^{-}} = \left(- \frac{\partial c}{\partial x} \right)_{x=0^{+}} + \left( \frac{Pe }{L} c \right)_{x=0^{+}} {\text{ for }} t \geq 0, \label{eq:nondimension_interface}$

(7)

where

$L= \frac{L_a}{L_p}, \varepsilon= \frac{D_p}{D_a}, Pe= \frac{V_a L_a}{D_a}, \eta=\frac{b^*}{c^*}, K_b=\frac{k_{\rm {on}} b^*}{ k_{\rm {off}}}, \delta= \frac{D_p}{L_p^2 k_{\rm {off}}}. \label{eqB:non_parameters}$

(8)

Here $L$ is the ratio of the artery and polymer coating thicknesses; $\varepsilon$ is the ratio of the diffusivities in the polymer coating and in the artery tissue; $\eta$ is the ratio of initial concentration of binding sites to initial drug concentration; $Pe$ is the Peclet number which determines the relative importance of drug transport by advection and by diffusion (if $Pe$ is small diffusion dominates over advection); $K_b$ gives the binding constant for drug to binding site ($K_b\gg 1$ and $K_b \ll 1$ correspond to strong retention and weak retention, respectively, of the drug by the binding site); and $\delta$ measures the ratio of the dissociation time scale to the diffusion time scale.

It is seen that for some drug/tissue systems of particular interest, the diffusion time scale is much longer than the two time scales associated with specific binding [27]. For example, taking $D_p=1.0\times 10^{-10}$ mm$^{2}$s$^{-1}$ (Table 3) as a representative diffusivity for a drug in polymers and the polymer thickness to be $L=0.01$ mm (Table 1), we calculate a typical diffusion time scale to be $L_p^{2}/D_p\approx 11.6$ days. However, the time scales associated with the binding reactions are frequently much shorter. For Rapamycin, taking $k_{\rm {off}} = 0.096$ min$^{-1}$ and $k_{\rm {on}} = 4.8 \times 10^{7}$ M$^{-1}$ min$^{-1}$ for the binding constants and $b^*= 3.3 \times 10^{-6}$ M for the initial concentration of binding sites [27] gives the binding time scales $1/k_{\rm {off}} \approx 10.4$ minutes and $1/(k_{\rm {on}} b^*) \approx 0.4$ s. Therefore in this paper, we confine our discussion to cases for which the diffusion time scale is much longer than the time scales associated with the binding reactions, so that:

$L_p^{2}/D_p\gg \max \left \{1/(k_{\rm {on}}b^{*}),1/k_{\rm {off}} \right \},$

which implies that $\delta \ll \min(K_{b},1)$. Hence, we neglect the term involving $\delta$ in the second equation in (6)₂, to obtain

$\eta a = K_{b} bc. \label{eq:reduced}$

(9)

The initial-boundary value problems (5), (6), (7) and (9) can be numerically integrated using the command pdepe in the mathematical software package MATLAB. The pdepe solver implements the method of lines to convert the partial differential equation to a set of ordinary differential equations (ODEs) using a second-order accurate spatial discretization. The resulting ODEs are integrated to obtain approximate solutions at various times. An implicit time-stepping finite difference algorithm is used with the time step determined automatically and adaptively by the ODE solver (ode15s in MATLAB).

---

2.2. Designing a stent coating system

In Tables 1, 2, 3, 4, and 5, values for the parameters appearing in the model above are shown. The point of view taken in the current analysis is that of a stent manufacturer who wishes to design a drug loaded polymer coating. Such a person would prefer to restrict the design space to a region where the performance of the stent could be robustly predicted by a simple model, if that restriction allows viable stents to be designed. We show that by restricting our attention to a system in which the polymer is monolithic and nondegradable, and that the drug is uniformly dispersed throughout the polymer bulk at a concentration below solubility, we are able to predict the performance of the system from a simple model. The drug delivery industry has extensive experience in designing monolithic polymeric devices.

The manufacturer wishes to design the system so that a sufficient amount of drug is released into the artery wall for a sufficient period to prevent restenosis. More precisely, the manufacturer wishes to design the system so that a significant proportion of the specific binding sites in the artery wall are occupied by the drug for a period of some months subsequent to the stent being implanted.

The task then is to identify a parameter regime for the governing equations that achieves the stated goal subject to the constraints. From (8), it is seen that there are five independent dimensionless parameters that can in principle be independently varied to tune the system. However, two of these parameters, $K_b$ and $Pe$, are largely determined by the nature of the drug, and since most modern stent systems use either sirolimus (or one of its close relatives) or paclitaxel, there is not much scope for varying these. The parameter $L$ may be varied by changing the thickness of the polymer. However, the parameters $\varepsilon=D_p/D_a$ and $\eta=b^*/c^*$ are probably the most convenient to use to optimise the system since the diffusivity and the drug-loading for the polymer are readily changed.

We first turn our attention to the selection of appropriate values for the parameter $\varepsilon$. A drug-eluting stent implanted in a coronary artery is required to release drug for a period of at least a few months subsequent to its deployment in order to prevent restenosis [7,14]. Hence, since we are assuming here that diffusion is the only mechanism for drug transport in the polymer coating, it is required that the drug diffusion time scale in the polymer, $L_p^2/D_p$, should be of the order of some weeks. For the sake of definiteness, we suppose that:

$L_p^2/D_p \sim 2 \mbox{ weeks.} \label{eq:Lp_Dp}$

(10)

Inspecting Table 1, it is seen that $L_p\sim 10$ $\mu$m, and (10) then implies that for the drug to release from the polymer over an appropriate time scale, the polymer should be fabricated so that:

$D_p \sim 10^{-10} \mbox{ mm}^2/\mbox{s} \mbox{ or smaller.} \label{eq:Dp_approx}$

(11)

It is noteworthy that the diffusivities $D_p$ displayed in Table 3 are $O(10^{-10})$ mm$^2$/s or smaller.

Inspecting the data in Table 2, it is seen that if $D_p\sim 10^{-10}$ mm$^2$/s, then:

$\varepsilon = \frac{D_p}{D_a} \sim \left\{ \begin{array}{ll} 10^{-6}&\mbox{for rapamycin,} \\ 10^{-4}&\mbox{for paclitaxel,} \\ 10^{-5}&\mbox{for heparin,} \\ 10^{-5}&\mbox{for dextran.} \end{array} \right. \label{eq:epsilon}$

(12)

If $\varepsilon$ is of order $10^{-4}$ or smaller, then it is clear from Table 5 that for commercially available stenting systems, we have:

$\varepsilon \ll 1/K_b \ll 1/L \ll 1. \label{eqB:epsilon}$

(13)

It would seem from this that $\varepsilon \rightarrow 0$ is a sensible limit to consider to analyse the behaviour of (5), (6) and (7). This is indeed the case, but care must be taken to ensure that combinations of $K_b$ and $L$ do not arise in the system which may interfere with the accuracy of the results. This follows from the fact that in order to preserve the accuracy of our results, no quantity may arise in our equations (6) that is comparable to or smaller than $\varepsilon$. This could arise because $1/K_b, 1/L$ are both small. In fact, we shall show that we require $\varepsilon \ll 1/L^2$ in the next section.

We now turn our attention to the selection of the parameter $\eta$, which corresponds to choosing a drug loading concentration for the polymer. After briefly investigating various limits for the governing equations, we propose the following:

$\eta/K_b = O(\varepsilon),$

with $\varepsilon \ll 1$ as explained above, or, in dimensional terms:

$c^*/K_{\rm {D}} = O\left(D_a / D_p\right), \text{ with } D_p \ll D_a,$

where $K_{\rm {D}}= k_{\rm {off}}/k_{\rm {on}}$ is the dissociation constant for the drug from its specific binding sites. It is noteworthy that with this choice the mass transfer Fourier number in the artery tissue will be $Fo=O(\frac{1}{\varepsilon L^2}) \gg 1$. Therefore we can assume that the drug concentration in the artery wall is approximately independent of space. In our second paper [19], the drug concentration is assumed to be constant in the intima for small timescales.

We now justify this choice by carrying out an asymptotic analysis of the governing equations in the limit $\varepsilon \rightarrow 0$, with $\eta/K_b = O(\varepsilon)$.

---

2.3. Preparing the polymer coating: $c^*/K_{\rm {D}} = O\left(D_a / D_p\right), D_p/D_a \ll 1$

Writing $\eta/K_b= \mu \varepsilon$ with $\mu =O(1)$, so that the second equation in (6)₂ becomes:

$\mu \varepsilon a= bc.$

Since $K_b \gg 1$ for most of the drugs of interest, we shall also make the choice $K_b= O(\varepsilon^{-1/2})$ here, which implies that $\eta= O(\varepsilon^{1/2})$. However, this choice is not particularly significant since most of the more important results shall be derived below depend only on the requirement that $\eta/ K_b = O(\varepsilon), \varepsilon \ll 1$. We write $\eta = \varepsilon^{1/2} \eta^*$ with $\eta^*=O(1)$.

There are two time scales to consider in the limit $\varepsilon \rightarrow 0$: a short time scale $t=O(\varepsilon)$, and a longer time scale $t=O(1)$.

---

2.3.1. Short time scale, $t=O(\varepsilon)$

In dimensional terms, this time scale is given by $t=O(L_p^2/D_a)$. Writing $t= \varepsilon \hat{t}$, it is found that $c_p \sim1$ in $-1<x<0, \hat{t}=O(1)$.

In $0<x<L, \hat{t}=O(1)$, we have $a=1+ O(\varepsilon^{1/2}), b, c= O(\varepsilon^{1/2})$, and we pose

$a\sim 1+ \varepsilon^{1/2} \hat{a}_{0}(x,\hat{t}), b\sim \varepsilon^{1/2} \hat{b}_{0}(x,\hat{t}), c \sim \varepsilon^{1/2} \hat{c}_{0}(x,\hat{t}) \mbox{ as } \varepsilon \rightarrow 0,$

to obtain

$\hat{a}_{0}(x,\hat{t}) = - \frac{\mu}{ \hat{c}_{0}(x,\hat{t})}, \hat{b}_{0}(x,\hat{t})=\frac{\mu}{ \hat{c}_{0}(x,\hat{t}) },$

and

$\frac{\partial \hat{c}_{0}}{\partial \hat{t}} = \frac{\partial^2 \hat{c}_{0}}{\partial x^2} - \frac{Pe}{L} \frac{\partial \hat{c}_{0}}{\partial x}, \; 0<x<L, \hat{t}>0. \label{eq:artery_short1}$

(14)

Recalling that $L \gg 1$, it is clear that the diffusion term dominates on the right hand side of equation (14); so that the drug traverses the artery wall on the diffusion time scale $\hat{t}= O(L^2)$, or $t=O(\varepsilon L^2)$, which corresponds in dimensional terms to $t= O(L_a^2/D_a)$, the diffusion time scale for the free drug in the artery wall. We deduce from this that in order for the asymptotic approximations to be valid, it is required that $\varepsilon \ll 1/L^2$, or, in dimensional terms:

$D_p \ll \frac{L_p^2}{L_a^2}D_a, \text{ or } \frac{L_p^2}{D_p} \gg \frac{L_a^2}{D_a}, \label{eqB:asymptotic_criterion}$

(15)

which implies that the drug diffusion time scale in the polymer must be much longer than the free drug diffusion time scale in the artery wall. For rapamycin, this implies that $L_p^2/ D_p \gg 1$ hours, and for paclitaxel it requires that $L_p^2/ D_p \gg 2$ days. However, we are insisting that $L_p^2/D_p \sim 2$ weeks here, as previously discussed, so these criteria are met.

The specification of the problem for $\hat{c}_{0}(x,\hat{t})$ requires the derivation of a boundary condition for $\hat{c}_{0}$ on $x=0$. This is obtained by noting that there is a boundary layer near the surface of the polymer where it interfaces with the artery tissue. This layer is located at $\hat{x}=O(1), \hat{x}<0$ where $x=\varepsilon^{1/2} \hat{x}$, and in $\hat{x}, \hat{t} =O(1)$, we pose $c_p \sim \hat{c}_{p{0}} (\hat{x}, \hat{t})$ to obtain

$\begin{align}& \frac{ \partial \hat{c}_{p 0}} {\partial \hat{t}}=\frac{\partial^2 \hat{c}_{p 0}}{\partial \hat{x}^2}, \; -\infty<\hat{x}<0, \hat{t}>0, \notag \\ &\hat{c}_{p0}(\hat{x}, \hat{t}) \rightarrow 1\;{\text{ as }} \hat{x} \rightarrow -\infty, \hat{t} \geq 0, \label{eq:artery_short2} \\&\hat{c}_{p0}(0^-,\hat{t})=0\;{\text{ for }} \hat{t}\geq 0, \notag \end{align}$

(16)

and this self-similar problem has solution

$\hat{c}_{p0}(\hat{x},\hat{t})= - \mbox{erf} \left( \frac{\hat{x}}{2 \sqrt{\hat{t}}} \right). \label{eq:short_cp0}$

(17)

The perfect sink boundary condition on $x=0^-$ in (16) is noteworthy because it implies that at leading order the problem for the drug concentration in the polymer decouples from that in the tissue. It also implies that in-vitro experimental release studies should adequately mimic the release behaviour from the polymerin on this time scale (and for longer times too, as we shall see). It follows from (17) that, in dimensional unscaled variables, the fraction of drug released from the polymer for short times is approximated by

$\frac{M(t)}{M(\infty)} \sim 2 \sqrt{\frac{D_p t}{\pi L_p^2}} \;{\text{ for }} t=O\left(L_p^2/D_a\right). \label{eqB:Mt_shorttime}$

(18)

The boundary condition for $\hat{c}_{0}$ on $x=0$ is obtained by matching the drug fluxes across $x=0$. (see equation (7)):

$\lim\limits_{\hat{x}\rightarrow 0^-} \left( - \frac{ \partial \hat{c}_{p 0}} {\partial \hat{x}} \right) = \lim\limits_{x\rightarrow 0^+} \left( - \frac{ \partial \hat{c}_{0}} {\partial x} + \frac{Pe}{L} \hat{c}_{0}\right),$

to obtain:

$-\frac{\partial \hat{c}_{0}}{\partial x}(0^+,\hat{t}) + \frac{Pe}{L} \hat{c}_{0}(0^+,\hat{t}) = \frac{1}{ \sqrt{\pi \hat{t}}} \;{\text{ for }} \hat{t} \geq 0. \label{eq:short_boundary_matching}$

(19)

Combining (14) and (19) yields a nonlinear initial boundary value problem which can in principle be solved for $\hat{c}_{0}$. However, the solution to (14) and (19) is not pursued here. In Figure 3, we plot some numerical solutions to (5-7) for $\varepsilon = 10^{-6}$ and $t=O(\varepsilon L^2)$ and for parameter values appropriate to rapamycin. It is clearly seen in this figure that the drug traverses the artery wall on this time scale, and that there is a sharp diffusion front tracking to the right.

---

2.3.2. Long time scale, $t=O(1)$

This is the time scale over which the drug empties from the polymer coating, and in dimensional terms, it is given by $t=O(L_p^2/D_p)$. In $-1<x<0, t=O(1)$, we pose $c_p \sim c_{p 0}(x,t)$ as $\varepsilon \rightarrow 0$, to obtain:

$\begin{align}& \frac{ \partial c_{p 0}} {\partial t}=\frac{\partial^2 c_{p 0}}{\partial x^2}, \; -1<x<0, t>0, \notag \\ & \frac{ \partial c_{p 0}} {\partial x} (-1,t)=0\;{\text{ for }} t\geq 0, \label{eq:polymer_long} \\ &c_{p0}(0^-, t)=0\;{\text{ for }} t \geq 0, \notag \\&c_{p0}(x,0)=1\; {\text{ for }}-1<x<0. \notag \end{align}$

(20)

Notice that on this time scale, we also have a perfect sink boundary condition for the drug on $x=0^-$. This is obtained by matching with the free drug concentration at $x=0^+$; see below. Hence, the leading order problem for the drug concentration in the polymer again decouples from the leading order problem in the tissue. Solving (20) yields:

$c_{p0}(x,t) = -\frac{4}{\pi}\sum\limits^{\infty}_{n=1}\frac{1}{2n-1}\sin \left( \frac{(2n-1)\pi x}{2} \right) \exp \left (-\frac{(2n-1)^2 \pi^2 t}{4} \right ). \label{eq:long_cp0}$

(21)

The fraction of the total drug released in $t=O(1)$ is easily calculated using this expression; combining the results of this calculation with (18) now gives (using dimensional variables):

$\frac{M(t)}{M(\infty)} \sim \left\{ \begin{array}{ll} 2 \sqrt{\frac{D_p t}{\pi L_p^2}}& {\text{ for }} t=O\left(L_p^2/D_a\right), \\ 1 - \frac{8} {{\pi ^2 }}\sum\limits_{n = 1}^\infty {\frac{1} {{(2n - 1)^2 }}\exp \left( { - \frac{{(2n - 1)^2 \pi ^2 D_pt}} {{4L_p^2 }}} \right)}& {\text{ for }} t=O\left(L_p^2/D_p\right), \end{array} \right. \label{eqB:Mt_short_long}$

(22)

In $0<x<L, t=O(1)$, we have $a,b=O(1), c=O(\varepsilon)$, and pose:

$a\sim a_{0}(x,t), b\sim b_{0}(x,t), c \sim \varepsilon c_{0}(x,t) \mbox{ as } \varepsilon \rightarrow 0,$

to obtain

$a_{0}(x,t) = \frac{c_{0}(x,t)}{\mu + c_{0}(x,t)}, b_{0}(x,t) =\frac{\mu}{\mu + c_{0}(x,t)}. \label{eq:long_ab0_c}$

(23)

Recalling that $\eta= O(\varepsilon^{1/2})$, as $\varepsilon \rightarrow 0$ we have that

$\begin{align}&\frac{ \partial^2 c_{0}} {\partial x^2} -\frac{Pe}{L} \frac{\partial c_{0}}{\partial x} =0, \; 0<x<L, t>0, \notag \\ &- \frac{ \partial c_{0}} {\partial x}(0^+,t) +\frac{Pe}{L} c_{0}(0^+,t) = \gamma(t)\;{\text{ for }} t\geq 0, \label{eq:artery_long} \\ &c_{0}(L, t)=0\;{\text{ for }} t \geq 0, \notag \end{align}$

(24)

where

$\gamma(t) = - \frac{ \partial c_{p0}} {\partial x}(0^-,t) = 2 \sum\limits^{\infty}_{n=1} \exp \left (-\frac{(2n-1)^2 \pi^2 t}{4} \right ), \label{eq;gamma}$

(25)

is the leading order flux of drug from the polymer into the tissue across $x=0$ for $t=O(1)$. These equations imply that the behaviour in the arterial tissue is approximately quasistatic, with the time dependence only entering via a time-varying flux of drug from the polymer coating, $\gamma(t)$. Integrating (24) leads to the following approximation

$c(x,t) \sim \frac{\varepsilon \gamma(t)L}{Pe}\left( 1- e^{-Pe(1-x/L)} \right), \label{eq:long_c0}$

(26)

so that

$\begin{aligned} &a(x,t) \sim \frac{\gamma(t)L \left( 1- e^{-Pe(1-x/L)} \right)}{\mu Pe +\gamma(t)L \left( 1- e^{-Pe(1-x/L)} \right)}, \\ & b(x,t) \sim \frac{\mu Pe}{\mu Pe +\gamma(t)L \left( 1- e^{-Pe(1-x/L)} \right)}; \label{eq:long_a0} \end{aligned}$

(27)

similar forms have recently been noted by Tzafriri et al. [24]. Hence, at the midpoint of the artery wall, we have

$a(L/2,t) \sim \frac{\gamma(t)L \left( 1- e^{-Pe/2} \right)}{\mu Pe +\gamma(t)L \left( 1- e^{-Pe/2} \right)}.$

Recalling that $1/\varepsilon \gg L \gg 1$, we have that

$a(L/2,t) \sim 1 \text{ for } \gamma(t) = O(1),$

so that we have approximately full occupancy of the specific binding sites at the centre of the artery wall on the time scale $\gamma(t)= O(1)$, or $t=O(1)$. In dimensional terms, this corresponds to $t = O(L_p^2/D_p)$, the diffusion time scale in the polymer which can be tuned by the manufacturer by making appropriate choices for $L_p$ and $D_p$. This suggests that the parameter regime we have chosen for the design of the coating system should produce satisfactory results since it predicts that a significant proportion of the specific binding sites in the artery wall will be occupied by the drug for a period of some months subsequent to the stent being implanted. In other words, it predicts that sufficient drug will be having its physiological effect in the artery tissue over the time scale of interest.

In Figure 4, we plot some numerical solutions to (5-7) for $\varepsilon = 10^{-6}$ and $t=O(1)$. It is confirmed in this figure that to a good approximation the behaviour is indeed quasistatic in the artery tissue, with the shape of the profiles scarcely changing for successive times. We also note that a significant proportion of the specific binding sites away from the outer wall of the artery are occupied for $\gamma(t)= O(1)$.

The average occupancy of the specific binding sites over the thickness of the artery wall is given by

$m_a(t)= \frac{1}{L} \int_0^L{a(x,t) dx}.$

In Figure 5 (a), we plot this quantity as a function of time for $\varepsilon= 10^{-6}$ and various values of $\eta$, with the other parameter values being appropriate for sirolimus, and $a(x,t)$ is calculated numerically from the initial boundary value problem (5-7). The rapid climb in the profiles near $t=0$ corresponds to the initial layer $t=O(\varepsilon L^2)$; this rapid behaviour has been observed in experimental studies [24]. We also note that a significant proportion of the binding sites are occupied for a period of a few months.

Using (27)₁, we obtain the useful approximation for $t=O(1)$

$m_a(t) \sim \frac{\gamma(t)L + \mu \ln \left( \frac{1}{1+ \frac{\gamma(t) L}{\mu Pe}\left(1 - e^{-Pe} \right) } \right)}{\mu Pe + \gamma(t) L} \text{ for } t\gg \varepsilon L^2. \label{eq:ma}$

(28)

In Figure 5 (b), we compare this result to the numerical profiles for various values of $\eta$ for the first month. There is a good agreement for $t=O(1)$, which in dimensional terms corresponds to $t=O(L_p^2/D_p)$, as expected. The sub window in Fig. 5 shows the relative errors between the asymptotic and numerical solutions, $(m_{a,\text{numerical}}-m_{a,\text{asymptotic}})/m_{a,\text{numerical}}$. The discrepancies increase with increasing time due to the fact that the asymptotic approximations are just to the leading order and cannot be expected to match the true solution exactly.

In addition, the approximation for the average of free drug over the thickness of the artery wall is calculated using (26)

$m_c(t)= \frac{1}{L} \int_0^L{c(x,t) dx} \sim \frac{\varepsilon\gamma(t)L}{Pe^2} \left(Pe -1 + e^{-Pe} \right) \text{ for } t\gg \varepsilon L^2, \label{eq:mc}$

(29)

then the average average mass of drug in the artery wall is

$m_T(t) = \eta m_a(t) + m_c(t). \label{eq:mT}$

(30)

The asymptotic analysis has justified the selection of parameters with $\eta/K_b = O(\varepsilon)$ and $\varepsilon \ll 1$, or in dimensional terms, $c^*/K_{\rm {D}} = O\left(D_a / D_p\right)$ and $D_p \ll D_a$. It implies that when preparing the polymer coating, the diffusion of drug in the polymer should be much slower than in the artery tissue and the ratio of initial drug concentration to the dissociation constant should be proportional to the ratio of the drug diffusivity in the artery tissue to the diffusivity in the polymer. Then we can expect that a significant proportion of the specific binding sites in the artery wall are occupied by the drug for a sufficient amount of time for the therapy to be effective.

---

2.4. Comparison with experimental data

We now compare the modelling results with in-vivo experimental data of drug elution from Cypher[24] and Xience V stents [20] implanted in porcine coronary arteries. We fit the expression (22)₂ for $t=O\left(L_p^2/D_p\right)$ to the experimental data using the method of least squares, with $D_p$ being the only unknown parameter. The nonlinear equation for $D_p$ arising from this fitting procedure was solved numerically using the mathematical package MAPLE, and the estimates for $D_p$ obtained from the fitting procedure are given in Table 3. The good agreement between the model results and the experimental data in Figure 6 demonstrates that diffusion is the dominant mechanism of release, at least for these particular stents. Therefore the simple diffusion model is sufficient here to describe drug release from the polymer coating. Figure 7 also shows a good agreement between the average mass of drug in the artery wall calculated by (30) and in-vivo experimental data. The fitting errors are estimated using the mean squared error (MSE) expression.

---

3. Conclusion

In this paper, we considered a simple model to provide an elementary description of drug release into artery tissue from an implanted stent. The model tracks the evolution of the concentration of both free drug and bound drug in the tissue. In the current work, we have only considered a one dimensional model, the structure of the artery wall was assumed to be homogeneous, and the effect of the non-specific binding site was neglected. Also a perfect sink condition was used at the distal end of the artery wall. In addition, we have only considered the case of a polymer that is monolithic and nondegradable, and where the drug is uniformly dispersed throughout the polymer bulk at a concentration below solubility. Although we have used many simplifying assumptions, as we demonstrate by comparing solutions to experimental data, the model is capable of providing some useful order of magnitude estimates for the key quantities of interest. When designing a drug-eluting stent system, we expect to have a sufficient amount of drug that is released into the artery tissue for a sufficient amount of time to prevent restenosis. A parameter regime is identified to optimise the system when preparing the polymer coating based on the model. It is shown that with the chosen parameter regime for the design of the coating system, a significant proportion of the specific binding sites in the artery wall are occupied by the drug for a period of some months subsequent to the stent being implanted. The model was evaluated by comparing with in-vivo experimental data and good agreement was found. In addition, we found that the mass transfer Fourier number in the artery tissue, $Fo=O(\frac{1}{\varepsilon L^2})$, is large. We therefore can assume a homogeneous distribution of drug when developing models of the interaction of the drug with the cellular processes leading to restenosis, the subject of the second paper in this sequence [19].

---

Acknowledgments

We gratefully acknowledge the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland (SFI) Investigator Award 12/IA/1683. Dr Vo also thanks the New Foundations Awards 2013 and 2014 from the Irish Research Council.

---