Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Analytical modeling of drug dynamics induced by eluting stents in the coronary multi-layered curved domain Michele d’Errico∗, Paolo Sammarco, Giuseppe Vairo Department of Civil Engineering and Computer Science Engineering (DICII), Università degli Studi di Roma “Tor Vergata”, via del Politecnico 1, Rome 00133, Italy

a r t i c l e

i n f o

Article history: Received 31 October 2014 Revised 18 June 2015 Accepted 19 June 2015 Available online 7 July 2015 Keywords: Multi-layered porous media Sturm–Liouville problem Drug-eluting stents Advective–diffusive pharmacokinetics Vessel curvature

a b s t r a c t Pharmacokinetics induced by drug eluting stents (DES) in coronary walls is modeled by means of a onedimensional multi-layered model, accounting for vessel curvature and non-homogeneous properties of the arterial tissues. The model includes diffusion mechanisms, advection effects related to plasma ﬁltration through the walls, and bio-chemical drug reactions. A non-classical Sturm–Liouville problem with discontinuous coeﬃcients is derived, whose closed-form analytical solution is obtained via an eigenfunction expansion. Soundness and consistency of the proposed approach are shown by numerical computations based on possible clinical treatments involving both hydrophilic and hydrophobic drugs. The inﬂuence of the main model parameters on drug delivery mechanisms is analyzed, highlighting the effects induced by vessel curvature and yielding comparative indications and useful insights into the concurring mechanisms governing the pharmacokinetics. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Coronary artery disease (CAD) is the most common type of heart disease, and it is one of the leading cause of death worldwide. As sketched in Fig. 1, coronary arteries (the left and the right one) are blood vessels (with an average diameter of the order of 3 − 5 mm) that branch off near the point where the aorta and the left ventricle meet, and they supply oxygen-rich blood to the entire heart muscle [1,2]. Healthy arteries are characterized by a well-deﬁned substructure, comprising of three concentric layers [1–3]: the tunica intima, the tunica media, and the tunica adventitia. The tunica intima is the innermost layer, comprising of a single layer of endothelial cells and a subendothelial layer mainly consisting of connective tissue and collagen ﬁbres. The outer boundary of the tunica intima is surrounded by an elastic tissue with fenestral pores known as the internal elastic lamina (IEL). The medial layer consists primarily of concentric sheets of smooth muscle cells (SMCs) and elastic connective tissue. The tunica media and the tunica adventitia are separated by another thin band of elastic ﬁbres known as the external elastic lamina (EEL). The outermost layer of the artery, the tunica adventitia, comprises of connective tissue ﬁbres and some capillaries. These ﬁbres blend into the surrounding connective tissues and aid in mechanically stabilizing the arteries. Due to the activation of the atherosclerotic process, CAD

∗

Corresponding author. Tel.: +393298737821. E-mail addresses: [email protected] (M. d’Errico), [email protected] (P. Sammarco), [email protected] (G. Vairo). http://dx.doi.org/10.1016/j.mbs.2015.06.016 0025-5564/© 2015 Elsevier Inc. All rights reserved.

may generally induce a pathologic coronary stenosis, that is a blockage or narrowing of coronary arteries. As a result, a limitation of the blood ﬂow to the heart is experienced, producing the possible death of myocardial cells because of the lack of oxygen. Thereby, a heart attack may occur, leading to damage and death of the heart muscle. In many cases, CAD and coronary stenoses can be effectively treated by non-surgical procedures based on coronary angioplasty (i.e., the vessel patency is restored by inﬂating a balloon, placed -in its deﬂated conﬁguration- at the blockage site via a catheter) and stents (i.e., tube-shaped metal frameworks that allow to keep the arteries permanently opened). Nevertheless, angioplasty and stenting techniques may induce tissue damage at the treatment site, triggering a physiological response that generally consists in an inﬂammatory immune response. It increases the risk of thrombotic events and can be accompanied by SMC proliferation in the arterial medial layer, that induces the development of a thick smooth-muscle tissue inside the lumen (neointimal hyperplasia). As a result, the recurrence of stenosis (namely, restenosis) can occur after few months from the clinical treatment [4]. Among available approaches for preventing pathological arterial restenosis after stenting and angioplasty coronary treatments, the most promising one is based on drug-eluting stents (DESs). They combine the mechanical properties of traditional stents with the release into the arterial wall of drugs generally characterized by antithrombotic and anti-proliferative features [3,5]. In particular, a therapeutic drug is loaded within a polymeric ﬁlm that coats the stent struts, and is delivered in the coronary wall after the stent implantation. As a matter of fact, drug-transport dynamics and thereby DES

80

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

Fig. 1. A synoptic sketch of heart, coronary arteries, stenosis, and stenting treatment through DES. 1: aorta; 2: right coronary artery; 3: left coronary artery; 4: pericardium; 5: myocardium; 6: endocardium; 7: drug diffusive transport; 8: drug advective transport; 9: drug-SMC binding reaction. SMC: smooth muscle cells; IEL: internal elastic lamina; EEL: external elastic lamina.

therapeutic effectiveness are strongly affected by several biomechanical factors, related to the stent geometrical features, to physicochemical properties of drug, coating and tissues, as well as to the tissue histological arrangement [3,5,6]. In order to contribute to a better understanding and characterization of dominant effects related to the previously-mentioned factors, many in-vivo and in-vitro experimental studies, as well as (pre-)clinical trials and follow-up analyses have been recently carried out (e.g., [7–17]). Furthermore, in the last few years many theoretical and computational models have been also developed, aiming to provide insights into the drug-release mechanisms, as well as to establish analysis tools and parametric indications for conception and design of effective DESs and therapeutic strategies. In this framework, the arterial wall is generally modeled either as a homogeneous porous medium or as a multi-layered structure. In some cases drug transport due to the plasma ﬁltration through the tissue (activated by physiological transmural pressure gradients) and chemical binding of drug with tissue proteins are also accounted for. Several computational approaches refer to simple one-dimensional formulations [18,19], as well as to more reﬁned bi-dimensional [20–25] and threedimensional [26–30] models. In the context of one-dimensional models, a multi-layered computational approach has been proposed by Lovich and Edelman [18], with the aim of investigating the propagation of a hydrophilic drug, such as heparin. By using a similar strategy and neglecting any effect induced by vessel curvature, Costantini et al. [19] described drug dynamics via a one-dimensional multi-layered model, giving numerical solutions that account for both diffusion and advection effects. They also compared different drugs, conﬁrming the clinical evidence that hydrophobic drugs generally lead to wider therapeutic windows than hydrophilic ones

[31,32]. By referring to the benchmark numerical solution proposed in [25], Vairo et al. [23] considered a consistent two-dimensional numerical model taking into account local hemodynamics, strut and coating geometric effects, as well as the anisotropic drug diffusion through the tissue. In that work an advection-diffusion model was formulated and implemented by using a ﬁnite-element technique, highlighting the inﬂuence of both atherosclerotic plaque and strut pattern. More reﬁned computational approaches were conceived by accounting for three-dimensional domains, either based on simpliﬁed geometry [26] or derived from the simulation of the stent expansion [27–30]. Although numerical models allow us to accurately simulate DES performances by integrating coupled physics, scales and complex geometries, simpliﬁed analytical approaches and closed-form solutions can give useful indications for conception, design and parametric optimization of such clinical devices, providing a synthetic understanding of dominant mechanisms, and therefore allowing to better fulﬁll clinical and therapeutic requirements. Moreover, consistent analytical solutions, when compared with suitable experiments, may allow the various physical parameters regulating pharmacokinetics to be properly estimated via an inverse approach. In the framework of simplifying hypotheses and considering onedimensional multi-layered straight domains, analytical closed-form solutions describing the drug transport due to diffusive mechanisms only have been recently established by Pontrelli and de Monte [33,34]. Previous approaches have been generalized by including effects induced by chemical reactions and plasma ﬁltration [34–36] or by modeling the solid-to-liquid phase-transition of the drug within the coating [34,37]. Nevertheless, to the best of authors’ knowledge,

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

81

Fig. 2. One-dimensional multi-layered curved domain and physics involved in pharmacokinetics model. Notation. β : consumption rate coeﬃcient in the tunica media; P1 : topcoat permeability; P4 : myocardium permeability at the outer boundary. Moreover, for the i-th layer, ri − ri−1 : thickness; ε i : porosity coeﬃcient; Di : drug diffusivity, ki : drug partition coeﬃcient; γ i (for i = 2, 3): drug hindrance coeﬃcient; Ki (for i = 2, 3): Darcy permeability.

available analytical solutions do not account for possible effects induced by the vessel curvature. In the present work a novel closed-form solution describing drug dynamics induced by coronary DESs is proposed, referring to a onedimensional multi-layered porous domain and considering a multiphysics approach. It allows to generalize models established in [33– 36], by accounting for non-homogeneous tissues properties and for vessel curvature. Drug diffusion and advection mechanisms induced by plasma transmural ﬁltration in the porous wall structure are modeled, describing also the inﬂuence of binding reactions with smooth muscle cells. The method of separation of variables is employed, resulting in a non-classical Sturm–Liouville eigenvalue problem with discontinuous coeﬃcients whose analytical closed-form solution is expressed as a series of eigenfunctions, orthogonal in the whole multi-layered domain. Several illustrative computations based on both hydrophilic and hydrophobic drugs are discussed, highlighting the inﬂuence on drug dynamics of vessel curvature, as well as showing the model capability to account in a parametric framework for several biophysical properties of tissues, coating and drug. 2. Mathematical model 2.1. Modeling assumptions Let the coronary artery domain be described by a hollow right circular cylinder, and assume that the drug is uniformly delivered along the radial direction only, from a tubular-like DES. Moreover, let us assume the coronary vessel to be embedded into the myocardium tissue, and the biophysical properties of coating, vessel and surrounding tissues to be axisymmetric with respect to the artery axis. Accordingly, the drug-release dynamics can be conveniently reduced to a one-dimensional (1-D) problem, where the radial coordinate r is the space variable (Fig. 2), such that r = 0 at the vessel axis a. The strut-coating system is taken directly connected to the coronary wall and, in agreement with the numerical evidence proposed in [23], any effect induced by local hemodynamics near the stent struts is not taken into account. As sketched in Fig. 2, the 1-D domain is multi-layered, consisting of a sequence of adjacent and concentric curved subdomains having different thickness and physical properties: the coating (in the following denoted by the index 1), the vessel layers (indexes 2 and 3) and the myocardium tissue (index 4). All regions are modeled as isotropic porous media and, disregarding the tunica intima, as well as internal and external elastic laminae, the arterial wall is herein reduced to a double-layered structure, compris-

