## Introduction

Arterial stenting, acting as a supporting scaffold, has revolutionized the treatment of coronary artery disease, as it reopens the occluded vessel and restores the normal flow of blood. The bare metal stents (*BMS ^{s}*), while revolutionary at the time, were soon rendered unsatisfactory due to their inability to prevent in-stent restenosis. The next wave of arterial stents coated by drug–the so-called drug-eluting stents (

*DES*) raised all sorts of questions by releasing antiproliferative agents in a controlled manner into the injured site to reduce restenosis rates.

^{s}^{1–9}Drug-eluting stents are now the primary choice of percutaneous coronary interventions (

*PCI*) in millions of patients, but questions regarding their longevity and safety still arise.

^{s}Some experimental studies have been carried out in the recent past, with the aim of quantifying the capability of this device to reduce in-stent restenosis after stent implantation.^{10–12} Lovich *et al.*^{10} studied the behaviour of heparin in explanted arteries and concluded that the presence of binding sites changes along the transmural direction, being higher in the endothelium and lower in the adventitia. Lovich and Edelman studied the effects of specific binding sites inside the arterial wall on the drug uptake,^{13} where the presence of specific binding site action was modelled using the reversible chemical reaction. Sakharov *et al.*^{14} disregarded the convective effects on the transport of free drug. Hwang *et al.*^{15} predicted the free as well as bound drug concentrations by solving for distribution of free drug and then using a multiplicative factor (partition approach) to predict the concentration of bound drug. Migliavacca *et al.*^{16} studied the drug release pattern in vascular wall from drug-eluting stents using a single species approach along with a partition coefficient approach to relate the free and the bound drug concentrations. Borghi *et al.*^{17} opined that the inclusion of reversible binding leads to delayed release and that the erosion of polymer affects the drug release from a single strut. Horner *et al.*^{18} considered a two-species drug delivery model including reversible binding sites, and their model predicted that a single species drug delivery model cannot accurately predict the distribution of bound drug. They also concluded that a two-species approach that includes reversible binding is the way forward for future stent-based drug delivery systems.

The main objective of this study was to advance the aforementioned work with a two-species model of drug transport eluted from a pair of struts, where the transport of free drug is governed by reaction-diffusion process and that of bound drug, assuming complete immobilization in the tissue, by reaction process. Following Tzafriri *et al.*,^{19} a second-order dynamic model that describes a saturating reversible binding process by treating bound drug as a dynamic variable has been taken into account to explore drug interaction with cells of the arterial wall. In most of the studies cited above, transient drug release has been modelled as a uniform release, which is unrealistic and not representative of actual stent-based delivery. Instead, a simple time-dependent Dirichlet boundary condition is applied on the surface of the struts.^{20–24} Arterial properties, such as porosity and tortuosity, dictate the transport of drugs within the arterial tissue.

When an endovascular drug-eluting stent is implanted, it has major impact on the structure of the arterial wall, eventually influencing the overall rates of diffusion through tissues.^{25} For diffusion in a porous material, the effective diffusion coefficient is assumed to depend on two factors viz porosity (a dimensionless parameter, which is the ratio of pore volume to the total material volume) and diffusion path tortuosity (ratio of the actual pore length to the distance between its ends; *i.e.* arc-chord ratio)^{26}—these parameters change the free diffusivity of the drug eluted from a pair of struts.^{27} Thus, the effects of porosity and tortuosity on the diffusivity of drug have also been accounted for in this present investigation.

## Governing equations and boundary conditions

The freely-transported drug, referred to as ‘free drug’, is denoted by *c _{f}* and the drug, completely immobilized in the tissue, is referred to as ‘bound drug’, which is denoted by

*c*. The interconversion of drug between the unbound plasma phase and the bound phase of tissue binding sites is controlled by a second-order chemical reaction. The computational domain is comprised of a long axial section of length

_{b}*L*and the wall thickness is taken to be 10 times the strut height (δ). The axis of symmetry is taken along the centerline of the artery (Fig. 1). The transport of free drug eluted from a pair of struts is governed by unsteady reaction-diffusion process (Eq. 1) and that of the bound drug is represented by unsteady reaction process (Eq. 2).

Symmetry boundary conditions for both the free and the bound drugs are applied at the proximal (Γ* _{ti}*) and the distal (Γ

*) walls (Eq. 3).*

_{to}^{28,29}Impermeable boundary condition for bound drug is assumed at the perivascular wall (Γ

*), lumen-tissue (Γ*

_{tp}*) and strut-tissue (Γ*

_{bt}*) interfaces (Eq. 4). For the free drug, perfect sink condition is imposed at the perivascular end (Eq. 5). Since a proper boundary condition for the free drug at lumen-tissue interface (Γ*

_{st}*) is not readily apparent, we considered two opposing extremes, either that flowing blood is extremely efficient at washing out mural-adhered drug, modelled as a zero-concentration interface condition, or that mural-adhered drug is insensitive to flowing blood, modelled as zero-flux boundary condition (Eq. 6).*

