1. Introduction

Arterial disease, in particular atherosclerosis, is a significant cause of mortality and morbidity in the western world [1]. A standard treatment for the acute form of the disease, in which an artery is almost entirely occluded by atheroma, is percutaneous balloon angioplasty in which the narrowed artery is opened by inflating a balloon inside the region. To prevent restenosis (re-narrowing of the artery) a drug eluting stent-a metal framework with a drug infused polymer layer on its surface-is often deployed. The stent mechanically holds the artery open, while the drug released from the surface prevents inflammation and cell proliferation from inducing restenosis.

While percutaneous balloon angioplasty with drug eluting stent implantation is an effective treatment for acute atherosclerosis it has long been suspected that there is scope to greatly increase the efficiency of the treatment [13]. These gains could be driven by patient specific modelling of the disease and its treatment allowing personalised treatments to be developed. However, there is a significant gap in the models that have been developed.

Modelling studies of arterial disease have focussed on two areas. These are: (1) Studies of the fluid mechanics of blood flow within the artery before and after a stent has been deployed, motivated by the observation that wall shear stress plays a critical role in the initiation of arterial disease (e.g. [6], [9], [12], [17]); and (2) Studies of the transport of drugs within the arterial wall motivated by a need to better understand the dose and spatial distribution of the drug delivered to the artery wall (e.g. [8], [13], [14], [18], [22], [25]). Both effects in combination have been considered as well (e.g. [3], [4]). In particular, the readers attention is drawn to Part Ⅰ of this paper, [21]. Less attention has been directed to the cellular response to treatment, in particular the growth response of cells to stent implantation and drug release. One exception to this is agent based lattice modelling of the cellular structure of the arterial wall [19,5] embedded in a multiscale model of restenosis. This discrete treatment of individual cells contrasts with the continuum approach presented here.

In order for a model to be used to optimise a drug eluting stent design it must be capable of modelling the process the stent is deployed to prevent: restenosis. Current models on their own are not enough because they do not directly address the question of restenosis. To directly model restenosis, these models must be coupled to models of the cellular processes occurring within the artery wall that lead to restenosis. That is the focus of this paper.

Experimental results which can be used to validate such a model are very sparse. The only relevant study is that by [7], who studied the neointimal thickening of a healthy rabbit iliac artery under implantation of a bare metal stent and a drug eluting stent. The results of the study are summarised in Figure 1.

In this paper we develop models of the biological response of the cells within the arterial walls to stent deployment and drug elution. We do so by adapting models already developed in the context of tumour growth, prompted by the fact that some of the drugs used in drug eluting stents (e.g. Paclitaxel) are also used as anti-tumour drugs [11]. Our aim in doing so is to show that such models are consistent with the (relatively sparse) experimental data available and to encourage experimentalists to carry out experiments that could be used to further validate and refine such models.

As figure 1 shows, there is an inflammation response to stent implantation leading to a thickening of the intimal layer of the artery. This response can however be substantially reduced by the release of drugs. Thus the mathematical model we develop must be able to describe cell proliferation via a thickening of the artery wall and its reduction by the action of drugs.

---

2. A cellular model

To develop a mathematical model of the cellular response of the artery wall to stent implantation and drug elution we adapt ideas from two models of tumour growth. The first is a model of the response of breast and ovarian cancer to Paclitaxel (as mentioned above, a drug also used in stents) by [16]. This model employs a pair of differential equations to describe the dynamics of quiescent, $Q$, and proliferative, $P$, populations of cells. In this treatment, $Q$-cells do not multiply, while $P$-cells multiply at some known rate. Thus, the growth of a tumour depends on both the rate at which $P$-cells multiply and the relative fractions of $P$- and $Q$-cells, as shown in Figure 2. The action of Paclitaxel is modelled as selectively killing cells in the proliferative phase.

The second model is of tumour cords: cylindrical layers of tumorous tissue surrounding blood vessels [2]. This model describes the effect of growth on the increasing thickness of the layer via a radial velocity field.

---

2.1. The proposed model

Our mathematical model takes the form

$\frac{\partial P}{\partial t} + \frac{\partial}{\partial x}(uP) = (\gamma - \alpha - \lambda_{P} -\mu_{P})P + (\beta + \eta - \mu_{Q})Q \label{eq:prolif_cellMass_Dim}$

(1)

$\frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}(uQ) = (\alpha)P - (\beta + \eta + \lambda_{Q} - \mu_{Q})Q \label{eq:qui_cellMass_Dim}$

(2)

$P + Q + E = 1$

(3)

where:

$P(x,t)$ = proliferating cell fraction