ing the tunica media (index 2) and the tunica adventitia (index 3). As a notation rule, mass volume-averaged drug concentration in the i-th layer is denoted in the following as ci (r, t), t being the time variable. It is worth pointing out that, although endothelium and intima layers can play an important role in early stages of coronary artery disease and restenosis process, due to the interventional stenting procedure (generally combined with an angioplasty treatment that pushes the stent against the arterial tissue, often inducing strut penetration into the intimal layer and, thereby, strut embedment within the tissue) and in agreement with well-established modeling approaches and evidence [6,22–25], the coating layer is herein assumed to be directly connected to the tunica media. Nevertheless, by adopting a strategy similar to that proposed by Pontrelli et al. [34], the model herein developed could be effectively generalized by considering other tissue layers, in order to describe, for instance, the possible inﬂuence of tunica intima, elastic laminae, as well as the possible presence of the atherosclerotic plaque. With the aim to establish a proper pharmacokinetics model of the drug transfer from the stent coating to the arterial wall, a suitable description of processes that take place inside coating, arterial wall and surrounding tissues is needed. Several models with very different levels of complexity are available [38]. Nevertheless, in this paper reference is made to well-established approaches that conveniently combine effectiveness and analytical feasibility/simplicity. The coating in a DES is made of a porous biocompatible polymeric matrix that encapsulates a therapeutic drug in solid phase. After the stent implantation into the arterial lumen, coating embeds the surrounding biological ﬂuids (namely, the blood plasma), leading to the drug dissolution. After such a dissolution process, the drug is able to diffuse freely inside the surrounding medium [14,39]. Since the diffusion process after the dissolution is generally much slower than the dissolution itself [14], it is possible to neglect the dissolution phase. Accordingly, the drug into the coating region is herein assumed to be completely dissolved in the ﬂuid phase and drug dynamics in the layer 1 to be governed by a diffusive process [21,23,33]. As regards the arterial wall layers, since their porous structure (related to the microscale arrangement of collagen and elastin ﬁbers, cells, and extracellular matrix [1,40]), and as a result of the physiological transmural pressure difference between the blood and the outer wall tissue, the liquid part of the blood (namely, the plasma) ﬁltrates through the wall, driving an advective drug mass transport mechanism. This process is really slow in comparison with the time variability of pressure gradients produced by the pulsatile nature of the blood ﬂow [41], and thereby its average description through a steady

82

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

approach can be conveniently adopted [42,43]. This occurrence can be justiﬁed by observing that the mean frequency of the pulsatile blood ﬂow is about 2 Hz [44], whereas -for an average transmural pressure difference of about 100 mmHg [41]- the plasma ﬁltrates through the arterial wall with an average velocity of 10−4 –10−5 mm/s [23,33], and thereby with a characteristic ﬁltration time (for a wall thickness of the order of 1–5 × 10−1 mm [1,2]) of about 103 – 104 s. Therefore, in agreement with a number of studies (e.g., [21– 23,27,33,35,36]), drug dynamics in layers 2 and 3 is herein modeled by considering both diffusion and advection mechanisms, with advective transport driven by the steady volume-averaged plasma ﬁltration velocity u(r) in the radial direction, and involving the concentration of the free drug in the ﬂuid phase. As a matter of fact, drug in the tissue is partially bound to speciﬁc receptors of the extracellular matrix and is partially dissolved into the plasma (free drug). Neglecting any transient and reversible binding process, and assuming the equilibrium between free and bound states, free drug concentration in the i-th layer can be expressed as [31] ci /(ki ε i ), where ε i is the porosity coeﬃcient, and ki is the tissue-to-free drug ratio, denoted also as partition coeﬃcient. For a highly-soluble drug (e.g., hydrophilic heparin-based drug) ki ∼ 1, whereas for hydrophobic drug (e.g., taxus-based drug) ki 1. As a therapeutic antiproliferative effect, drug reacts through a metabolic process with receptors on the surface of smooth muscle cells (SMC), and this process is signiﬁcant only inside the media layer [45]. Following Creel et al. [31] and in agreement with the simpliﬁed approach employed in [22,35,36], this mechanism is herein modeled by considering, for the advection–diffusion problem in the layer 2, a consumption term linear in the drug concentration c2 and governed by a ﬁrst order consumption rate coeﬃcient β . Finally, since no signiﬁcant pressure gradient occurs inside the myocardium layer, drug transport driven by advective mechanisms can be considered as fully negligible in layer 4, wherein drug dynamics is thereby described by a diffusive process only. The effective physical properties affecting the drug dynamics in the considered porous multi-layered curved domain are treated as piecewise functions, assumed to be constant in each layer. 2.2. Governing equations Due to the previously-stated assumptions and with reference to the notation introduced in Fig. 2, mass conservation law leads to the following governing equations, wherein vessel curvature is explicitly accounted for:

∂ c1 1 ∂ ∂ c1 + −D1 r = 0, ∂t r ∂r ∂r c2 γ2 ∂ c2 1 ∂ ∂ c2 + +ru −D2 r + β c2 = 0, ∂t r ∂r ∂r k2 ε2 c3 γ3 ∂ c3 1 ∂ ∂ c3 + +ru −D3 r = 0, ∂t r ∂r ∂r k3 ε3 ∂ c4 1 ∂ ∂ c4 −D4 r = 0, + ∂t r ∂r ∂r

r0 < r ≤ r1 r1 < r ≤ r2

drug eluting rate), the interface between coating and media layer at r = r1 is treated as a thin permeable and homogeneous membrane. In agreement with other well-established models (e.g., [19,21,23,33]), it is described following the Kedem–Katchalsky approach, deduced on the basis of the linear thermodynamics of irreversible processes [46]. In particular, the through-the-membrane transport, generally stimulated by the hydrostatic pressure jump pm and the osmotic pressure jump across the membrane (i.e., at r = r1 ), can be described by practical Kedem–Katchalsky equations for diluted and well-mixed binary non-electrolyte solutions [46–49]:

Jv = L p ( pm − σ ), Js = ω + (1 − σ )cs Jv ,

r3 < r ≤ r4