_{bt}^{30–31}Instead of modelling uniform release of drug from a pair of struts, we assumed simple time-dependent release kinetics with a flux condition (Eq. 7).

^{13,20}Therefore, the governing equations representing the transport of free as well as bound drugs together with their appropriate boundary conditions are as follows:

*z*and

*r*coordinates are along the axial and the radial directions respectively.

*k*and

_{a}*k*are the rates of association and dissociation constants respectively.

_{d}*B*is the net tissue binding capacity.

_{M}*D*and

_{c}*c*

_{0}are the diffusion coefficient and initial concentration in coating respectively.

Here, *D _{T}*, the true diffusivity of the free drug can be written as:

^{27,32}

*D*and

_{free}*D*are the coefficients of free and effective diffusivity respectively.

_{eff}*R*(=

_{d}*k*/

_{d}*k*) is the equilibrium dissociation constant.

_{a}We now introduce nondimensional variables as

*V*is the transmural filtration velocity.

_{y}Under these assumptions, the equations (1–7) take their respective nondimensional forms (dropping tilde) as

*Pe*[= (

_{T}*V*δ)/(

_{y}*D*)] and

_{T}*Da*[= (

*k*δ

_{a}B_{M}^{2})/(

*D*)] are the Peclet number and Damköhler number in the tissue respectively. Here, ε

_{T}_{1}(=

*R*/

_{d}*c*

_{0}) and ε

_{2}(=

*c*

_{0}/

*B*) are two scaling parameters.

_{M}*Pe*{= [

_{c}*V*(

_{y}*h*

^{2}/δ)]/

*D*} is the Peclet number in the coating of struts and

_{c}*h*is the thickness of the coating of the strut.

## Solution procedure

The governing equations (8–9) representing the transport of free and bound drugs along with the boundary conditions (10–14) are solved numerically by finite-difference scheme. Forward-time centered-space discretization technique is leveraged in this explicit scheme. Following, we describe our finite-difference scheme in more detail. We denote *x _{j}* =

*j*δ

*x*,

*z*=

_{i}*i*δ

*z*,

*t*=

^{n}*n*δ

*t*, where

*n*refers to the time level and δ

*t*is the time increment. Here, δ

*x*and δ

*z*stand for the step sizes along the radial and the axial directions respectively. Steady state is achieved when the convergence criterion for concentration was 10

^{−7}for both drug forms.

The finite-difference approximation of (8) for the transport of free drug is given by

Likewise, the discretized version of Eq. 9 may be written as

## Result and discussion

For the purpose of numerical computation of the desired quantities of major physiological significance by using the baseline values included in Table 1^{7,13,23,28,33–41}, solutions are computed through the generation of grids, with a size of 301 × 101 for δ*t* = 0.00001. The simulation concerning the grid independence study was performed for the purpose of examining the error associated with the grid sizes used, and is presented in Figure 2a. One may notice from this figure that the transmural variation of normalized free drug concentration concerning the three grid sizes viz 151 × 51, 301 × 101, 601 × 201 overlap one another at *t* = 100 for *Pe _{T}* = 2 and

*Da*= 40. Thus, the grid independence study in the present context of numerical simulation has its own importance to establish the correctness of the results obtained. Figure 2b displays the results of normalized free drug concentration for different time steps δ

*t*viz 0.0001, 0.00001, 0.00005. Here too, the concentration profiles get overlapped for three distinct time steps, which also establish the correctness of the time step used.

Description | Parameter | Value | Reference |
---|---|---|---|

Strut dimension, m | δ | 0.0001 | ^{[13]} |

Strut coating thickness, m | h | 5.0 × 10^{−5} | ^{[33]} |

Mean wall thickness, m | Aw (= 10δ) | 0.001 | ^{[34]} |

Interstrut distance, m | Δ (= 7δ) | 0.0007 | ^{[35]} |

Transmural filtration velocity, m s^{−1} | V_{y} | 4 × 10^{−8} | ^{[36]} |

Porosity of the arterial wall | ε | 0.787 | ^{[37]} |

Tortuosity of the arterial wall | τ | 1.333 | ^{[37]} |

Coating drug diffusivity, m^{2} s^{−1} | D_{c} | 1 × 10^{−15} | ^{[38]} |

Coefficient of free diffusivity, m^{2} s^{−1} | D_{free} | 3.65 × 10^{−12} | [^{28}^{,}^{39}] |

Coefficient of effective diffusivity, m^{2} s^{−1} | D_{eff} | 2.15 × 10^{−12} | Our study |