$Q(x,t)$ = quiescent cell fraction

$u(x,t)$ = cell velocity

$\alpha$ = transition rate from $P$ to $Q$

$\beta$ = transition rate from $Q$ to $P$

$\gamma$ = growth rate of cycling cells

$\lambda_{P}$ = proliferative cell loss rate

$\lambda_{Q}$ = quiescent cell loss rate

$\eta(t)$ = transition rate from $Q$ to $P$ induced by inflammation

$\mu_{P}(x,t)$ = effectivity of the drug against $P$-cells

$\mu_{Q}(x,t)$ = effectivity of the drug against $Q$-cells

$E$ = volume fraction of extracellular fluid

Equation (1) describes the volume fraction of the arterial wall occupied by cells in the proliferative phase. The left-hand side is a material derivative, as we are concerned with describing the time rate of change of material elements of cellular tissue which are subjected to a velocity field which varies in both space and time. As a consequence of this, the cell velocity is a function of both space and time. The right-hand side describes an increase in the number of proliferative cells by cell division at a rate $\gamma$, and by quiescent cells transitioning to proliferative at a rate $\beta$ which is enhanced by $\eta$ (the inflammation rate due to stent implantation) and reduced by $\mu_Q$ due to eluted drug blocking the transition. The proliferative fraction also decreases due to transitions from the proliferative to quiescent phase at a rate $\alpha$ and due to proliferative cell death at a natural rate of $\lambda_P$ and a drug induced rate of $\mu_P$.

Similarly, equation (2) describes the rate of change of the volume fraction of quiescent cells. The left hand side is again a material derivative. The right hand side describes an increase in the number of quiescent cells due to proliferative cells transitioning to quiescence at a rate $\alpha$ as well as a decrease in the number of quiescent cells due to cell death at a rate $\lambda_Q$ and transitioning to the proliferative state at a rate $\beta+\eta-\mu_Q$, describing a base rate enhanced by inflammation and reduced by the action of drugs.

The final equation describes the relationship between cells in the various proliferative phases and the extra-cellular fluid surrounding the cells. We may consider the volume fraction of cells to be $V = 1 - E$, noting that the volume fraction implicitly describes the radial velocity, $u$, of cells throughout the arterial wall. The radial velocity arises from the accumulation of new cells and is used to describe the model the change in thickness of the arterial wall.

We take $x$ to be a spatial coordinate through the intimal layer, with $x = 0$ corresponding to outermost extent of the intima which is in contact with the lumen and $x = L$ corresponding to the innermost part. For simplicity we assume that the intimal layer is thin and neglect curvature correction to the advection terms, following [15]. The system of equations requires one boundary condition on $u$ which we take to be $u(x = 0,t) = 0$. As the inner boundary remains fixed by the stent, the size of the domain, i.e. the thickness of the vessel wall, satisfies a Stefan condition in velocity, i.e.

$\dfrac{\text{d}L}{\text{d}t} = u(L).$

(4)

To solve these equations we take initial conditions for $P = P_0$, $Q = Q_0$, $u = 0$, $L = L_0$ at $t = 0$.

It should be noted that we restrict our model to the initial post-implantation period, i.e. that before the contraction of the vascular thickness. Noting that it is the early post-implantation phase which has the largest implications for patient care and the dynamics of cell proliferation at this stage are dominated by the inflammatory response and pharmaceutical effects, we thus consider it reasonable to investigate only this phase of a stent life-cycle.

---

3. Parameter estimation

The difficulties presented by this model lie not in its solution but rather in the fixing of the parameters. It is hoped that, by elucidating a path through which these parameters may be fixed from experimental data, we may guide experimentalists to provide exactly this type of data. As it stands now, the fixing of parameters is a non-trivial task as there exists a paucity of experimental data of the requisite type. We may first simplify the above model in order to reduce the number of undetermined constants which must be fixed.

In order to simplify the notation we notice that the growth and death rates of the cycling cells may be combined linearly, as we presume that the death rate is in general smaller than the growth rate. We thus combine them, forming the effective growth rate, $\gamma'$:

$\gamma' = \gamma - \lambda_{P}$

(5)

We further linearly combine the rates governing the $Q \to P$ transition, such that $\psi$ is the net transition rate:

$\psi = \beta + \eta$

(6)

Thus, we have reduced our system to:

$\frac{\partial P}{\partial t} + \frac{\partial}{\partial x}(uP) = (\gamma' - \alpha - \mu_{P})P + (\psi - \mu_{Q})Q \label{eq:reducedP_model}$

(7)

$\frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}(uQ) = \alpha P - (\psi + \lambda_{Q} - \mu_{Q})Q \label{eq:reducedQ_model}$