where u = u(r) is the volume-averaged plasma ﬁltration velocity, and ci = ci (r, t ) is the unknown drug-concentration function for the i-th layer, with r ∈ [ri−1 , ri ]. Moreover, in Eq. (1) and referring to the i-th layer, Di denotes the effective drug diffusivity, and γ i (with 0 < γ i ≤ 1) is the hindrance coeﬃcient accounting for drug–tissue frictional effects. The following boundary conditions are enforced. Concentration ﬂux at the strut-coating interface (at r = r0 ) is assumed to be zero. Furthermore, in order to model the possible presence of a topcoat covering the polymeric coating (generally employed to control the

(3)

where Jv is the volume ﬂux, Js the solute ﬂux, c¯s the average value of the solute concentration cs across the membrane, Lp the membrane hydraulic conductivity, ω the solute diffusive permeability coeﬃcient of the membrane, and σ the membrane reﬂection coeﬃcient for a given solute, accounting for membrane selectivity. In the present model, a null hydrostatic pressure jump at the topcoat interface is considered (namely, pm = 0), and the topcoat layer is assumed to behave as an unselective membrane (namely, σ = 0), that is as permeable for both solvent (blood plasma) and solute (drug). As a result, Jv = 0 and the through-the-membrane drug transport mechanisms are governed by the osmotic pressure jump only. Therefore, by expressing in Eq. (3) via the van’t Hoff’s formula [47–49] (i.e., = RT cs , where R is the gas constant, T the absolute temperature, and cs the solute concentration jump), the drug ﬂux at r = r1 results to be proportional, by means of the permeability coeﬃcient P1 = ωRT, to the jump between coating and tunica media of the free ﬂuid-phase drug concentration, i.e. cs = [c1 /(k1 ε1 ) − c2 /(k2 ε2 )]r=r1 . A similar condition is enforced at the outer boundary for r = r4 , by introducing a myocardium permeability coeﬃcient P4 and by assuming a zero free ﬂuid-phase drug concentration for r ≥ r4 . Finally, drug ﬂux is assumed to be continuous at the coating-media interface and at each tissue-tissue interface, and the free ﬂuid-phase drug concentration is assumed to be continuous through tissue–tissue interfaces. As a result, boundary conditions for the differential problem in Eq. (1) read as:

∂ c1 =0 r = r0 ∂r c c2 γ2 c1 ∂ c2 ∂ c1 ∂ c1 2 −u , D1 = P1 D2 = D1 − r = r1 ∂r k2 ε2 ∂r ∂r k2 ε2 k1 ε1 D1

D2

c2 γ2 c3 γ3 ∂ c2 ∂ c3 −u −u = D3 , ∂r k2 ε2 ∂r k3 ε3

D3

c3 γ3 ∂ c3 ∂ c4 −u , = D4 ∂r k3 ε3 ∂r

D4

∂ c4 c4 = −P4 ∂r k4 ε4

(1)

r2 < r ≤ r3

(2)

c2 c3 = k2 ε2 k3 ε3

c3 c4 = k3 ε3 k4 ε4

r = r2 r = r3 r = r4 (4)

In agreement with previously-stated assumptions, at the initial time t = 0 the drug is considered as entirely contained in the coating region with a uniform concentration C0 and completely dissolved in the ﬂuid phase. Therefore, the initial conditions are expressed as:

c1 (r, 0) = C0 ,

r0 < r ≤ r1

ci (r, 0) = 0,

ri−1 < r ≤ ri ,

i = 2, 3, 4.

(5)

Let χ be the dimensionless curvature parameter, χ = (r3 − r1 )/r1 , i.e. the ratio between the wall thickness (r3 − r1 ) and the lumen radius r1 . As a notation rule, problem deﬁned by Eqs. (1)–(5) will be denoted in the following as Pχ . When curvature effects are neglected (namely, when χ → 0+ ) the differential problem Pχ reduces to P0 ,

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

governed by:

K∗ =

∂ c1 ∂ ∂ c1 + −D1 = 0, ∂t ∂r ∂r c2 γ2 ∂ c2 ∂ ∂ c2 + +u −D2 + β c2 = 0, ∂t ∂r ∂r k2 ε2 ∂ c3 ∂ ∂ c3 c3 γ3 −D3 = 0, + +u ∂t ∂r ∂r k3 ε3 ∂ c4 ∂ ∂ c4 + −D4 = 0, ∂t ∂r ∂r

r0 < r ≤ r1

83

K2 K3 ln (1 + χ ) K3 ln

r2 r1

+ K2 ln

r3 .

(11)

r2

Finally, since Eq. (7), the radial ﬁltration velocity results in:

u(r) =

r1 < r ≤ r2 (6)

r2 < r ≤ r3

with boundary and initial conditions that formally read as in Eqs. (4) and (5), respectively. It is worth pointing out that the problem P0 recovers the model formulated in [36] when a double-layered approach is considered, and models proposed in [33–35] when advection and drug-SMC reaction effects are disregarded.

for r1 ≤ r ≤ r3

(12)

where the uniform ﬁltration velocity

u¯ = K ∗ r3 < r ≤ r4

r1 u¯ χ , ln (1 + χ ) r

p μ (r3 − r1 )

(13)

satisﬁes the problem (7) when curvature effects are neglected (namely, for χ → 0+ ), and thereby it has to be employed in P0 [19,21,23,33,36]. If instead curvature effects must to be accounted for, Eq. (12) allows to generalize consistently the treatment of advection contributions in Pχ . 3. Closed-form solution

2.3. Plasma ﬁltration description In order to evaluate advective terms in Eqs. (1)–(6), the volumeaveraged plasma ﬁltration velocity u in the coronary wall has to be determined. As previously-observed, and in agreement with other well-established studies (e.g., [19,21,23,33,36]), the pulsatile nature of the blood ﬂow can be neglected, due to its smaller temporal scales than those relevant to the drug transfer mechanisms. Therefore, let a steady transmural pressure jump p = pb − pe be considered, where pb is an average measure of the blood pressure inside the vessel lumen and pe is the external pressure. Moreover, let the pressure ﬁeld p(r) and the transmural velocity ﬁeld u(r) be assumed to be independent on the drug dynamics and described as piecewise functions, that is u(r) = ui (r) and p(r) = pi (r) (with i = 2, 3) for ri−1 ≤ r ≤ ri . By enforcing a plasma incompressibility condition (i.e., prescribing a solenoidal constraint on u(r) in a planar polar coordinate system) and by treating the volume-averaged velocity u as the 1-D Darcy-type steady transmural ﬂux [42,43], the radial plasma ﬁltration for the i-th layer is described by:

1 d(rui ) =0 r dr Ki dpi ui (r) = − μ dr

with

ri−1 ≤ r ≤ ri ,

ri−1 ≤ r ≤ ri ,

p2 = pb K3 dp3 K2 dp2 = μ dr μ dr

p3 = pe

at r = r2

(9)

pressure distributions in medial (p2 ) and adventitia (p3 ) layers result from Eq. (8) as [50]:

r

K p 1 ln + pb K2 ln (1 + χ) r r K∗ p 3 p3 (r) = ln + pe K3 ln (1 + χ) r

p2 (r) =

with

r1 ≤ r ≤ r2 ,

with

(14)

r2 < r ≤ r3

∂ c4 ∂ 2 c4 D4 ∂ c4 + = D4 , r3 < r ≤ r4 ∂t r ∂r ∂ r2 ∂ c1 =0 r = r0 D1 ∂r c1 ∂ c2 α2 c2 ∂ c1 ∂ c1 c2 − = D1 , D1 = P1 D2 − r = r1 ∂r r ∂r ∂r k1 ε1 k2 ε2 D2

∂ c2 α2 c2 ∂ c3 α3 c3 − = D3 − , ∂r r ∂r r

D3

∂ c3 α3 c3 ∂ c4 − = D4 , ∂r r ∂r

c2 c3 = k2 ε2 k3 ε3

c3 c4 = k3 ε3 k4 ε4

∂ c4 c4 = P4 ∂r k4 ε4

r = r2 r = r3 r = r4 (15)

where the following advection coeﬃcients are introduced:

αi =

u¯ χ r1 γi , ln (1 + χ ) ki εi

i = 2, 3.

(16)

The method of separation of variables [51] is applied, and thereby the solution of Pχ can be written in the form:

at r = r3

∗

∂ c3 ∂ 2 c3 D3 − α3 ∂ c3 + = D3 , ∂t r ∂r ∂ r2

(8)

at r = r1 and

r0 < r ≤ r1

∂ c2 ∂ 2 c2 D2 − α2 ∂ c2 = D2 − β c2 , r1 < r ≤ r2 + ∂t r ∂r ∂ r2

−D4

where μ is the plasma viscosity and Ki is the Darcy’s permeability (assumed to be constant) for the i-th layer. Therefore, by prescribing pressure values at r = r1 and r = r3 , and by enforcing continuity of both pressure and Darcy-type ﬂux at the interface for r = r2 , that is

p2 = p3

∂ c1 ∂ 2 c1 D1 ∂ c1 = D1 , + ∂t r ∂r ∂ r2

(7)

or equivalently by

1 dpi d 2 pi = 0 with + r dr dr2

Given the plasma ﬁltration velocity described by Eq. (12), the governing Eq. (1) and the boundary conditions (4) respectively become:

r2 ≤ r ≤ r3

where K∗ denotes the equivalent Darcy permeability, deﬁned as:

(10)

ci (r, t ) = Ri (r) Ti (t ),

i = 1, . . . , 4.

(17)

It is immediate to prove that the time-dependent solutions are expressed as:

Ti (t ) = e−Di λi t , 2

i = 1, . . . , 4,

(18)

where λi is the separation constant in the i-th layer, whereas the space-dependent solutions Ri (r) have to satisfy a non-classical Sturm–Liouville problem. By imposing the physical requirement T1 (t ) = Ti (t ) (for i = 2, 3, 4) we obtain λ2i = D1 λ2 /Di , where λ = λ1 . Moreover, the space-dependent concentration proﬁles Ri (r)

84

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

(b)

(c)

(d)

Fig. 3. Fraction (M2 + M3 )/M of HBD and TBD in the vessel wall vs. time. Inﬂuence of the plasma ﬁltration through the arterial wall, and inﬂuence of thickness and permeability of the myocardium layer. (a) r4 − r3 = 5 × 10−2 cm, P4 = 10−4 cm/s. (b) r4 − r3 = 5 × 10−2 cm, P4 = 0. (c) r4 − r3 = 1 cm, P4 = 10−4 cm/s. (d) r4 − r3 = 1 cm, P4 = 0. Other model parameters are set as in Table 1.

satisfying Pχ result in (see Appendix A):

R1 (r) = A1 J0 (δ1 r) + B1 Y0 (δ1 r),

R2 (r) =

rs2 [A2 Iν2 (δ2 r) + B2 I−ν2 (δ2 r)], rs2 [A2 Jν2 (δ2 r) + B2 J−ν2 (δ2 r)],

¯ for λ ≤ λ

, ¯ for λ > λ

(19)

R3 (r) = rs3 [A3 Jν3 (δ3 r) + B3 J−ν3 (δ3 r)], where Ai and Bi are integration constants, J0 and Y0 are 0-order Bessel functions of the ﬁrst and second kind, and Jq and Iq indicate q-order classical and modiﬁed Bessel functions of the ﬁrst kind, respectively, with q ∈ R (R being the space of the real numbers). In Eq. (19) the following quantities have been introduced:

si = ν i =

2Di

,

with i = 2, 3,

λ¯ =

ci (r, t ) =

C k e−σk t Rki (r),

for ri−1 < r ≤ ri ,

(21)

k

R4 (r) = A4 J0 (δ4 r) + B4 Y0 (δ4 r),

αi

homogeneous, leading to an eigenvalue problem whose solution is an inﬁnite series of ordered eigenvalues λk . For each eigenvalue λk the integration constants Aki and Bki are determined in cascade from Eq. (15), as dependent on an arbitrary constant Ck only. Hence, the solution of Pχ can be written as (for each i = 1, . . . , 4):

β D1

,

D1 λ2 − β D1 . δi (λ) = λ , with i = 1, 3, 4, δ2 (λ) = Di D2

(20)

Coeﬃcients Ai and Bi (i = 1, . . . , 4) are determined by imposing boundary conditions (15). The resulting 8 × 8 algebraic system is

where σk = [λk ]2 D1 , and Rki (r) (for i = 1, . . . , 4) are given by Eq. (19), with λ = λk . The set of functions Rki (r) (i = 1, . . . , 4) deﬁnes the k-th piece-wise eigenfunction Rk (r), corresponding to the eigenvalue λk and orthogonal in r ∈ [r0 , r4 ] to Rm (r) for any m = k. As reported in Appendix B for the sake of compactness, constants Ck can be determined by enforcing the orthogonality property of the piece-wise eigenfunctions Rk (r) and making use of the initial conditions (5). To this aim, these latter can be expressed as:

c1 (r, 0) =

C k Rk1 (r) = C0 ,

r 0 < r ≤ r1

k

ci (r, 0) =

C k Rki (r) = 0,

ri−1 < r ≤ ri ,

i = 2, 3, 4.

(22)

k

It is worth pointing out that, the proposed solution can be directly related to the Péclet dimensionless ratio (relating advective and

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

85

(a)

(b)

(c)

(d)

Fig. 4. Fraction RD of HBD and TBD reacting with smooth muscle cells in the media layer vs. time. Inﬂuence of the plasma ﬁltration through the arterial wall, and inﬂuence of thickness and permeability of the myocardium layer. (a) r4 − r3 = 5 × 10−2 cm, P4 = 10−4 cm/s. (b) r4 − r3 = 5 × 10−2 cm, P4 = 0. (c) r4 − r3 = 1 cm, P4 = 10−4 cm/s. (d) r4 − r3 = 1 cm, P4 = 0. Other model parameters are set as in Table 1.

diffusive effects) and to the ﬁrst and second Damköhler numbers (the ﬁrst relating reaction rate and advective effects, the second one relating reaction rate and diffusive effects) [36,52]. When the curvature inﬂuence is neglected, similar arguments allow to prove that the differential problem P0 admits the solution form in Eq. (21), with Rki (r) described, for the sake of brevity, in Appendix C.

•

•

Mi (t ) = 2π •

4. Applicative examples In order to verify soundness and consistency of the proposed closed-form solution, different cases have been considered for computational applications, referring to different types of drug and different values of model parameters. In detail, the proposed model has been applied by referring to a hydrophilic heparin-based drug (HBD) and to a hydrophobic taxus-based drug (TBD) [20,31,54]. Reference values herein employed for model parameters, and summarized in Table 1, have been chosen in agreement with well-established studies [2,19,21–23,25,27,31,33,34,36,53]. It is worth observing that values herein adopted for lumen radius and layer thicknesses are typical of a moderately thickened coronary artery [2,25] and correspond to χ = 1/3. Applicative examples are numerically computed by truncating the series in Eq. (21) at the ﬁrst 50 terms, and numerical results are presented in the following by referring to:

drug concentration proﬁles in the tissue, normalized by the initial value C0 in the coating; drug mass Mi (t) per unit length in the i-th layer, normalized by the total drug mass M per unit length at t = 0, with

ri ri−1

r ci (r, t ) dr,

M = C0 π (r12 − r02 ) ;

(23)

drug mass per unit length undergoing the consumption process due to binding reactions with smooth muscle cells in the medial layer, normalized by the total drug mass M per unit length at t = 0:

RD(t ) =

2π M

t 0

r2 r1

β r c2 (r, τ ) drdτ .

(24)

Referring to the solution of the problem Pχ , and with the aim to elucidate the inﬂuence of the advection-based drug transport inside the arterial wall, some calculations have been carried out taking into account or not (i.e., by assuming u¯ = 0 in Eq. (12)) plasma ﬁltration mechanisms. In detail, Fig. 3 shows results concerning the time evolution of the drug fraction in the vessel wall, i.e. (M2 + M3 )/M; Fig. 4 depicts the reacting drug amount RD(t); Fig. 5 illustrates the drug fraction M4 (t)/M in the myocardium region. Those ﬁgures depict also results computed by considering different values (with respect to the reference ones in Table 1) of the myocardium permeability at the outer boundary (namely, P4 ) and of the myocardium thickness (r4 − r3 ).

86

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96 Table 1 Reference values of model parameters employed in applicative examples. HBD: heparinbased drug. TBD: taxus-based drug. Subdomain

1. Coating

2. Media

3. Adventitia

4. Myocardium

Parameter

Reference value

Lumen radius, r1 [mm] Blood plasma viscosity, μ [g/cm/s] Transmural pressure difference, p [mmHg] Thickness, (r1 − r0 ) [μm] Porosity, ε 1 HBD diffusivity, D1 [cm2 /s] TBD diffusivity, D1 [cm2 /s] HBD partition coeﬃcient, k1 TBD partition coeﬃcient, k1 Topcoat permeability, P1 [cm/s] Thickness, (r2 − r1 ) [cm] Porosity, ε 2 Darcy permeability, K2 [cm2 ] HBD diffusivity, D2 [cm2 /s] TBD diffusivity, D2 [cm2 /s] HBD partition coeﬃcient, k2 TBD partition coeﬃcient, k2 Drug hindrance coeﬃcient, γ 2 Drug consumption rate, β [1/s] Thickness, (r3 − r2 ) [cm] Porosity, ε 3 Darcy permeability, K3 [cm2 ] HBD diffusivity, D3 [cm2 /s] TBD diffusivity, D3 [cm2 /s] HBD partition coeﬃcient, k3 TBD partition coeﬃcient, k3 Drug hindrance coeﬃcient, γ 3 Thickness, (r4 − r3 ) [cm] Porosity, ε 4 HBD diffusivity, D4 [cm2 /s] TBD diffusivity, D4 [cm2 /s] HBD partition coeﬃcient, k4 TBD partition coeﬃcient, k4 Outer-boundary permeability, P4 [cm/s]

1.5a , b , c 0.72 × 10−2 c , d , e , f 100.0c , f , g 5.0e , d , h , i , j 0.1c , d , h , j 10−10 d , h , i , j 10−10 d , f 1.0d , h , j 1.0c , d , f 10−6 d , f , h , i , j , k 4.0 × 10−2 a , f 0.61c , d , e , f , h , j 2.0 × 10−4 c , d , e 7.7 × 10−8 d , f , j 2.6 × 10−8 d , f 1.0d , f , h , i , j , l 20.0c , d , e , f , l 1.0c , d , e , f , l 10−4 f , k , l , m 10−2 a , f 0.85f , j , l 2.0 × 10−4 c , d , e 12 × 10−8 f , j , m 4.0 × 10−8 f , j , m 1.0d , f , h , i , j , l 20.0c , d , e , f , l 1.0c , d , e , f , l 5.0×10−2 f 0.61a , f 7.7 × 10−8 f , m 2.6 × 10−8 f , m 1.0f , l 20.0f , l 10−4 f

a

Waller et al. (1992) [2], Mongrain et al. (2007) [25], c Vairo et al. (2010) [23], d Zunino (2004) [21], e Migliavacca et al. (2007) [27], f Costantini et al. (2008) [19], g Meyer et al. (1996) [41], h Grassi et al. (2009) [22], i Pontrelli and de Monte (2007) [33], j Pontrelli et al. (2013) [34], k Pontrelli and de Monte (2009) [36], l Creel et al. (2000) [31], m Ai and Vafai (2006) [53]. b

Table 2 Root mean square errors ( ) over the time interval [0, 3.6 × 104 ]s (corresponding to the ﬁrst 10 h after the stent implantation) of the P0 -based solutions with respect to those of Pχ , for χ = 1/3, in terms of the drug fractions (M2 + M3 )/M (index M23), M4 /M (index M4) and RD (index RD). Inﬂuence of the main model parameters. Model parameters not speciﬁed, as well as the reference model, are set by considering values summarized in Table 1. HBD Reference model D1 = 10−11 P1 = 10−8 β = 2.5 × 10−4 r1 − r0 = 5 × 10−3 r4 − r3 = 1.0 P4 = 0.0

[cm2 /s] [cm/s] [1/s] [cm] [cm] [cm/s]

TBD

M23

M4

RD

M23

M4

RD

0.0094

0.0114

0.0117

0.0012

0.0061

0.0068

0.0170 0.0198 0.0131 0.0901 0.0230 0.0081

0.0922 0.0338 0.0120 0.1280 0.0075 0.0160

0.0197 0.0421 0.0508 0.0898 0.0226 0.0079

0.0104 0.0012 0.0063 0.0145 0.0057 0.0012

0.0102 0.0057 0.0028 0.0020 0.0220 0.0074

0.0284 0.0108 0.0035 0.0312 0.0121 0.0071

In order to furnish further indications about the sensitivity of simulated pharmacokinetics on model parameters, Figs. 6–10 show results in terms of the drug fraction in the arterial wall (M2 + M3 )/M and of the reacting-drug fraction RD, computed for different values of topcoat permeability P1 , coating diffusivity D1 , coating thickness

r1 − r0 (keeping constant the initial concentration C0 ), reaction rate coeﬃcient β , and equivalent Darcy’s permeability K∗ , respectively. As a notation rule, in the following the peak value of the drug amount inside the arterial tissue (i.e., the peak value computed for the time evolution of (M2 + M3 )/M) will be denoted as DP; the time

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

87

(a)

(b)

(c)

(d)

Fig. 5. Fraction M4 /M of HBD and TBD in the myocardium region vs. time. Inﬂuence of the plasma ﬁltration through the arterial wall, and inﬂuence of thickness and permeability of the myocardium layer. (a) r4 − r3 = 5 × 10−2 cm, P4 = 10−4 cm/s. (b) r4 − r3 = 5 × 10−2 cm, P4 = 0. (c) r4 − r3 = 1 cm, P4 = 10−4 cm/s. (d) r4 − r3 = 1 cm, P4 = 0. Other model parameters are set as in Table 1.

interval in which the fraction (M2 + M3 )/M is greater than 10% will be referred to as the therapeutic window and denoted as TW; ﬁnally, the asymptotic value (for t → +∞) of the reacting drug fraction in the medial layer will be indicated as RD∞ . The inﬂuence of curvature effects is analyzed in Figs. 11–15, wherein drug concentration proﬁles based on the solution of Pχ and those obtained via P0 (namely, for χ → 0+ ) are compared at t = 104 s for different values of some model parameters. Finally, by considering the time interval [0, 3.6 × 104 ] s (corresponding to the ﬁrst 10 h after the stent implantation) and different values for the main model parameters, the inﬂuence of the vessel curvature on integral measures introduced in Eqs. (23) and (24) is focused in Table 2, wherein the root mean square errors ( ) of the P0 -based solutions with respect to those of Pχ are summarized. 5. Discussion Results depicted in Fig. 3 clearly highlight that advective contributions can play a signiﬁcant role in the case of hydrophilic drugs (such as HBD), reducing the residence time in the arterial wall in comparison with the no-plasma-ﬁltration case. Accordingly, in the case of heparin-based drugs a purely diffusive approach (as adopted for instance in [33,34]), that does not account for advective drug transport, can give some inaccurate indications. On the contrary, advective transport does not induce signiﬁcant effects on the release dynamics of a taxus-based drug (such the paclitaxel), due to the hydrophobic character of the drug, related to high values of the

partition coeﬃcient and thereby to small amounts of free drug. Moreover, results in Fig. 3 show that for both drug types the peaks (DPs) of the drug amount in the arterial tissue are almost similar, but the effective therapeutic window (TW) is wider in TBD-based simulations than for HBD, fully in agreement with modeling results proposed in [19,21] and with well-established clinical studies [31,32,55]. In detail, proposed results suggest that TBD tends to accumulate more easily then HBD in the arterial wall, because its tendency to bind with the tissue and its lower diffusivity. Such an occurrence supports also the clinical evidence that, although heparin-coated stents are reported to have lower rates of subacute thrombosis, they exhibit a worse impact on restenosis events than paclitaxel-based DESs [55]. Results proposed in Fig. 3 also show the combined inﬂuence of myocardium permeability and thickness. In detail, the variation of P4 and (r4 − r3 ) with respect to their reference values slightly affects the drug residence in the arterial wall, and only in the post-peak phase. In the limit case P4 = 0, since the drug ﬂux at r = r4 is fully prevented, a change in sign of the concentration gradient occurs inside the myocardium region at a certain time. As a result, the drug begins to diffuse towards the vessel lumen until the whole amount of the drug is exhausted by bio-chemical reactions. Such a diffusive reverse ﬂux is delayed as the myocardium layer is thick. Thereby, the reverse-ﬂux effect tends to be negligible when the myocardium thickness is large enough and it fully disappears when myocardium permeability is different from zero. The greater inﬂuence of advection effects for HBD than TBD clearly appears also by analyzing results relevant to the amount of the

88

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

(b)

(c)

(d)

Fig. 6. Fraction (M2 + M3 )/M in the vessel wall and fraction RD vs. time for different values of the topcoat permeability P1 . Other model parameters are set as in Table 1. (a,b) HBD. (c,d) TBD.

reacting drug in the media layer (Fig. 4). In fact, when reference values are considered (P4 = 10−4 cm/s and r4 − r3 = 5 × 10−2 cm) and referring to HBD-based results, the asymptotic amount of the drug consumed by binding reactions with smooth muscle cells (RD∞ ) is about 35% when advection transport is considered, and 85% when ﬁltration effects are disregarded. Such a difference is basically related to the fact that plasma ﬁltration forces the drug to leave the vessel wall and to concentrate in the myocardium region (Fig. 5), leading to high values of outgoing ﬂuxes. On the contrary, for the TBD case the asymptotic fraction of the reacting drug is about 93–95%, and such a percentage is slightly affected by advection contributions. In this case, only a small drug amount is transferred by dominant diffusion mechanisms to the myocardium layer and thereby is lost out (Fig. 5). By reducing myocardium permeability up to the limit value P4 = 0, a different time evolution of the reacting drug amount is numerically experienced with respect to the previous case. In particular, the fractions of the consumed drug for both HBD and TBD seem to tend asymptotically to the unit value, with a dynamics that is almost unaffected (respectively, highly affected) by the plasma ﬁltration for the TBD case (respectively, for HBD). Such an inﬂuence of the myocardium permeability is strongly mitigated when a thick myocardium region (r4 − r3 = 1 cm) is considered. It is worth pointing out that, within the limitations of the proposed approach, results obtained by considering a thin myocardium layer could be referred to the portion of the coronary artery facing at the pericardium side (see Fig. 1). Similarly, results relevant to a thick myocardium region could be considered as representative of the portion of the coronary artery facing at the endocardium side.

In order to investigate the inﬂuence of coating properties on drug release dynamics, different values of permeability P1 at the interface between coat and arterial tissue (generally controlled by a suitable thin topcoat layer) and of drug diffusivity D1 in the coating (directly related to the coating material features) have been considered. As regards the inﬂuence of the topcoat permeability, results presented in Fig. 6 show that P1 mainly affects peak values (DPs) of drug concentrations in the coronary tissue and allows to regulate the time extent of the therapeutic window. In particular, a P1 -reduction leads to the reduction in DPs (up to about 50% for HBD and 40% for TBD, passing from P1 = 10−6 cm/s to P1 = 10−8 cm/s), producing their shifting towards higher time values and inducing also wider time intervals in which a signiﬁcant drug amount resides in the tissue. In fact, referring to TBD (respectively, HBD), the therapeutic window is about 9.3 h (respectively, 8 h) for P1 = 10−8 cm/s, and it becomes about 7 h (respectively, 5.5 h) for P1 = 10−6 cm/s. Proposed results seem also to indicate that for values of P1 greater than a given level (10−6 cm/s, for the cases under investigation) no signiﬁcant effect on drug dynamics inside the tissue wall is experienced. Moreover, P1 does not affect the asymptotic values of the drug amount RD reacting with smooth muscle cells but, due to the inﬂuence on the drug-release rate from the coating, it slightly modiﬁes the time evolution of RD. As the results depicted in Fig. 7 conﬁrm, a similar inﬂuence on the drug dynamics is observed by referring to the drug diffusivity D1 in the coating. In detail, RD∞ is practically not affected by D1 , but a reduction of the drug diffusivity in the coating induces a slower reacting process and thereby a better in-time distribution of the drug inside the arterial tissue, consisting in a reduction of the concentration peaks (DPs)

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

89

(a)

(b)

(c)

(d)

Fig. 7. Fraction (M2 + M3 )/M in the vessel wall and fraction RD vs. time for different values of the drug diffusivity D1 in the coating. Other model parameters are set as in Table 1. (a,b) HBD. (c,d) TBD.

(a)

(b)

Fig. 8. (a) Fraction (M2 + M3 )/M in the vessel wall and (b) fraction RD vs. time for different values of the coating thickness (r1 − r0 ), and considering HBD and TBD. Other model parameters are set as in Table 1.

and in an enlargement of TW. Referring to TBD (respectively, HBD), and passing from D1 = 10−9 cm2 /s to D1 = 10−11 cm2 /s, DP reduces of about 60% (respectively, 68%), and TW increases from 7 to 10 h (respectively, from 5.5 to 8.6 h). With reference to Fig. 8, a signiﬁcant inﬂuence has been experienced on the drug dynamics inside the arterial tissue when the coating thickness has been varied (keeping constant the initial concentration C0 ). As a matter of fact, by increasing the coating thickness with respect to its reference value, a slower drug release occurs for both

HBD and TBD. As a result, the time evolution of drug fractions in media and adventitia layers are more homogeneous, with larger values of the therapeutic window and smaller peaks. Moreover, a more gradual drug consumption process is observed, characterized by asymptotic values that seem to be comparable with those relevant to the case of a smaller coating thickness, but that tend to be reached in a time interval larger than 30 h. The inﬂuence of the reaction rate coeﬃcient β is investigated in Fig. 9. Proposed results show that β affects drug dynamics of

90

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

(b)

(c)

(d)

Fig. 9. Fraction (M2 + M3 )/M in the vessel wall and fraction RD vs. time for different values of the reaction rate coeﬃcient β . Other model parameters are set as in Table 1. (a,b) HBD. (c,d) TBD.

TBD more than HBD. For TBD, an increase of the reaction rate from 5 × 10−5 /s to 2.5 × 10−4 /s leads to a DP reduction of about 30% (respectively, 20% for HBD) and to a TW reduction of about 80% (respectively, 50% for HBD). On the contrary, for the cases under investigation, an increase in β allows to increase the reacting drug fraction up to RD∞ = 93% for TBD and RD∞ = 63% for HBD. Therefore, proposed results clearly highlight that the variation of β and of the coating thickness produce some counteracting effects on the DES therapeutic effectiveness: when the therapeutic window (i.e., the duration of the therapeutic action) increases then the reacting drug fraction (i.e., the antiproliferative effectiveness) reduces, and viceversa. Accordingly, fully in agreement with clinical evidence (e.g., [6–9,12,31,32,55]) and other modeling approaches (e.g., [18–30,34– 36]), applicative examples herein discussed suggest that a suitable therapeutic effectiveness can be obtained by choosing coating material and drug type such that drug diffusivity D1 and topcoat permeability P1 are small enough, as well as by prescribing values for coating thickness and drug reaction rate within certain optimal ranges. As a matter of fact, the distribution of the drug inside the arterial wall is also strongly affected by the properties of the tissue constituting the wall. Even if arterial features cannot be considered as design parameters for a DES, it can be useful to elucidate the inﬂuence of the equivalent Darcy permeability K∗ of the coronary tissue (related to patient-speciﬁc features), because it directly occurs in advective mechanisms related to the plasma ﬁltration inside the wall. Results proposed in Fig. 10 show that, for values herein adopted and within physiological ranges, K∗ practically does not affect the time evolution

of the drug in the arterial wall. Nevertheless, as a result of different rates of the ﬁltration process, K∗ may induce an effect on the asymptotic values of the reacting drug amount. In detail, for HBD (for which advection mechanisms can play a signiﬁcant role), the reduction of K∗ from 2 × 10−14 to 8 × 10−15 cm2 induces an increment in RD∞ of about 68% (less than 2% for TBD). Proposed approach allows to establish the inﬂuence of vessel curvature on drug dynamics. As shown in Figs. 11–15 at t = 104 s, local concentration proﬁles in the tissue computed by solving the problem P0 may exhibit signiﬁcant differences with respect to those relevant to Pχ , depending on drug type and generally resulting in higher curvature inﬂuence for HBD. In detail, drug concentration patterns related to P0 describe upper bounds for Pχ solutions, with differences that are strongly related to the values employed for model parameters. For instance, the analysis of Fig. 11 suggests that the myocardium permeability P4 has a signiﬁcant effect (which is mitigated for a thick myocardium layer) mainly for the HBD case and localized inside the myocardium region (i.e., for r > 0.05 cm). Nevertheless, as also conﬁrmed by results in Table 2, an increase in P4 tends to reduce discrepancies among P0 and Pχ solutions for the drug amount in coronary tissue layers, and tends to enhance curvature effects for the reacting drug fraction in the tunica media. Moreover, by considering a reaction rate coeﬃcient β equal to 10−4 /s, the P0 solution for HBD overestimates the Pχ one up to about 30% for HBD in both adventitia and myocardium layers, and up to about 10% in the medial region for TBD (Fig. 15). Similarly, by setting the coating thickness equal to 5 × 10−3 cm, the P0 -based drug concentration proﬁle computed at t = 104 s for HBD (Fig. 14) overestimates the Pχ one up to about 60% in the

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

91

(a)

(b)

(c)

(d)

Fig. 10. Fraction (M2 + M3 )/M in the vessel wall and fraction RD vs. time for different values of the equivalent Darcy’s permeability K∗ for the vessel wall. Other model parameters are set as in Table 1. (a,b) HBD. (c,d) TBD.

medial layer, and 70% in both adventitia and myocardium regions (up to about 15% in medial layer for TBD). Nevertheless, dealing with the integral quantities deﬁned in Eqs. (23) and (24), the differences among study cases accounting or not for curvature effects signiﬁcantly reduce. Such an evidence, highlighted in Table 2, indicates that effects related to vessel curvature can be considered as basically negligible if one aims to estimate quantities directly useful for clinical applications (i.e., for designing a therapeutic strategy and/or for giving an estimate of DES global performance). On the contrary, curvature effects may signiﬁcantly affect local and sublocal drug transfer mechanisms, that highly depend on both local concentration gradients and possible discontinuities of total drug concentration at tissue interfaces. Results presented in Figs. 11–15 reveal that the proposed approach is able to recover these possible concentration jumps at the internal tissue interfaces, as a result of the continuity conditions on the free-drug concentrations (i.e., related to the drug amount dissolved into the plasma and not bound to tissue receptors) and of the modeled inhomogeneous tissue properties (see Section 2.1 and Eq. (4)). In particular, concentration discontinuities at the media-adventitia and adventitia-myocardium interfaces may be signiﬁcant when not-negligible advection effects are experienced (i.e., for hydrophilic drugs). Figs. 11–15 also highlight really different drug distributions when HBD and TBD cases are compared. Such an occurrence conﬁrms again the different inﬂuence of concurring drug transfer mechanisms for these two types of drug, resulting in a TBD dynamics mainly driven by dominant diffusive mechanisms, and in not-negligible advection effects for HBD.

Some limitations of the modeling assumptions adopted in this study can be found. In detail, the model is formulated in the framework of a one-dimensional approach, along the radial direction of the coronary wall, disregarding any circumferential and transversal (i.e., along the vessel axis) effect. Moreover, the proposed closed-form solution has been obtained under the assumption of an embedded conﬁguration of stent struts, that is disregarding the inﬂuence of possible plaque occurrence and intima layer, as well as disregarding any effect related to internal and external elastic laminae belonging to the arterial substructure. Moreover, drug is assumed to be perfectly dissolved in the ﬂuid phase inside the coating and characterized by a uniform concentration at the initial time. Thereby, solid-to-liquid phase transition of the drug in the coating layer, as well as possible initial concentration gradients inside the polymeric substrate (usually associated to the production technology and/or to a speciﬁc requirement for controlling the time evolution of the release process) have been neglected. Drug reaction with smooth muscle cells have been treated by a simpliﬁed linear approach, as well as effects associated to local hemodynamics near stent struts, “washing” mechanisms associated with blood ﬂow in the lumen, and the pulsatile nature of blood ﬂow and transmural pressure gradients have been also disregarded. Nevertheless, in agreement with other modeling and numerical studies (e.g., [20–30,33–36]), the present assumptions can be retained as simply enough and effective in order to deduce via a closedform approach useful clinical indications and insights into the drug transfer mechanisms and pharmacokinetics. In fact, although the application of such an approach in a clinical/experimental context is beyond the purposes of the present work, proposed results -obtained

92

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

(b)

(c)

(d)

Fig. 11. Drug concentration proﬁles in the tissue, normalized by C0 , for HBD and TBD at t = 104 s. Comparison among solutions of Pχ (for χ = 1/3) and P0 (χ → 0+ ). Inﬂuence of thickness and permeability of the myocardium layer. (a) r4 − r3 = 5 × 10−2 cm, P4 = 10−4 cm/s. (b) r4 − r3 = 5 × 10−2 cm, P4 = 0. (c) r4 − r3 = 1 cm, P4 = 10−4 cm/s. (d) r4 − r3 = 1 cm, P4 = 0. Other model parameters are set as in Table 1.

(a)

(b)

Fig. 12. Drug concentration proﬁles in the tissue, normalized by C0 , for (a) HBD and (b) TBD at t = 104 s. Comparison among solutions of Pχ (for χ = 1/3) and P0 (χ → 0+ ). Inﬂuence of the topcoat permeability P1 [cm/s]. Other model parameters are set as in Table 1.

by considering different sets of values for model parameters within typical physiological and clinical ranges- reveal the model capable to properly describe well-established evidence from in-vivo and invitro experimental analyses, as well as from (pre-)clinical trials and follow-up studies (e.g., [7–17,31,32,55]). Further quantitative comparisons among present results and available experimental data have

not been performed, due to the wide patient-dependent variability of many model parameters. In fact, they often are not available and/or are not measured when usual clinical trials and follow-up procedures are considered. Accordingly, a proper quantitative model validation, although performed on the basis of integral measures and analyses allowable by typical clinical trials, could be accomplished only

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

93

(b)

Fig. 13. Drug concentration proﬁles in the tissue, normalized by C0 , for (a) HBD and (b) TBD at t = 104 s. Comparison among solutions of Pχ (for χ = 1/3) and P0 (χ → 0+ ). Inﬂuence of the drug diffusivity D1 [cm2 /s] in the coating. Other model parameters are set as in Table 1.

(a)

(b)

Fig. 14. Drug concentration proﬁles in the tissue, normalized by C0 , for (a) HBD and (b) TBD at t = 10 s. Comparison among solutions of Pχ (for χ = 1/3) and P0 (χ → 0+ ). Inﬂuence of the coating thickness r1 − r0 [cm]. Other model parameters are set as in Table 1. 4

if proper experimental measures (via in-vitro, in-vivo or ex-vivo procedures) of the main involved physiological parameters are suitably performed. With the aim to enhance the proposed approach, trying to overcome some of the previously-discussed limitations, future works will be devoted to generalize present closed-form solution in agreement with the strategy proposed by Pontrelli et al. [34], by considering a more reﬁned multi-layered domain including also tunica intima, elastic laminae, as well as the possible presence of the atherosclerotic plaque. Moreover, following the model recently proposed in [34,37], the solid-to-liquid phase-transition dynamics in the coating layer will be coupled to the propagation model of the drug in the arterial wall, possibly accounting also for dominant non-linearities affecting drugSMC reaction process. Moreover, some circumferential effects and the possible inﬂuence of local hemodynamics will be modeled by removing the present assumption of a uniformly-distributed polymeric coating, and by considering, under suitable periodicity constraints, a bi-dimensional formulation that allows for a more realistic description of the discrete distribution of stent struts. 6. Concluding remarks In the present work a novel analytical closed-form solution describing pharmacokinetics induced by coronary drug-eluting stents is proposed. The model is based on a one-dimensional multilayered

approach and is formulated by including effects related to vessel curvature. It accounts for non-homogeneous tissue properties, diffusive mechanisms, advection effects induced by plasma transmural ﬁltration inside the vessel wall, binding reactions of drug with smooth muscle cells, as well as for the possible occurrence of a topcoat layer at the coating-tissue interface and of the surrounding myocardium tissue. The closed-form solution has been established by using the method of separation of variables, resulting in a non-classical modiﬁed Sturm–Liouville eigenvalue problem solved in terms of a series of piece-wise eigenfunctions. The obtained solution has been employed to perform simpliﬁed computational experiments, referring to study cases relevant to heparin-based (hydrophilic) and taxus-based (hydrophobic) drugs. Proposed results clearly show that present analytical solution is able to describe the inﬂuence of the drug type on the transport dynamics. Moreover, computed results highlight that advective contributions play a crucial role in the transport mechanisms governing the dynamics of hydrophilic drugs, whereas their inﬂuence for hydrophobic drugs can be retained as marginal. The model enables to quantify also the inﬂuence of vessel curvature. In particular, discussed results have shown that disregarding curvature effects can induce signiﬁcant discrepancies in local drug concentration proﬁles within the tissues. The proper descriptions of these effects, although they do not mainly affect integral measures of the drug amount inside the multi-layered domain, can give useful

94

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

(a)

(b)

Fig. 15. Drug concentration proﬁles in the tissue, normalized by C0 , for (a) HBD and (b) TBD at t = 104 s. Comparison among solutions of Pχ (for χ = 1/3) and P0 (χ → 0+ ). Inﬂuence of the reaction rate coeﬃcient β [1/s]. Other model parameters are set as in Table 1.

insights into the concurring mechanisms that locally govern the drug transfer dynamics and pharmacokinetics. Applicative examples herein discussed prove that the proposed closed-form solution allows for parametric investigations aiming to quantify the inﬂuence of many biophysical properties of coating, drug and tissues on the clinical effectiveness of the therapeutic treatment, opening towards the possibility of conception and design of novel DES devices, as well as of novel optimization approaches and therapeutic strategies.

Authors declare that they have no competing interests. Acknowledgment The authors express their gratitude to Professor Franco Maceri for fruitful discussions on this paper and to Simone Michele for useful support in some numerical tasks. This work was developed within the framework of Lagrange Laboratory, a European French–Italian research group. Present research study was partially supported by Italian Minister of University and Research, MIUR (PRIN, grant number F11J12000210001). Appendix A. Reduction to the Bessel equation Eq. (14) combined with the solution form in Eq. (17) lead to spacedependent differential equations, whose general structure is:

dy d2 y r2 2 + η r + ω r2 y = 0, dr dr

(A.1)

where y = y(r), and η and ω are real constants. By substituting the solution form y(r) = rs g(r) (with s ∈ R) and by assuming s = (1 − η)/2, Eq. (A.1) becomes

dg d2 g +r + (ω r2 − ν 2 ) g = 0, dr dr2

(A.2)

where ν = (1 − η)/2. By performing the variable change z = δ r, with |ω|, Eq. (A.2) reduces to the classic Bessel equation [56]:

δ=

d2 g 1 dg + + z dz dz2

1−

ν2 z2

g = 0.

(A.3)

When ν is not an integer number (as for the case under investigation), the solution of Eq. (A.3) is:

g(r) =

y(r) =

rs [A Jν (δ r) + B J−ν (δ r)],

if

rs [A Iν (δ r) + B I−ν (δ r)],

if

with s = ν = (1 − η)/2 and δ =

ω>0 ω≤0

(A.5)

|ω|.

Appendix B. Coeﬃcients Ck and orthogonality of eigenfunctions in Pχ

Competing interests

r2

where A and B are integration constants, and Jq (r) and Iq (r) indicate classic and modiﬁed ﬁrst-kind Bessel functions, respectively, of order q ∈ R. Thereby, the general solution of Eq. (A.1) results in:

A Jν (δ r) + B J−ν (δ r),

ω>0 A Iν (δ r) + B I−ν (δ r), if ω ≤ 0 if

(A.4)

Coeﬃcients Ck describing the closed-form solution of the problem Pχ are determined by combining the initial conditions (5) and the orthogonality property for the piece-wise eigenfunctions Rk (r) deﬁned in terms of Rki (r) (with i = 1, . . . , 4). It is immediate to prove that the following identity is satisﬁed for the i-th layer:

ri

ri−1

n n m Rm i Li [Ri ] − Ri Li [Ri ] dr = 0,

(B.1)

where Li denotes the spatial differential operator for the governing equation that applies in ri−1 < r ≤ ri . In detail, by combining Eqs. (14), (17), and (18) and by imposing T1 (t ) = Ti (t ) (for i = 2, 3, 4), Li results from:

Li [R(r)] =

d dR(r) qi (r) dr dr

+ δi2 qi (r) R(r),

(B.2)

where

q1 (r) = q4 (r) = r, qi (r) = r

Di −αi Di

,

(B.3)

for i = 2, 3,

and quantities α i (with i = 2, 3) and δi = δi (λ) (with i = 1, . . . , 4) deﬁned in Eqs. (16) and (20), respectively. Accordingly, Eq. (B.1) leads (for each i = 1, . . . , 4) to the following relationship:

D1 2 (λ − λ2m ) Di n

ri

n qi Rm i Ri dr =

ri−1

qi Rni

dRm dRni i − Rm i dr dr

ri .

(B.4)

ri−1

Let the quantity Imn be deﬁned as:

Imn =

4 i=1

Qi

k1 ε1 ki εi

ri

n qi Rm i Ri dr,

ri−1

(B.5)

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

with

the differential operators introduced in Eq. (B.2) reduce to:

Q1 = Q4 =

α2 −D2

1 , r1

Q2 = r1

D2

,

Q3 =

r α2D−D2 r α3D−D3 2 3 1 2 ·

r2

r3

·

r

α2 −D2 D2

1

α3 −D3

· r2

r2

D3

I

= 0,

for m = n

= 0,

for m = n

1 . r3

(B.6)

.

(B.7)

Due to the initial conditions expressed by Eq. (22), the following relationship holds: 4 i=1

k1 ε1 Qi ki εi

ri

C0 qi Rm i ci (r, 0)dr = r1 ri−1

Ck =

C0 r1 Ikk

r1

r0

r0

rRm 1 dr.

(B.8)

(B.9)

The method of separation of variables [51] can be employed to solve the problem P0 , to which Pχ reduces when curvature effects are disregarded. The solution has the same form as in Eq. (21), with the functions Rki (r) (with i = 1, . . . , 4), that deﬁne the piece-wise eigenfunctions Rk (r) for the corresponding Sturm–Liouville problem, expressed as:

Rk1 (r) = Ak1 sin (λk r) + Bk1 cos (λk r),

⎧ k ak r k bk r ¯2 if λk ≤ λ ⎨ A2 e 2 + B2 e 2 , k α2 R2 (r) = ⎩ e 2D2 r Ak2 cos δˆ2k r + Bk2 sin δˆ2k r , if λk > λ¯ 2 ⎧ k ak r k bk r ¯3 if λk ≤ λ ⎨ A3 e 3 + B3 e 3 , k α3 R3 (r) = ⎩ e 2D3 r Ak3 cos δˆ3k r + Bk3 sin δˆ3k r , if λk > λ¯ 3 + Bk4 cos

(C.1)

D1 k λr , D4

with λk the corresponding eigenvalue and with:

δ = ˆk 2

akj , bkj =

λ¯ 2 =

4D2 [(λk )2 D1 − β ] − α22

αj 2D j

2D2 ± i δˆ kj ,

α22 + 4D2 β 4D1 D2

,

δ = ˆk 3

α j = u¯ λ¯ 3 =

4D1 D3 (λk )2 − α32

γj

k jε j

,

(C.3)

d2 R(r) D1 2 + λ R(r), D4 dr2

with λ = λ1 . Therefore, constants Ck are univocally identiﬁed as:

Ck =

C0 Ikk

r1 r0

Rk1 (r) dr

(C.4)

q1 (r) = q4 (r) = 1,

qi (r) = e

α

− Di r i

, with i = 2, 3,

(C.5)

and

Q1 =

1 , r1

Q4 = e

α

α2

r

Q2 = e D2 1 ,

α

α

[ D3 r2 + D2 (r1 −r2 )] 3

2

, (C.6)

α

[ D3 (r2 −r3 )+ D2 (r1 −r2 )] 3

Q3 = e

2

.

References

Appendix C. Solution of problem P0

D1 k λr D4

L4 [R(r)] =

r1

r Rk1 dr.

Rk4 (r) = Ak4 sin

d2 R(r) α3 dR(r) D1 2 L3 [R(r)] = + − λ R(r), D3 dr D3 dr2

Ikk being formally deﬁned as in Eq. (B.5), with:

Finally, by combining Eqs. (B.7) and (B.8), coeﬃcients Ck (inﬁnite in number) result in:

d2 R(r) + λ2 R(r), dr2 d2 R(r) α2 dR(r) D1 2 β + L2 [R(r)] = − λ − R(r), D2 dr D2 D2 dr2 L1 [R(r)] =

,

Therefore, by using the boundary conditions (15), the orthogonality property for the piece-wise eigenfunctions Rk (r) of the Sturm– Liouville problem associated to Pχ can be expressed as: mn

95

2D3 ,

α32 4D1 D3

j = 2, 3,

, (C.2)

,

i being the imaginary unit and the advection coeﬃcients α j ( j = 2, 3) being now expressed in terms of the uniform ﬁltration velocity u¯ (see Eq. (12)). By using same arguments as in Appendix B, the values of the unknown coeﬃcients Ck in the solution form (21) specialized to P0 are obtained by combining the initial conditions (22) and the orthogonality property of the piece-wise eigenfunctions Rk (r). In this case,

[1] A.C. Guyton, J.E. Hall, Textbook of Medical Physiology, 12th ed., Saunders/Elsevier, Philadelphia, 2011. [2] B.F. Waller, C.M. Orr, J.D. Slack, C.A. Pinkerton, J.V. Tasselm, T. Peters, Anatomy, histology, and pathology of coronary arteries: a review relevant to new interventional and imaging techniques-Part I, Clin. Cardiol. 15 (1992) 451–457. [3] B.M. O’Connell, T.M. McGloughlin, M.T. Walsh, Factors that affect mass transport from drug eluting stents into the artery wall, Biomed. Eng. Online 9 (2010) 15. [4] R.E. Kuntz, D.S. Baim, Deﬁning coronary restenosis. Newer clinical and angiographic paradigms, Circulation 88 (3) (1993) 1310–1323. [5] J.R. Weiser, W.M. Saltzman, Controlled release for local delivery of drugs: barriers and models, J. Control. Release 190 (2014) 664–673. [6] C. Yang, H.M. Burt, Drug-eluting stents: factors governing local pharmacokinetics, Adv. Drug Deliv. Rev. 58 (2006) 402–411. [7] J.E. Sousa, M.A. Costa, A.C. Abizaid, B.J. Rensing, A.S. Abizaid, L.F. Tanajura, K. Kozuma, G. van Langenhove, A.G.M.R. Sousa, R. Falotico, J. Jaeger, J.J. Popma, P.W. Serruys, Sustained suppression of neointimal proliferation by sirolimus-eluting stents one-year angiographic and intravascular ultrasound follow-up, Circulation 104 (17) (2001) 2007–2011. [8] E. Ube, S. Silber, K.E. Hauptmann, Six and twelve month results from a randomized, double blind trial on a slow release pacilaxel eluting stent for de novo coronary lesions, Circulation 107 (2003) 38–42. [9] G. Nakazawa, A.V. Finn, E. Ladich, F. Ribichini, L. Coleman, F.D. Kolodgie, R. Virmani, Drug-eluting stent safety: ﬁndings from preclinical studies, Expert Rev. Cardiovasc. Ther. 6 (10) (2008) 1379–1391. [10] Y. Huang, S.S. Venkatraman, F.Y.C. Boey, E.M. Lahti, P.R. Umashankar, M. Mohanty, S. Arumugam, L. Khanolkar, S. Vaishnav, In vitro and in vivo performance of a dual drug-eluting stent (DDES), Biomaterials 31 (15) (2010) 4382–4391. [11] A. Schwarzmaier-D’Assie, N. Nyolczas, R. Hemetsberger, C. Strehblow, J. Matiasek, S. Farhan, Z. Petrasi, K. Huber, J. Wojta, D. Glogar, C. Plass, M. Gyöngyösi, R. Karnik, Comparison of short-and long-term results of drug-eluting vs. bare metal stenting in the porcine internal carotid artery, J. Endovasc. Ther. 18 (4) (2011) 547–558. [12] X. Ma, S. Oyamada, F. Gao, T. Wu, M.P. Robich, H. Wu, X. Wang, B. Buchholz, S. McCarthy, Z. Gu, C.F. Bianchi, F.W. Sellke, R. Laham, Paclitaxel/sirolimus combination coated drug-eluting stent: in vitro and in vivo drug release studies, J. Pharm. Biomed. Anal. 54 (4) (2011) 807–811. [13] W. Khan, S. Farah, A. Nyska, A.J. Domb, Carrier free rapamycin loaded drug eluting stent: In vitro and in vivo evaluation, J. Control. Release 168 (1) (2013) 70–76. [14] A. Seidlitz, S. Nagel, B. Semmling, K. Sternberg, H. Kroemer, W. Weitschies, In vitro dissolution testing of drug-eluting stents, Curr. Pharm. Biotechnol. 14 (1) (2013) 67–75. [15] B. Semmling, S. Nagel, K. Sternberg, W. Weitschies, A. Seidlitz, Impact of different tissue-simulating hydrogel compartments on in vitro release and distribution from drug-eluting stents, Eur. J. Pharm. Biopharm. 87 (3) (2014) 570–578. [16] A. Habib, A.V. Finn, Endothelialization of drug eluting stents and its impact on dual anti-platelet therapy duration, Pharmacol. Res. 93 (2015) 22–27. [17] Y. Liu, L. Gao, Y. Song, L. Chen, Q. Xue, J. Tian, Y. Wang, Y. Chen, Eﬃcacy and safety of limus-eluting versus paclitaxel-eluting coronary artery stents in patients with diabetes mellitus: A meta-analysis, Int. J. Cardiol. 184 (2015) 680–691. [18] M.A. Lovich, E.R. Edelman, Computational simulations of local vascular heparin deposition and distribution, Am J. Physiol. Heart Circ. 271 (1996) H2014–H2024.

96

M. d’Errico et al. / Mathematical Biosciences 267 (2015) 79–96

[19] S. Costantini, F. Maceri, G. Vairo, Un modello del rilascio di farmaco in stent coronarici (in Italian), Proceedings of the XVII National Congress of Computational Mechanics Group (GIMC), 2008. [20] C.W. Hwang, D. Wu, E.R. Edelman, Physiological transport forces govern drug distribution for stent-based delivery, Circulation 104 (2001) 600–605. [21] P. Zunino, Multidimensional pharmacokinetic models applied to the design of drug-eluting stents, Cardiovasc. Eng. 4 (2004) 181–191. [22] M. Grassi, G. Pontrelli, L. Teresi, G. Grassi, L. Comel, A. Ferluga, L. Galasso, Novel design of drug delivery in stented arteries: a numerical comparative study, Math. Biosci. Eng. 6 (3) (2009) 493–508. [23] G. Vairo, M. Cioﬃ, R. Cottone, G. Dubini, F. Migliavacca, Drug release from coronary eluting stents: a multidomain approach, J. Biomech. 43 (2010) 1580–1589. [24] X. Zhu, D.W. Pack, R.D. Braatz, Modelling intravascular delivery from drug-eluting stents with biodurable coating: investigation of anisotropic vascular drug diffusivity and arterial drug distribution, Comput. Methods Biomech. Biomed. Engin. 17 (3) (2014) 187–198. [25] R. Mongrain, I. Faik, R.L. Leask, J. Rodés-Cabau, E. Larose, O.F. Bertrand, Effects of diffusion coeﬃcients and struts apposition using numerical simulations for drug eluting coronary stents, J. Biomech. Eng. 129 (2007) 7333–7342. [26] D.R. Hose, A.J. Narracott, B. Griﬃths, S. Mahmood, J. Gunn, D. Sweeney, P.V. Lawford, A thermal analogy for modeling drug elution from cardiovascular stents, Comput. Methods Biomech. Biomed. Engin. 7 (2004) 257–264. [27] F. Migliavacca, F. Gervaso, M. Prosi, P. Zunino, S. Minisini, L. Formaggia, G. Dubini, Expansion and drug elution model of a coronary stent, Comput. Methods Biomech. Biomed. Engin. 10 (2007) 63–73. [28] E. Cutrì, P. Zunino, S. Morlacchi, C. Chiastra, F. Migliavacca, Drug delivery patterns for different stenting techniques in coronary bifurcations: a comparative computational study, Biomech. Model. Mechanobiol. 12 (2013) 657–669. [29] P. Zunino, C. D’Angelo, L. Petrini, C. Vergara, C. Capelli, F. Migliavacca, Numerical simulation of drug eluting coronary stents: mechanics, ﬂuid dynamics and drug release, Comput. Methods Appl. Mech. Eng. 198 (2009) 3633–3644. [30] L. Cattaneo, C. Chiastra, E. Cutrì, F. Migliavacca, S. Morlacchi, P. Zunino, An immersed boundary method for drug release applied to drug eluting stents dedicated to arterial bifurcations, Numer. Math. Adv. Appl. 198 (2013) 401–409. [31] C.J. Creel, M.A. Lovich, E.R. Edelman, Arterial paclitaxel distribution and deposition, Circ. Res. 86 (2000) 879–884. [32] I. Sheiban, G.P. Ballari, C. Moretti, W.G. Marra, E. Meliga, P. Omedè, F. Sciuto, G.P. Trevi, Paclitaxel-eluting stents for the treatment of complex coronary lesions: immediate and 12-month results, J. Cardiovas. Med. 8 (2007) 582–588. [33] G. Pontrelli, F. de Monte, Mass diffusion through two-layer porous media: an application to the drug-eluting stent, Int. J. Heat Mass Transf. 50 (2007) 3658–3669. [34] G. Pontrelli, A.D. Mascio, F. de Monte, Local mass non-equilibrium dynamics in multi-layered porous media: application to drug eluting stent, Int. J. Heat Mass Transf. 66 (2013) 844–854. [35] G. Pontrelli, F. de Monte, A multi-layer porous wall model for coronary drugeluting stents, Int. J. Heat Mass Transf. 53 (2010) 3629–3637.

[36] G. Pontrelli, F. de Monte, Modeling of mass dynamics in arterial drug-eluting stents, J. Porous Media 12 (2009) 19–28. [37] F. Bozsak, J.M. Chomaz, A.I. Barakat, G. Pontrelli, On the role of phase change in modelling drug-eluting stents, Biomedical Technology, Lecture Notes in Applied and Computational Mechanics, 74, Springer, 2015, pp. 69–80. [38] P. Macheras, A. Iliadis, Modeling in Biopharmaceutics, Pharmacokinetics, and Pharmacodynamics. Homogeneous and Heterogeneous approaches, Springer, New York, 2006. [39] G.A. Truskey, F. Yuan, D.F. Katz, Transport phenomena in biological system, 2nd ed., Pearson Prent. Hall, Upper Saddle River, 2010. [40] F. Maceri, M. Marino, G. Vairo, Age-dependent arterial mechanics via a multiscale elastic approach, Int. J. Comput. Meth. Eng. Sci. Mech. 14 (2013) 141–151. [41] G. Meyer, R. Merval, A. Tedgui, Effects of pressure stretch and convection on lowdensity lipoprotein and albumin uptake in the rabbit aortic wall, Circ. Res. 79 (1996) 532–540. [42] J. Bear, Y. Bachmat, Introduction to Modeling of Transport Phenomena in Porous Media, Kluwer Academic Publishers, Dordrecht, 1990. [43] A.-R.A. Khaled, K. Vafai, The role of porous media in modeling ﬂow and heat transfer in biological tissues, Int. J. Heat Mass Transf. 46 (2003) 4989–5003. [44] H. Park, E. Yeom, S.J. Seo, J.H. Lim, S.J. Lee, Measurement of real pulsatile blood ﬂow using X-ray PIV technique with CO2 microbubbles, Sci. Rep. 5 (2015) 8840. [45] N. Yang, K. Vafai, Low-density lipoprotein (LDL) transport in an artery - a simpliﬁed analytical solution, Int. J. Heat Mass Transf. 51 (3-4) (2008) 497–505. [46] O. Kedem, A. Katchalsky, Thermodynamic analysis of the permeability of biological membranes to non-electrolytes, Biochim. Biophys. Acta 27 (1958) 229–246. [47] M. Kargol, A more general form of Kedem and Katchalsky’s practical equations, J. Biol. Phys. 22 (1996) 15–26. [48] M. Kargol, A. Kargol, Mechanistic equations for membrane substance transport and their identity with Kedem–Katchalsky equations, Biophys. Chem. 103 (2003) 117–127. [49] M. Jarzynska, M. Pietruszka, Derivation of practical Kedem–Katchalsky equations for membrane substance transport, Old New Concepts Phys. 5 (2008) 459–474. [50] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, Springer, New York, 1999. [51] W.E. Boyce, R.C. DiPrima, Elementary Differential Equations and Boundary Value Problems, Wiley & Sons Inc., New York, 1969. [52] H.D. Baehr, K. Stephan, Heat and Mass Transfer, 2nd ed., Springer, Berlin, 2006. [53] L. Ai, K. Vafai, A coupling model for macromolecules transport in a stenosed arterial wall, Int. J. Heat Mass Transf. 49 (2006) 1568–1591. [54] M.A. Lovich, M. Philbrook, S. Sawyer, E. Weselcouch, E.R. Edelman, Arterial heparin deposition: role of diffusion, convection, and extravascular space, Am. J. Physiol. 275 (1998) H2236–H2242. [55] M. Shapiro, S. Hanon, D. Misra, Drug-eluting stents, Hosp. Physician 40 (8) (2004) 11–20. [56] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press Inc., New York, 1980.