True diffusivity of the free drug, m^{2} s^{−1} | D_{T} | 24 × 10^{−12} | Our study |

Initial drug concentration in the coating, mol m^{−3} | c_{0} | 1 | ^{[23]} |

Net tissue binding capacity, mol m^{−3} | B_{M} | 0.01 | ^{[7]} |

Association rate constant, m^{3} mol^{−1} s^{−1} | k_{a} | 10 | ^{[40]} |

Dissociation rate constant, s^{−1} | k_{d} | 0.01 | ^{[41]} |

Equilibrium dissociation constant, mol m^{−3} | R_{d} | 0.001 | Our study |

Dimensionless Peclet number in the coating | Pe_{c} | 100 | Our study |

Dimensionless Peclet number in the tissue | Pe_{T} | 2 | Our study |

Dimensionless Damköhlar number in the tissue | Da | 40 | Our study |

Dimensionless scaling parameter | ε_{1} | 0.001 | Our study |

Dimensionless scaling parameter | ε_{2} | 100 | Our study |

Distributions of mean free as well as bound drug concentrations for different values of scaling parameter ε_{1}(= *R _{d}*/

*c*

_{0}) are presented in Figure 3a and 3b respectively, and the same for different values of ε

_{2}(=

*c*

_{0}/

*B*), which are depicted in Figure 4a and 4b respectively. Evidently, ε

_{M}_{1}, depending on

*R*(=

_{d}*k*/

_{d}*k*), decreases with a decrease in the dissociation rate constant

_{a}*k*and also with an increase in the association rate constant

_{d}*k*. However, ε

_{a}_{2}increases with decreasing

*B*(keeping

_{M}*c*

_{0}as fixed). Figure 3a shows that the normalized mean free drug concentration decreases with decreasing ε

_{1}for

*Pe*= 2,

_{T}*D*= 40, ε

_{a}_{2}= 100, up to a certain time and, thereafter, no significant changes occurred. It may be justified in the sense that, as ε

_{1}decreases, the rate of reversible binding (

*k*) decreases and/or the rate of forward binding increases, which is lowering the mean concentration of free drug.

_{d}An interesting phenomenon may be observed when the transported free drug is insensitive to flowing blood [*i.e.* modelled as zero-flux lumen-tissue interface condition, the mean concentration of free drug is always higher (zoom inset) as compared to zero-concentration interface condition and approaches a quasi-steady state faster]. Figure 3b shows how the rates of forward as well as reversible binding affect the mean concentration of bound drug within the arterial tissue. It is to be observed that the mean concentration is increased with the decrease of ε_{1}, owing to the increase in the rate of forward binding and/or to the decrease in the rate of reversible binding. Here too, the mean concentration of bound drug appears to be higher in case of zero-flux lumen-tissue interface condition, and saturation of binding sites takes place very rapidly.

Effects of ε_{2} (*i.e.* net tissue binding potential on the mean concentrations of free as well as bound drug) are displayed in Figure 4a and 4b respectively. It may be recalled that ε_{2} increases with decreasing binding potential. The results of these figures indicate that the mean concentration of free drug increases with decreasing binding potential up. to ε_{2} = 100, but the concentration reaches a quasi-steady state for weaker binding capacity (ε_{2} = 1,000) than the other cases. As anticipated, the mean concentration of free drug appears to be higher for zero-flux lumen-tissue interface condition (zoom inset) (Fig. 4a). Similar observations can be made for normalized bound drug concentration but with a distinction. In the case of that free drug for ε_{2} ≤ 100, the quasi-equilibrium is not fully established until *t* = 500. On the other hand, in the case that bound drug for ε_{2} ≥ 10, the quasi-equilibrium is attained very rapidly (Fig. 4b). Furthermore, the effects of zero-flux lumen-tissue interface condition on the spatiotemporal distribution of free and bound drug can be visualized in Figure 5, panels a and b.

The time-evolution of the normalized free and bound drug concentrations at positions *P*_{1} and *P*_{2} (near the left and the perivascular boundaries in Fig. 1) are shown in Figure 6, panels a, b and c, d respectively. The drug distribution profiles indicate that quasi-equilibrium has not yet been fully established until *t* = 500 for small values of *Da* and also that the normalized free drug concentration at *P*_{1} and *P*_{2} decreases with increasing *Da* (Fig. 6a and 6b). Moreover, the reduced local free drug concentration at *P*_{2} with increasing *Da* further reduces the drug dissipation at the perivascular boundary. Similar results on the distribution of bound drug have been shown in Figure 6c and 6d. All of the above observations are in good agreement with those of Hwang *et al.*^{15} who studied the distribution of free drug only. The lack of drug far away from the strut (*i.e.* at *P*_{2}) can be a serious factor for in-stent restenosis occurrence, as high drug concentration at the perivascular end is more important in suppressing restenosis.^{42–44}