(8)

The challenge now lies in determining the behaviour of the drug effectiveness rates, $\mu_{P}$ and $\mu_{Q}$. These two constants sufficient to describe the behaviour of a drug which acts to control inflammation. This implies that we may determine the other constants through experiments on bare-metal stents, i.e. stents where no drug is present. This control-group approach has in practice been taken in experiments (cf. [7], [23]). Our proposed method of parameter-fixing then begins by considering the case where $\mu_{P} = \mu_{Q} = 0$.

In this section we consider two slightly different steady states. The first steady state corresponds to the uninflamed case. This corresponds to a steady state in which $P$ and $Q$ are constants and the velocity field $u$ is zero (technically we only require $\tfrac{\partial u}{\partial x} = 0$). In this case the thickness of the artery is not changing. The second steady state we discuss occurs in the inflamed case. In this steady state $P$ and $Q$ are constants, however the thickness of the artery is increasing so we do not have a constraint on $u$.

Firstly, we consider the steady-state case where inflammation is not present. As there must be a steady-state proliferative fraction and thickness towards which the vasculature tends after the implantation of a stent, the implanted device represents a perturbation from a steady-state. At this steady-state, the rate of the $Q \to P$ transition is very small, i.e. $\psi \sim 0$ as $t \to \infty$. This is a necessary physical condition to maintain a constant thickness of the artery in the absence of any external inflammatory or pharmaceutical effects. By summing equations (7 - 8) and assuming constant porosity, i.e. that $E$ is a constant, we obtain:

$V \frac{\partial u}{\partial x} = \gamma' P - \lambda_{Q} Q \label{eq:dudx}$

(9)

Where $V = 1-E = P+Q$ and, by the constant porosity assumption is a constant. As the vessel wall is assumed to have constant thickness at steady state, $\frac{\partial u}{\partial x} = 0$, and we find that there exists a steady-state $P$ as well as $L$, i.e.:

