Introduction
The drug-eluting stent (DES) is now a common treatment for atherosclerosis. Placed into the narrowed artery, it serves to reopen the artery and slowly release drug to inhibit cell proliferation. The DES consists of three parts, namely the stent platform, a polymer coating that binds the drug to the stent and releases drug, and the drug. Success of the DES is usually associated with the effective delivery of potent therapeutics to the target site, at a sufficient concentration for a sufficient time and in a biologically active state. Although DESs are now the primary choice for treatment of atherosclerosis, questions still arise in regard to longevity and safety.1
To investigate the drug release mechanism, various models have been developed, including the one-dimensional model,2–4 the two-dimensional or axi-symmetric model and the three-dimensional model.5–13 The drug release process is also affected by some physical properties, namely the geometrical design of the stent itself, the mechanical properties of the stent’s material and the chemical properties of the drug.14,15 A number of experimental and numerical studies on DES have been published in recent years to address such issues as longevity and safety.3,11,16–22
Since binding of the ligand to the cell surface receptors has been amenable to direct experimental investigation for roughly the past three decades,23,24 a sincere attempt has also been made to take into account the binding event that occurs on the cell surface. However, no such attempt has been made to consider the aspects of trafficking process, cell signaling or receptor mobility. The first mathematical model developed to investigate binding kinetics used a constant partition coefficient.25 But, it was too simple to predict the distribution of drug or to consider retention of the drug in the arterial wall. Some authors used a chemical reaction to explain the binding kinetics and also to relate the association and dissociation coefficients by equilibrium coefficients.26,27 A nonlinear reversible chemical reaction is generally accepted to predict the two types of drug forms.28
Vo and Meere29 considered a nonlinear reversible binding model to describe the release of heparin-binding growth factors from an affinity-based delivery system. Groh et al30 adopted a nonlinear reaction model to describe the binding of drug to binding sites within cells and also investigated the chemotherapeutic interaction with the microenvironment of cells in tumour drug delivery. A nonlinear saturable binding model was developed to describe the drug binding in arterial tissue sites, which included two phases of drug: the free drug and the bound drug.18,31
Different binding models namely, nonlinear saturable reversible binding, nonlinear saturable irreversible binding, linear reversible binding and linear irreversible binding, are demonstrated by McGinty and Pontrelli.32 However, it is well established that when the drug binds to the specific receptor (SR), there is also occurrence of non-specific binding caused by the trapping of drug in the extracellular matrix (ECM). Tzafriri et al33 investigated two types of binding–specific binding and non-specific binding–and both of the binding kinetics are modelled by using a nonlinear reversible chemical reaction process. Very recently, McGinty and Pontrelli34 investigated the importance of modelling the specific and non-specific saturable binding in the arterial wall as separate phases.
Very little of this type of work has been carried out considering the specific and non-specific binding of drug,33,34 and some significant aspects remain unanswered. Tzafriri et al33 considered only the transmural transport of eluted drug with constant release kinetics, whereas, McGinty and Pontrelli34 opined that when the ratios of SR to ECM binding site density and of SR to ECM unbinding rate are both small, drugs such as sirolimus and paclitaxel have same modes of action.
Usually, modelling with well-apposed DES represents very early stage after the implantation. The scientific interest for embedment of strut came from the need for accurate and quantitative evaluation of the vessel wall and stent interaction. The degree of stent embedment could be one of the surrogate parameters of the vessel wall-stent interaction.35 In terms of shear stress, the deeper the struts are embedded, the less disturbed the shear stress will be.36 Moreover, deeper penetration into the vessel wall increases direct contact and drug delivery, and reduces recirculation zones, which, in turn, decreases the area of the endothelisation exposed to disturbed flow, thus increasing the probability of endothelisation of adjacent tissues; ultimately, an antiproliferative and anticoagulant flow environment is established, which is the optimal condition for clinical success.35
There seems to be a down side, however, in the sense that embedded struts denote penetration of the cutting edge of the struts through vessel wall, implying larger injury of the vessel that can trigger neointimal hyperplasia.37 Even if the embedment of struts is small, a larger width of strut could contribute to larger amount of vessel wall-stent interaction. Hence, relationship between this injury (i.e. the degree of embedment and neointimal hyperplasia) will be the topic of further study.38 On that theoretical basis, we should not expect an excess of neointima, despite the embedment and injury to the vessel wall.
Keeping all the relevances in mind, the present investigation describes three-phase drug transport phenomena within the arterial wall, namely free drug, specific binding of drug with the SR and non-specific binding of drug in the ECM site from a half-embedded DES. The transport of free drug has been modelled by reaction-diffusion process, whereas the bound drug has been modelled by nonlinear reversible saturable chemical reaction. The arterial wall has been considered as a single homogeneous layer with identical diffusivity property. In this study, the therapeutic domain of length L consists of three struts with centre to centre distance of 0.7 mm (Fig. 1).39,40 At the strut surfaces, a time-dependent drug release kinetic in exponential manner has been accounted for McGinty et al.41 The governing equations together with their physiological realistic boundary conditions are solved numerically in an explicit manner.
Governing equations and boundary conditions
A three-phase model of drug transport has been modelled in the tissue, namely free drug,16 specific binding (SR) and non-specific binding (ECM).42,43 The transport of free drug has been modelled by reaction-diffusion process and binding of drug by nonlinear saturable reversible chemical reaction. Now the dimensionless governing equations for the model considered may be written as
∂cf∂t=1ε2∂2cf∂r2+1r∂cf∂r+ε2∂2cf∂z2−∂cSR∂t−∂cECM∂t,
∂cSR∂t=DaSR{cf(cSRmax−cSR)−cSRkeqSR},
∂cECM∂t=DaECM{cf(cECMmax−cECM)−cECMkeqECM},
where cf, cSR and cECM, normalised by initial strut drug concentration (c0), are the normalised concentration of free drug, specific and non-specific binding drug respectively. r is the dimensionless radial coordinate and z is the dimensionless axial coordinate, normalised by the thickness of the arterial wall (δW) and diameter of the strut (d) respectively. Here, cimax is the dimensionless local density of binding site, normalised by initial strut drug concentration(c0). Dai(kfic0d2Di) and keqi(c0kfikri) are the Damköhler number and equilibrium constant respectively, in which kfi and kri are the association and dissociation reaction rate respectively. Here and in the sequel, i=SR, ECM. Here, t(t¯Dtd2) is the dimensionless time and ε=(δWd).In the tissue, symmetry boundary conditions are applied on the proximal (Γw,in) and the distal (Γw,out) walls:
∂cf∂z=0 on Γw,in and Γw,out
At the perivascular wall, a perfect sink condition for free drug is applied:
cf = 0 at Γwu
A time-dependent drug release kinetic has been used at the strut surfaces:41
cf = e−φt on Γws,
where φ=λd2Dt, λ is the rate of drug release from strut surfaces.At the blood-tissue interface (Γbw), it is assumed that the transported drug is insensitive to luminal flow, which is modelled as zero-flux condition given by
∂cf∂r=0 on Γbw.
Radial coordinate transformation
In order to avoid interpolation error, we transform the therapeutic domain (Fig. 1) into a rectangular one, by making use of the following radial coordinate transformation:
ξ=1+r−RtlΓwu−Rtl where Rtl=Γbw∪Γws.
Considering this transformation, the irregular domain transforms into [0,L]×1,2. Now using the transformation (8) the equations (1)–(7) are transformed as follows:
∂cf∂t=1ε2[∂2cf∂ξ2{1R2+ε2(∂ξ∂z)2} +∂cf∂ξ{1R{R(ξ−1)+Rtl}+ε2∂2ξ∂z2+ε2R2dRdz((ξ−1)dRdz+dRtldz)}] +∂2cf∂z2+2∂ξ∂z∂2cf∂ξ∂z−∂cSR∂t−∂cSR∂t,
∂cSR∂t=DaSR{cf(cSRmax−cSR)−cSRkeqSR},
∂cECM∂t=DaECM{cf(cECMmax−cECM)−cECMkeqECM},
∂cf∂z=0 on z=0 and z=L,
cf = 0 at ξ = 2,
cf = e−φt at ξ = 1 and z ∈ Γws
[∂cf∂ξ]ξ=1=0 for z∈Γbw
Method of solution
The governing equations with the set of initial and boundary conditions are solved numerically by finite-difference scheme. Forward-time centred-space discretisation technique has been made use of in the explicit numerical scheme. Let us describe our finite-difference method in more detail: We denote zi=iδz, tn=nδt, ξk=kδξ, where n refers to the time level, δt the time increment and δz is the space step size along the axial direction. Here, δξ denote space step sizes along the radial direction.
The finite difference approximations of (9–11) are as follows:
cfi,kn+1−cfi,knδt=1ε2[diffcf]i,kn−cSRi,kn+1−cSRi,knδt−cECMi,kn+1−cECMi,knδt,
cSRi,kn+1−cSRi,knδt=[reacSR]i,kn,
cECMi,kn+1−cECMi,knδt=[reacECM]i,kn,
where the expressions for [diffcf]i,kn, [reacSR]i,kn, [reacECM]i,kn are included in the Supplementary Term S1.Now, the discretised equations are solved numerically in an explicit manner. No standard package has been used in the present simulation, rather the computational code has been successfully programmed using FORTRAN language.
Results and discussion
For the purpose of numerical computation of the quantities of significance, the computational domain has been confined to a finite nondimensional length of 25 and solutions are computed with the grid size of 501×61 by making use of the input values presented in Table 1.11,16,31,33,34,39–41,44–46 The simulation concerning grid-independent study was performed for examining the error associated with the grid sizes used and displaced in Figure 2. The transmural variations of normalised bound drug concentration in SR phase at an axial distance of z=11 for three distinct grid sizes almost overlap on one another, which clearly establishes the correctness of the grid sizes used.
Table 1Plausible values of involved parameters (dimensional)
Description | Parameter | Value | Reference |
---|
Diameter of the strut (cm) | d | 0.02 | [31,44] |
Interstrut distance (cm) | zd | 0.06 | [39] |
Wall thickness (cm) | δw | 0.1 | [16,45] |
Free drug diffusion coefficient (cm2·s−1) | Dt | 10−9 | [40,46] |
Initial strut drug concentration (mol·cm−3) | c0 | 10−6 | [11] |
Association rate constant in ECM-bound phase((mol·cm−3·s)−1) | kfECM | 2×106 | [33,34] |
Dissociation rate constant in ECM-bound phase(s−1) | krECM | 5.22×10−6 | [33,34] |
Local density in ECM binding site(mol·cm−3) | c¯ECMmax | 3.63×10−7 | [33,34] |
Association rate constant in SR-bound phase((mol·cm−3·s−1)−1) | kfSR | 8×108 | [33,34] |
Dissociation rate constant in SR-bound phase(s−1) | krSR | 1.6×10−4 | [33,34] |
Local density in SR binding site(mol·cm−3) | c¯SRmax | 3.3×10−9 | [33,34] |
Drug release rate(s−1) | λ | 10−5 | [41] |
Transmural variations of free (NDCFREE), bound in ECM phase (NDCECM) and bound in SR phase (NDCSR) at three distinct nondimensional times have been depicted in Figure 3, panels a–c, respectively. Drug enters the arterial wall at ξ=1 in the free phase, and is rapidly bound to both ECM and SR sites. It is worthy to note that the penetration length of both free and ECM-bound drug increases with increasing time, with saturation of binding sites ultimately taking place. Moreover, the SR-bound drug is absorbed at the adventitial boundary (ξ=2) with the increase of time. The free drug (NDCFREE) and ECM-bound drug(NDCECM) penetrate half of the thickness of the tissue at t=0.5, and although the NDCFREE and NDCECM profiles are similar, drug concentrations within ECM-bound phase are an order of magnitude greater than for the SR-bound phase. Saturation of SR-bound drug takes place with increasing time (Fig. 3c). All the above observations are in conformity with Tzafriri et al33 McGinty and Pontrelli,34 although the embedment of stent has been ignored in the latter study.
Temporal variations of normalised mean drug concentrations are presented in Figure 4. Evidently, the concentrations of free and ECM-bound drug rise from zero to a maximum value and then slowly decay with time; however, the concentration of SR-bound drug steadily increases with time considered. It is worthy to note that the peak of NDCSR occurs later than that of NDCFREE and also more drug is contained within the ECM-bound phase than the other two phases, which are in good agreement with those of Tzafriri et al33 and McGinty and Pontrelli.34 By choosing kfSR=0=krSR, we convert the model considered into the single bound phase model. Here, too, we see that the mean concentration of free drug in case of the single bound phase model is all time higher than that of the double bound phase model, which is due to less conversion of free drug into bound form in the single bound phase model (Fig. 5).
The influence of Damköhler number on the normalised concentrations of free drug, ECM-bound and SR-bound drug in the arterial tissue over a stipulated period of time is displayed in Figure 6a–c respectively. Evidently, the mean concentration of free drug decays rapidly and also the mean concentrations of ECM- and SR-bound drug attain their respective steady state for smaller Damköhler number. Spatial distribution of free, ECM-bound and SR-bound drug concentration can be visualised through Figure 7a–c respectively. These figures clearly establish the impact of time-dependent release kinetics on the transport and retention of drug in stent-based delivery.
Conclusion and study limitation
In this numerical study, we proposed two-dimensional axi-symmetric models of drug transport eluted from a DES. Following Tzafriri et al33 and McGinty and Pontrelli,34 we considered a three-phase model, capable of predicting the time-dependent delivery of free drug and its ECM- and SR-binding in arterial tissue with the binding site actions modelled using a nonlinear reversible saturable chemical reaction. The transport of free drug was modelled as an unsteady reaction-diffusion process, while both the bound drugs were modelled as a reaction process only. The simulated results predict the saturation of binding sites that takes place with increasing time and SR-bound drug is totally absorbed at the adventitial boundary. This study also highlights the fact that the concentration of free drug is always higher when the model reduces to two-phase, in which a single bound phase is considered. Another important observation noted is that the concentration of free drug decays rapidly and the earlier saturation of binding sites take place in case of smaller Damköhler number.
As with the case of any computational model study, our study was based on a number of assumptions made regarding selected parameters (Table 1) and the boundary conditions that were assigned due to non-availability of data in the literature. Although not ideal, derivation of all these parameters from human tissues may not be feasible and hence our assumption may be considered as an approximation towards quantifying arterial pharmacokinetics in stent-based delivery.
Future direction
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. For realistic modelling, the inclusion of luminal flow along with its pulsatility would certainly predict the delivery system one step closer to the real situation. Future directions may also include heterogeneous tissue composition comprising of healthy tissue and regions of fibrous, fibro-fatty, calcified and necrotic core lesions, with varying diffusivities. The effect of porosity and tortuosity of the arterial tissue may not be ruled out for further research.
It is well-established that the presence of stent coating together with its design greatly influence the distribution and retention of drug within the vessel wall but,33,34 for simplicity, the stent-coating system is not taken into account in the present study, which is in agreement with the numerical evidence proposed in.8,26,45 Following McGinty et al,47 Zhu and Braatz,48 the effects of variable porosity, percolation and biodegradability of the stent coating may be distinct topics of future research, and which we intend to include in our further studies.
Abbreviations
- DES:
drug-eluting stent
- SR:
specific receptor
- ECM:
extra cellular matrix
- NDC:
normalised drug concentration
- NDCFREE:
normalised free drug concentration
- NDCSR:
normalised bound drug concentration in SR phase
- NDCECM:
normalised bound drug concentration in ECM phase
Declarations
Acknowledgement
The authors are grateful to the learned reviewers for careful consideration of the manuscript and for valuable suggestions. Partial financial support for carrying out this investigation was provided by the Special Assistance Programme (SAP-III) sponsored by the University Grants Commission, New Delhi, India (Grant No. F.510/3/DRS-III/2015 (SAP-III)).
Conflict of interest
The authors have no conflict of interests related to this publication.
Authors’ contributions
Designing and performing the research as well as writing the paper (APM and PKM).