The influences of diffusivity (Peclet number) and time-dependent release kinetics of a drug on the temporal variation of free as well as bound drug concentrations within the arterial tissue are portrayed in Figure 7, panels a, b and c, d respectively. Figure 7a and 7b indicate that the free drug concentrations at both *P*_{1} and *P*_{2} do not attain quasi-equilibrium state until *t* = 500, due to the time-dependent release-kinetics. However, if one assumes constant release of drug from a well apposed strut, a quasi-steady state has been fully established at the positions considered. On the contrary, the bound drug concentration attains a quasi-steady state irrespective of time-dependent and constant release of drug. The effect of release kinetics on the spatial distribution of free as well as bound drug can also be visualized clearly in Figure 8, panels a–d, in which heterogeneous distribution and retention of drug are found to be observed throughout the domain and higher concentration of free drug is observed in case of time-dependent release of drug from struts.

Figure 9a shows the concentration profile at a height of one strut thickness for free drug concentration depending on interstrut distance for *Pe _{T}* = 2 and

*Da*= 40 at

*t*= 500. A single peak is noted when the struts are placed one-strut width apart. As the interstrut distance increases, the peak concentration falls and distinct peaks over each strut are observed. The distributions of bound drug concentration are analogous to those of the distributions of free drug qualitatively (Fig. 9b).

When a stent is implanted within an artery, it has a major impact on the structures of the wall and subsequently alters the transport and retention of the drug within the arterial tissue. Any porous material under compression will demonstrate smaller pore ratio than that of its relaxed state, which in turn influences the tortuosity. Arterial property such as porosity (ε) dictates the transport and retention of drug eluted from a pair of struts, as depicted in Figure 10a and 10b respectively. It is observed from these figures that the mean concentrations of free as well as bound drug increases with increasing ε. However, the enhancement is maximum when ε is allowed to change from 0.2 to 0.78, up to *t* = 250 for both drug forms. It is interesting to note that the mean concentration of bound drug attains a quasi-steady state for larger ε (*i.e.* saturation of binding sites takes places very rapidly).

The above observation may be justified in the sense that as the porosity increases, the effective as well as true diffusivity do increase, which eventually leads to expedition of the diffusion process. In a porous media, diffusion takes place in confined tortuous pores and its progression is impeded as the tortuosity increases. Our simulation also demonstrates the fact that a decrease in the mean concentration of free drug with increasing tortuosity (τ) (Fig. 11a) (*i.e.* an inverse relationship between free drug concentration and tortuosity is revealed). A similar pattern is also observed for bound drug (Fig. 11b). The above observation may be justified in the sense that as the tortuosity increases so too does the effective distance over which diffusion has to take place (*i.e.* the progression of diffusion eventually lowering the mean concentration of both drug forms is impeded).

## Conclusions and future work

In this numerical study, a two-dimensional axisymmetric model of drug transport eluted from drug-eluting stents has been developed. We consider a two-species model, capable of predicting the delivery of drug and its retention in the arterial tissue where the transport of free drug is governed by reaction-diffusion process and that of bound drug, assuming complete immobilization in the tissue, by reaction process. Following Tzafriri *et al.*^{19} a second-order dynamic model that describes a saturating reversible binding process has been taken into account. A time-dependent release kinetic is applied on the surface of the struts.^{20,23} As the arterial structures such as porosity and tortuosity affect the diffusion process, their influence on the transport and retention of drug eluted from a pair of struts has been investigated successfully.

Though experimental studies can provide information on release kinetics and other histological information, computational studies can provide detailed predictions of the drug distribution and its retention within the arterial tissue. With the rapid ascent of stent-based drug delivery in the treatment of vascular disease, many important issues concerning drug delivery and its retention in the arterial tissue need to be addressed. Indeed, arterial vessels with disease are the target sites for endovascular intervention, and quantifying the pharmacokinetics for this scenario is needed towards full appreciation of drug-eluting stents and like technologies.

## Abbreviations

- BMS:
bare metal stent

- CAD:
coronary artery disease

- DES:
drug-eluting stent

- ISR:
in-stent restenosis

- PCI:
percutaneous coronary intervention

## Declarations

### Acknowledgement

The authors gratefully acknowledge the careful scrutiny and suggestions of the learned reviewers. The authors gratefully acknowledge the partial financial support from Special Assistance Programme (SAP-III) sponsored by University Grants Commission (UGC), New Delhi, India (Grant No. F.510/3/DRS-III/2015(SAP-I)).

### Conflict of interest

The authors have no conflict of interests related to this publication.

### Authors’ contributions

Formulating the problem (RS, PKM); running the in-house code (RS); writing the manuscript (RS, PKM).