$P_{eq} = \frac{\lambda_{Q} V}{\gamma' + \lambda_{Q}}.$

(10)

We may assume that during the inflammatory phase, $\lambda_{Q} \ll \psi$, i.e. the dynamics are dominated by the transition to proliferation as opposed to necrosis. This is a necessary condition for an increase in the thickness of the vessel wall. It may also be shown by considering the volume fraction of cells that when $\psi$ is relevant and $\lambda_{Q}$ is assumed to be negligible, the rate of change of $P$ is governed by the following ordinary differential equation:

$\frac{\partial P}{\partial t} = -\frac{\gamma'}{V}P^{2} + (\gamma' - \alpha - \psi)P + \psi V \label{eq:dpdt_para}$

(11)

This equation is derived by taking equation 7, using the equation $V = P+Q$ to eliminate $Q$, using equation 9 to eliminate the spatial derivative in $u$ and, as discussed above, setting $\lambda_{Q} = \mu_{P} = \mu_{Q} = 0$. Finally we assume that $P$ is uniform across the thickness of the artery wall so that its spatial derivative is zero.

On the other hand, in the case with no inflammation, i.e. prior to implantation of the stent or after significant time has elapsed, $\psi$ is effectively zero and equation (11) reduces to:

$\frac{\partial P}{\partial t} = -\frac{\gamma'}{V}P^{2} - (\gamma' - \alpha)P \label{eq:dpdt_normal}$

(12)

We see the two possible graphs in Figure 3, where the higher of the two parabolae is the inflamed case (eq. 11) and the lower is the uninflamed case (eq. 12). Thus, there exist two steady state values (where $\frac{\partial P}{\partial t} = 0$) for the uninflamed case: one when $P = 0$ and the other when it is very small but positive, as previously discussed (i.e. $P_{eq,healthy}$). In order to enforce this, $\gamma'$ and $\alpha$ should be very similar, although $\alpha$ must be slightly larger. If these two constants are the same, the expression for steady state behaviour in the absence of external forcing reduces to a simple parabola, with a double root at $P = 0$.

The values of $P_{eq,\mathit{infl}}$ and $V$ may be estimated from experimental data. Using this, we may estimate the value of $\psi$ during the inflammatory phase through the distance between the two solutions, indicated with a double-headed arrow in Figure 3. As labelled, the distance here is $\psi(V-P)$. For known values of $V$ and $P$, $\psi$ may therefore be estimated.

From conservation of volume, it follows that:

$\frac{\partial u}{\partial x} = \frac{d L}{d t}\frac{1}{L}$

(13)

Recalling equation (9) and that we have shown immediately above that $\lambda_{Q} \ll \psi$ in the inflammatory limit allows us to write:

$\gamma' = \frac{V}{LP}\frac{dL}{dt} \label{eq:gamma_prime}$

(14)

We then assume that $\frac{dL}{dt}$ is constant on the initial, highly inflamed response (cf.Figure 1). Using a fit through the data in this region to approximate $\frac{dL}{dt}$, $\gamma'$ may be determined as well. As a first estimate, it was assumed that $\alpha$ takes the same value, but it was found in practice that this leads to significant issues with simulation. To explain this, consider equation (14), which may be re-written as:

$\frac{dL}{dt} = \gamma' \frac{LP}{V}$

(15)

We see that the rate at which the thickness, $L$, of the layer changes is proportional to both the proliferative fraction, $P$, and the thickness of the cell layer. If $P$ does not go to zero quickly enough, these terms balance each other and cause $\frac{dL}{dt}$ to achieve a non-zero, positive steady state, resulting in unchecked growth of the artery. A value of $\alpha$ slightly larger than that of $\gamma'$ is necessary to prevent this. In practice, the value of $\alpha$ may be determined from equation (12).

---

3.1. Inclusion of drug effects

The equations governing the growth of the cell wall in response to the implantation of a stent in the absence of anti-proliferative agent may be broken down onto two temporal domains. On the first domain which follows the implantation of the sten, inflammatory effects are relevant, i.e. $\eta > 0$ This directly implies that the $P \to Q$ transition rate is relevant. After some time has passed the dynamics of the problem change, and on the second domain, the inflammatory response has died away, i.e. $\eta = 0$. Here, the dynamics of the system are dominated by the return of the P-fraction to its normal value. We shall assume for this section that $\eta$ takes some constant value on this first domain, and therefore so does $\psi$.

We begin by considering the equation on the first domain, where the inflammatory $P \to Q$ response is present, i.e. $\psi \neq 0$. In this case, equation (11) may be rewritten in terms of its two roots, denoted $P_{+}$ and $P_{-}$, with $P_{+}$ being the positive root. This yields:

$\frac{\partial P}{\partial t} = -\frac{\gamma'}{V}(P-P_{+})(P+P_{-})$

(16)

We now rescale in the region where $P \sim P_{+}$, yielding:

$\frac{\partial \bar{P}}{\partial \bar{t}} = -(\bar{P}-1)(\bar{P}+\rho) \label{eq:rescaled_infP}$

(17)

where:

$\begin{align} P & = P_{+}\bar{P} \\ t & = \dfrac{V}{\gamma' P_{+}}\bar{t} \nonumber \end{align}$

(18)

and

$\rho = \frac{P_{+}}{P_{-}}$

(19)

We proceed by expanding the right-hand side of equation (17) in a Taylor series about $\bar{P} = 1$. To leading order, we have:

$\frac{\partial \bar{P}}{\partial \bar{t}} = -(\rho + 1)(\bar{P}-1) + \mathcal{O}(\bar{P}^{2})$

(20)

that is:

$\begin{array}{l r} \bar{P}(\bar{t}) \sim 1 + (\bar{P}_{0}-1)e^{-(\rho + 1)\bar{t}} \ \ {\text{as }}\ \bar{P} \to 1 \end{array}$

(21)

We require an initial value of $P$, i.e. $P_{0}$, in order to find the particular solution. Recall that this equation describes the initial response of the vasculature to the implantation of the stent, thus at $\bar{t} = 0$, $P$ and $Q$ take their 'normal', i.e. unstimulated values. We have already discussed that, in the absence of stimulus, $P$ takes a very small positive value. We shall here assume that this value is small enough to set $P_{0} = 0$ and, upon applying this initial condition, we obtain the particular solution describing the fraction of P-cells on the first domain, denoted $\bar{P}_{1}$, i.e. that where $\psi$ is relevant.

$\bar{P}_{1}(\bar{t}) = 1 - e^{-(\rho + 1)\bar{t}} \label{eq:P1_tau}$

(22)

On the second temporal domain, we consider the case where $\psi = 0$, resulting in the following, simplified form of equation (11):

$\frac{\partial P}{\partial t} = P \left(-\frac{\gamma'}{V}P - \delta \right)$

(23)

where $\delta$ is defined as:

$\delta = \alpha - \gamma'$

(24)

In this case, we rescale as follows:

$\begin{align} P & = \dfrac{\delta V}{\gamma'} \bar{P} \\ t & = \dfrac{1}{\delta}\bar{t} \nonumber \end{align}$

(25)

resulting in the rescaled version of the equation on this domain:

$\frac{\partial \bar{P}}{\partial \bar{t}} = -\bar{P}(\bar{P} + 1) = -\bar{P} + \mathcal{O}(\bar{P}^{2})$

(26)

working to the same order of accuracy as above. The general solution on the second domain, $\bar{P}_{2}$, is then:

$\bar{P}_{2}(\bar{t}) = P_{m}e^{-t} \label{eq:P2}$

(27)

where $P_{m}$ acts as the initial condition on this domain, and is the maximum value achieved on the first domain as seen in Figure 4. This value is the steady state to which equation (22) tends as $\bar{t} \to \infty$. In effect, the proliferative fraction on the first domain very quickly reaches a steady state value where is remains until $\psi$ becomes zero. At this point, $P$ begins to decline as above, tending back to its original value, as shown in Figure 4, where the solution on the first time domain is shown in red and on the second in blue.

It is notable that after the inflammatory response (in the sense of a positive value for $\psi$) is present, inflammation is still observed in the form of residual $P$-cells. The $P$-fraction takes several weeks to return to its healthy steady-state value, leading to difficulties in describing exactly the nature of the inflammatory response.

Of particular interest is the thickness of the arterial wall, or perhaps some portion thereof. From conservation of volume it follows that we may describe this thickness through a Stefan condition. We have re-written this condition such that it depends on the proliferative fraction directly instead of the cell velocity at the boundary (eq. 14). We may now write it in dimensionless form as:

$\frac{dL}{dt} = \zeta L(t) P(t) \label{eq:dldt}$

(28)

where the constant, $\zeta$, depends on the relevant time- and P-scales. Noting that we have dropped the use of overbars for rescaled $P$ and $t$, and are using the convention of letting square braces denote a scale, $\zeta$ is:

$\zeta = \frac{\gamma'[t][P]}{V}$

(29)

The general solution to equation (28) is:

$L(t) = L_{i}e^{\zeta \int_{0}^{t}P(\tau)d\tau}$

(30)

It is then a simple matter to show that the change in length with respect to time is defined as follows on the first time domain:

$L_{1}(t) = L_{0}\exp\left ({\zeta \frac{(1 + \rho)t + e^{-(\rho + 1)t} - 1}{1 - \rho}}\right ) \label{eq:L1}$

(31)

where $L_{0}$ is the initial thickness before implantation; and as follows on the second domain:

$L_{2}(t) = L_{m}\exp\left({\zeta(\ln(e^{t}(P_{m} + 1) - P_{m}) - t)}\right)$

(32)

where $L_{m}$ is the thickness at the end of the inflammatory phase. The above may be simplified to yield:

$L_{2}(t) = L_{m}(P_{m}(1 - e^{-t}) + 1)^{\zeta} \label{eq:L2}$

(33)

Thus, we see that the presence of the $\psi$ response causes the thickness of the vessel wall to increase without limit. Considering the second domain, however, it is clear that the subsequent decline in the proliferative fraction to its undisturbed value causes the thickness to tend towards some steady limit, i.e.:

$L_{inf} = L_{m}(P_{m} + 1)^{\zeta} \label{eq:L_inf}$

(34)

As with the proliferative fraction, the two solutions are shown together in Figure 5. The reader will note the presence of an inflection point between the two domains as the thickening response changes nature in keeping with the inception of decline in the $P$-fraction.

We are ultimately concerned with fitting a descriptive model for the vascular response after the implantation of a stent which considers the effects of drug elution from the stent into the vasculature. The problem is now to determine the drug effectiveness constants $\mu_{P}$ and $\mu_{Q}$. In practice, the determination of these constants is largely through a comparison with the analytical results above and requires that quality experimental data of both the bare-metal and drug eluting cases exists.

We begin by assuming that the rate constants which have been determined earlier still hold, i.e. that the effect of the drug is to suppress the un-medicated proliferative behaviour of the vasculature. Thus, we are able to compartmentalise our earlier work. We use the drug delivery model to describe the concentration of drug present in the intima, noting that we consider herein only the case with binding to the cells and restricted diffusion in the stent coating. The task of coupling the models then falls to including the effects of the drug as a function of concentration in the cell proliferation model.

---

3.2. Consideration of drug effectivity

We now return to the drug effectivity constants, $\mu_{P}$ and $\mu_{Q}$. The effects of the drug have been considered by some authors (e.g. [16]) to be proportional to the concentration of the drug. See [20] for a different perspective on this in which the efficacy of the drug is taken to depend on stair step regions of saturation. In order to enable analytical results, we shall assume a smooth mapping. For example, the effectiveness of the drug may be defined as:

$\mu = Ce^{-\frac{1-a_{t}}{\tau}} \label{eq:mu_P}$

(35)

where $C$ and $\tau$ are constants of proportionality and $a_{t}$ is the total drug concentration, scaled with the maximum possible concentration which depends on the amount of drug present in the stent and the elution rate such that $0 \leq a_{t} \leq 1$, shown in Figure 6, following [2].

Of analytical difficulty is the fact that effectiveness of the drug is taken as a function of the total drug concentration in the intima. This concentration is, in turn, a function of both position and time. According to [21], the mass transfer Fourier number is very large when considering the intima, and so we make the approximation that the concentration of drug is approximately independent of position, and thus somewhat reduce the complexity of the problem. We also note from numerical studies that the concentration of drug appears to vary slowly in time as well as position in the intima, due largely to the restricted diffusion out of the stent coating [21]. Thus, in order to make analytic progress, we shall assume that the concentration of drug in the intima may, for small timescales, be assumed to be constant.

Two of the main drugs with which drug eluting stents are loaded are Paclitaxel and Sirolimus (Rapamycin) with many of the others, such as Zotarolimus and Everolimus being analogues [10]. The effect of these drugs is to somehow inhibit the proliferation of cells in order to prevent restenosis, although it should be noted that at higher concentrations they may induce cell death, as in the case of Paclitaxel when used for cancer treatment [24].

The mode of action of Paclitaxel is to polymerise the $\alpha$- and $\beta$-units of tubulin, this preventing cells in the $G2$ phase from entering the $M$ phase. Sirolimus is a macrocyclic antibiotic with potent immunosuppressive properties, which works by binding to specific cytosolic proteins and blocking cell proliferation. It is further a strong inhibitor of inflammation and is not toxic to cells in low doses.

---

3.3. A growth-inhibiting model

It then makes sense to consider some general anti-proliferative agent as acting exclusively on $P$-cells as a first approximation, following [16]. Thus, a model of drug action which acts on the cells in the proliferative phase only was considered. Drug action was considered to be a function of the total (i.e. free plus bound) drug concentration in the intima. Consider our general system for the cell proliferation with the inclusion of the drug effects, $\mu_{P}$:

$\frac{\partial P}{\partial t} + \frac{\partial}{\partial x}(uP) = (\gamma' - \alpha - \mu_{P})P + \psi Q \label{eq:dpdt_mup}$

(36)

$\frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}(uQ) = \alpha P - \psi Q$

(37)

Referring to equation (36), it may be seen that the effect of the drug, $\mu_{P}$, is to inhibit the net growth of the vasculature by reducing the effective value of $\gamma'$, thus reducing the growth rate of $P$-cells. Consider that we may group terms as per:

$\frac{\partial P}{\partial t} + \frac{\partial}{\partial x}(uP) = [(\gamma' - \mu_{P}) - \alpha ]P + \psi Q$

(38)

This should hopefully make it clear that the effect of $\mu_{P}$ may change character somewhat depending on its magnitude. For $\mu_{P} < \gamma'$, we see it has the effect of reducing the growth rate observed in the intima. However, when $\mu_{P} > \gamma'$, the growth rate of cells becomes negative. Effectively, the drug at this concentration appears to be killing cells off, rather than simply inhibiting their growth.

It has been suggested by [18] that there exists some therapeutic range of delivered concentration for anti-proliferative agents in drug-eluting stents, below which the drugs are ineffective and above which they are toxic. The point at which the growth rate of the cells becomes negative may be understood to be the upper bound of the therapeutic range. We shall see that this point has some interesting mathematical ramifications for our model as well.

The intention is to describe which values are permissible or plausible for $\mu_{P}$ to take. Thus, we assume that $\mu_{P}$ may be approximated as constant at discrete points in time and that we may then approximate the behaviour of the system in response to this. The analysis may be extended to the more general case where $\mu_{P}$ is not constant, nor is $\psi$.

We may write the change in the $P$-fraction with respect to time in a similar fashion to equation (11), but we include the effects of the anti-proliferative drug here, yielding:

$\frac{\partial P}{\partial t} = -\left (\frac{\gamma' - \mu_{P}}{V}\right ) P^{2} + (\gamma' - \alpha - \psi - \mu_{P})P + \psi V \label{eq:dpdt_kill}$

(39)

The origin of the quadratic term is in equation (9), which describes $u$ in terms of $P$ and permits us to eliminate $u$ from the second term on the left-hand side. Importantly, we have that $\frac{\partial u}{\partial x} \propto P$ and so a positive prefactor on the quadratic term (once it has been moved to the right-hand side) would imply a contraction of the cellular region. This is something we wish to avoid as in this model such a scenario has no stable, real solutions, as can be seen in the case of the bold curve shown in Figure 7.

We note that the right-hand side considered as a function of $P$ forms a parabola which opens downwards. In order to agree with experimental results, we require that $\frac{\partial P}{\partial t}$ remains positive (i.e. the proliferative fraction and thickness continue to increase) during the initial, $\psi$-positive phase.

We thus require that in order for physically valid solutions to exist:

$\mu_{P} \leq \gamma' \label{eq:muP_phys_valid}$

(40)

We now seek the value of the positive steady state, recalling that the drug-free model had been developed in order to yield a steady state in the first domain at $P = 0.96V$. From experiments [7] we seek a steady state at $P = 0.29V$. Application of the quadratic formula to equation (39) yields a positive root at:

$P_{+} = \frac{-(\psi + \mu_{P})+\sqrt{(\psi + \mu_{P})^{2} + 4\psi \left (\dfrac{\gamma'-\mu_{P}}{V} \right )}}{2(\gamma'-\mu_{P})/V} \label{eq:Proot_kill}$

(41)

We expect that $\mu_{P}$ must be close in size to $\gamma'$, so we define:

$\gamma' - \mu_{P} = \epsilon \ll 1$

(42)

We exploit the smallness of $\epsilon$ and perform a binomial expansion of equation (41), yielding:

$P_{+} = \frac{\psi}{\psi + \mu_{P}} + \mathcal{O}(\epsilon^{2})$

(43)

Using the estimation procedure outlined in Fig. 3 applied to data from Ref. [7] we can estimate $\psi \approx 23\gamma'$. In addition, as discussed in Section 3 $\mu_{P} \leq \gamma'$. Thus,

$P_{+} \gtrsim \frac{23\gamma'}{(23 + 1)\gamma'} = \frac{23}{24} \label{eq:p_plus_crit}$

(44)

Effectively, the drug may suppress the growth rate enough to keep $P$ constant without inducing a significant increase in the intimal thickness. However, the existence of the critical value in equation (44) justifies the physical requirement given in equation (40), i.e. a growth-inhibiting drug cannot go beyond a certain point of effectivity without leading to a decline in the wall thickness.

---

3.4. A transition-blocking model

We now turn our attention to an alternate model of drug action, wherein the drug is taken to be effective against $Q$-cells, namely in preventing their transition to the proliferative phase. This is of mathematical interest as we have not found any result which would prohibit large-scale inhibition of the rate constant, $\psi$.

A physical justification for this model may be presented as well. Recall that the individual phases of the overall proliferative phase were not considered in detail. Rather, the entire phase was considered together. However, the effect of some drugs is to block the transition into the mitotic portion of the proliferative cycle, wherein cell division occurs. Blocking the $Q \to P$ transition is effectively the same, from the perspective of this model, as blocking the $G1 \to M$ transition.

Consider the modified system of equations where we consider the drug to act against $Q$-cells through transition blocking rather than against $P$-cells through growth inhibition, hence the consideration of $\mu_{Q}$ instead of $\mu_{P}$:

$\frac{\partial P}{\partial x} + \frac{\partial}{\partial t}(uP) = (\gamma' - \alpha)P + (\psi - \mu_{Q})Q$

(45)

$\frac{\partial Q}{\partial x} + \frac{\partial}{\partial t}(uQ) = \alpha P - (\psi - \mu_{Q})Q$

(46)

where $\mu_{Q}$ is defined in the same fashion as $\mu_{P}$, i.e.:

$\mu_{Q} = \psi_{0}e^{-\frac{1 - a_{t}}{\tau_{Q}}} \label{eq:mu_q}$

(47)

where $\tau_{Q}$ is a constant and $\psi_{0}$ is the activated value of $\psi$ such that $0 < \mu_{Q} < \psi$. This is a necessary condition, as values greater than $\psi$ no longer represent blocking of the $Q \to P$ transition, but rather lead to a boundless proliferation of $Q$-cells.

Considering equation (47), we may redefine $\psi$ as follows, such that it considers the effect of the drug:

$\psi' = \psi e^{-\frac{a_{t}}{\tau_{Q}}} \label{eq:psi_prime}$

(48)

Equation (48) clearly takes a value of $\psi_{0}$ during the inflammatory phase if no drug is present, and of 0 after the inflammatory phase, as seen in Figure 8.

Thus, we may write the change in the proliferative fraction, following equation (11) as:

$\frac{\partial P}{\partial t} = -\frac{\gamma'}{V}P^{2} + (\gamma' - \alpha - \psi')P + \psi' V$

(49)

In this case, we may easily compute the positive steady state value of $P$ to be:

$P_{+} = \frac{-\psi'+\sqrt{\psi'^{2} + 4\psi' \left (\dfrac{\gamma'}{V} \right )}}{2\gamma'/V} \label{eq:Proot_block}$

(50)

We note that all rate constants remain non-negative at all times, and thus solutions which yield a non-negative steady state exist for all allowable values of $\psi'$. Thus, it is reasonable to use a $\mu_{Q}$-model to induce the requisite $P$-fraction at some given time. Again, this result agrees with what one would expect as the presence of the blocking drug does not kill off cells with an associated possibility of negative $\frac{dL}{dt}$, but rather slows the $Q \to P$ transition leading to - at its most severe - no change in the intimal thickness.

We have shown here a method by which the rates both $\gamma'$ and $\alpha$, i.e. the effective growth rate and the $P \to Q$ transition rate, may be determined from experiments involving bare metal stents. We have further outlined a method by which the effective $Q \to P$ transition rate, $\psi'$, as well as the effectivities of the drug on both the quiescent and proliferative phases may be determined, albeit under fairly strong assumptions about their behaviour. In order for this method to be put into practice, it is required that experimental data be developed which yields a detailed and reliable view of both the proliferative fraction and the thickness at various times throughout the course of treatment, both with and without the presence of some anti-proliferative agent.

---

4. Discussion

The cell proliferation model was based on a consideration of three phases: quiescent and proliferative cells and extracellular fluid. Constant porosity was assumed, eliminating the need for detailed consideration of the extracellular fluid. It was also assumed that the four sub-phases of the proliferative phase could reasonably be combined. The drug was taken to act independently on the proliferative and quiescent phases and to be a function of the concentration at any point in the intima. Another important assumption was that the various rates associated with the cell proliferation combined linearly, i.e. the effects of the drug could be added directly to the underlying rates of $Q \to P$ transition, cell mitosis, etc.

We have shown that it is possible to obtain results with this model which appear to correspond to the experimental results over the initial inflammatory response period. There does not exist, however, sufficient experimental data to fix the parameters and perform a dimensional simulation of the response. A method to estimate the various rate constants from easily-obtained experimental data is presented herein in the hopes that it motivates further experimental research and drives development of this field of research.

With regards to the experimental results, the proliferative fraction and thickness at a reasonable temporal resolution should be considered the bare minimum experimental requirement for the determination of the parameters of this model. Should it be possible to determine any of the rate constants in equations (1 - 2) experimentally, this would certainly be the preferred option. Of particular importance is the determination of the $Q \to P$ transition rate, $\psi$, and the effectivities of the anti-proliferative and anti-inflammatory agents $\mu_{P}$ and $\mu_{Q}$.

Also notably lacking from experimental results is the response in the period immediately following implantation. A simple stepped function was assumed for this response, $\psi'$, in the interest of making analytical progress, but it is reasonably unlikely that this assumption is accurate in practice. It is suspected, rather, that it takes the form of some decaying function of time. The method outlined above, while suitable to determine the other constants is insufficient to determine the behaviour of this rate. Progress may likely be made on the determination of $\psi'$ by solving the inverse problem, but this will again require a reasonable amount of experimental data, in particular for several points in time.

With improved experimental data, it is hoped that significant progress may be made on this topic. It is already apparent that this model wields significant predictive power. The ultimate outcome of this area of research should be in greatly improved predictive capabilities with respect to the biophysical and biochemical responses of the body to a drug-eluting stent. Such information would, it is hoped, lead to an ability to develop designer stents or, at the very least, improve currents designs with an eye towards reducing the rates of restenosis which are currently observed.

---

5. Conclusion

Drug eluting stents have been very successful in preventing restenosis following treatment of atherosclerosis by percutaneous balloon angioplasty. However there is scope for improvement, and a holy grail of the research community has been to use modelling to optimise the design of these stents to maximise their effectiveness. To use modelling to predict whether a new design of stent will prevent restenosis it is essential that the model can describe the processes leading to restenosis. Research in this area has focussed on the mechanical and chemical aspects of the problem: calculating fluid mechanical stresses on the artery, and modelling the transport of drugs through the artery wall. To complete the picture a model of the biology of restenosis is also required. We have shown that mathematical models of cancer biology can be adapted to describe the cellular processes of restenosis. Future work should focus on collecting experimental data which can be used to refine the model, and incorporating mechanical effects into the model.

---

Acknowledgments

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

---

Appendix A. Tables of equations

---