Multiscale Boundary Conditions for Non-Fickian Diffusion Applied to Drug-Eluting Stents

. We propose a multiscale model for the study of blood ﬂow, plasma ﬁltration, and non-Fickian mass transport by a coronary drug-eluting stent. We derive an analytic solution of a one-dimensional model for non-Fickian diﬀusion and replace the problem of drug transport in the coating of the stent with suitable boundary conditions. In this way, the computational cost of numerical simulations in realistic three-dimensional stent geometries is reduced. Numerical experiments to show the eﬀectiveness of the method are also presented.


Introduction.
When an artery becomes obstructed due to the presence of an atherosclerotic plaque, it is necessary to perform an endovascular medical procedure called angioplasty.During this procedure, a medical device consisting of a small wire mesh tube, called a stent, is inserted and guided to the place of the obstruction with the aid of a catheter.Then, the catheter is removed and the stent stays, permanently acting as a supporting scaffold (Figure 1).
New generation drug-eluting stents (DES) are coated with a polymeric membrane containing antiproliferative drugs (paclitaxel, rapamycin, everolimus, etc.) in order to interrupt the cell cycle of smooth muscle cells (SMC) [31].In this way the proliferation of SMC can be regulated in situ and avoid the renarrowing of the lumen (restenosis) due to the accumulation of tissue around the stent.
In order to account for drug diffusion from the polymeric coating of the DES into the arterial wall, many mathematical models have been proposed [22].These models provide a low-cost alternative to research done in laboratories, where many experiments must be carried out to understand the different variables that can influence the process of the drug delivery, i.e., the geometrical design of the device, the type of polymer used in the coating, the type of drug embedded in the polymer, the initial loading of drug, etc.For brief reviews on models for drug deposition on the arterial wall and for three-dimensional models for drug transport in stents, we refer the reader to [1] and [17], respectively.Lower-dimensional models can provide some idea of drug kinetics in the DES and drug deposition in the arterial wall, but the fluid dynamics of the blood around the stent is fully three-dimensional [8,9,33]; therefore, in order to capture more accurately the interactions between drug diffusion and hemodynamics, it is important to set up mathematical models in three dimensions, considering realistic stent geometries.The computational cost of numerical simulations for a three-dimensional stent is very high, since it is necessary to consider many elements in order to generate a mesh that correctly represents the coating of the stent.Thus, a large number of degrees of freedom is needed to numerically solve the system of partial differential equations (PDEs).
Motivated by the fact that the polymeric coating of the DES is very thin, it is reasonable to approximate its surface by a uniform slab of finite thickness.Then, by calculating an analytical solution of the one-dimensional diffusion problem in the slab, it is possible to replace the drug diffusion problem in the coating of the stent with a Robyn-type boundary condition.Some multiscale models have been proposed in the literature [8,9,10,27,32], but to the best of our knowledge all of them neglect the influence of the viscoelastic properties of the polymer in the transport problem.The aim of this work is to deduce suitable boundary conditions for a mathematical multiscale model of drug delivery in a coronary DES on complex three-dimensional geometries that takes into consideration the viscoelastic properties of the polymeric coating of the stent.Therefore, our main focus is on DES with coatings where the polymer can affect the transport process.
In section 2 we introduce the complete model for the study of drug release in the coating of the DES, the lumen, and the arterial wall.A semianalytic solution of a one-dimensional problem for non-Fickian drug diffusion is described in section 3, while in section 4 a multiscale boundary condition to approximate the outgoing flux of drug into the lumen and the wall is presented.In section 5 we introduce the model for blood flow and plasma filtration.The complete multiscale model for drug transport is introduced in section 6.In section 7 we present the numerical scheme used to  approximate the solution of the multiscale model.Numerical results are presented in section 8, including some plots to exhibit the behavior of the mathematical model.Finally, section 9 is devoted to discussion and conclusions.

Modeling drug delivery.
A healthy arterial wall is a complex structure with different homogeneous layers (Figure 2).When atherosclerosis occurs, a lipid fibrous plaque is formed in the arterial wall due to the accumulation of a high level of low-density lipoprotein (LDL) cholesterol in the lumen and to a series of biochemical reactions that ends with the formation of the atheroma.This lipid accumulation causes the intima to thicken and change its physical properties, which may influence the drug deposition in the arterial wall [26].In this work, we consider as regions of interest the intima with the presence of a soft lipid plaque and the media, separated by the internal elastic lamina (IEL).For brief reviews on the modeling of LDL accumulation on the arterial wall, we refer the reader to [20,21].
Figure 3 shows a truncated portion of a stented coronary artery.In this domain we denote by Ω w = Ω in ∪ Ω p ∪ Ω me the arterial wall, where Ω in is the subregion of the wall that represents the intima layer, Ω p the atherosclerotic plaque, and Ω me the where Γ l,w is the interface between the lumen and the wall, Γ w is the artificial wall boundary, Γ c,w is the interface between the wall and the coating of the stent, Γ p,w is the interface between the plaque and the wall, Γ iel represents the IEL, and Γ adv is the interface between Ω w and the adventitia.
Analogously, the lumen boundary is defined as , where Γ c,l is the interface between the coating of the stent and the lumen, Γ in is the inflow boundary, and Γ out is the outflow.
The boundary of the polymeric coating of the stent is represented as ∂Ω c = Γ c,l ∪ Γ c,w ∪ Γ str , where Γ str is the interface between the polymeric matrix and the stent strut.The exterior unit normal vectors on ∂Ω w , ∂Ω l , Ω p , and ∂Ω c are denoted, respectively, by η w , η l , η p , and η c .
We neglect the presence of the endothelium as a physical barrier at the interface between the lumen and the wall.This is motivated by the fact that after the implantation of the stent the endothelium is injured and loses its ability to regulate the transport of molecules.
In order to model the drug diffusion problem we make the following assumptions: 1.The coating of the DES is made of a porous polymeric matrix that encapsulates a therapeutic drug.2. The initial concentration of the drug exceeds the drug solubility; therefore the entire drug is present from the start in the dissolved state.3. The transport of drug out of the polymer occurs by non-Fickian diffusion.4. Drug transfer does not affect blood and plasma flows. 5.The transport of drug in the tissue occurs by diffusion, advection, and saturable binding.In most polymers the plasticizing effect produced by a penetrant causes a change in its internal structure that may affect the transport process [6].In this case, many authors agree [6,7,11,16] that the classical Fick's law does not include this phenomenon.Therefore, a modified flux must be considered, given by where c denotes the concentration of the penetrant, J In this work we consider that the non-Fickian diffusion coefficient D v is constant.But, as shown in [13], a functional expression in terms of the concentration can be obtained for D v when the non-Fickian flux is interpreted as a convective field.
In order to study the viscoelastic behavior of the polymeric coating of the DES in the drug desorption process, we consider a mechanical model (Figure 4) which consists of a Maxwell fluid model in parallel with a free spring.The relaxation modulus can be calculated as [4,Chapter 5] where E 1 and E 0 denote the Young modulus of the spring elements of the Maxwell fluid arm and the free spring, respectively.The parameter τ 1 = μ1 E1 is the relaxation time of the Maxwell fluid arm, where μ 1 represents the viscosity of the dashpot element.From a practical point of view, this simple model is well adapted to simulating a realistic polymer behavior by fitting the parameters to the values obtained from experimental tests.
Let c c denote the concentration of drug in the coating.The stress associated to the desorption and exerted by the polymer is defined as where α > 0 is the proportionality constant between the strain of the polymer and the concentration, i.e, ε(t) = αc c (t).For a more complex model for non-Fickian transport in a stent, we refer the reader to [17].
The evolution of drug diffusion in the polymeric matrix is described by the following equation: where D c is the constant diffusion coefficient of the drug.
The concentration of the drug in the lumen, denoted by c l , can be tracked by the following equation: where D l denotes the diffusion coefficient of the drug in the lumen and u l is the velocity vector field in the lumen.It is known that drug binding to suitable sites of the extracellular matrix may strongly influence the deposition of the drug in the arterial wall [1,3,18,23,24,29].Therefore, in order to include this phenomenon in the model, we consider a nonlinear saturable reversible binding model [29].Let b denote the concentration of drug which is bound to nonspecific general extracellular matrix sites, and c w the concentration of drug which is free to diffuse and is transported by plasma.Then, the equations for drug diffusion in the arterial wall are given by where D w denotes the diffusion coefficients of the drug in the wall; u w the velocity vector field in the wall; k f,w and k b,w , respectively, the forward and backward rate constants; and b max,w the local density of binding sites.
The concentration of free drug c p and bound drug b p in the plaque are described by the following system of equations: where D p , u p , k f,p , k b,p , and b max,p are the same parameters defined above but now applied in the plaque.
The previous system of equations is closed with initial conditions and boundary conditions As interface conditions, we consider the continuity of the drug concentrations and the continuity of the normal components of the fluxes, For the transport across the IEL we impose the continuity of the fluxes and allow for a possible concentration jump, where η IEL is the exterior unit normal to the IEL, c in the concentration of free drug in the intima, D in the diffusion coefficient in the intima, c me the concentration of free drug in the media, D me the diffusion coefficient in the media, and P IEL the permeability of the IEL.
Assuming that most of the drug delivery occurs in the direction normal to the stent coating, we can calculate an analytical solution for a one-dimensional diffusion problem in the coating and obtain a Robyn boundary condition to approximate the problem (2), (8) with interface conditions (21) and (22).In this way, we significantly reduce the computational cost of modeling the drug delivery process over realistic three-dimensional stent geometries, as will be shown in section 8.

Analytical solution of the one-dimensional problem.
In what follows we calculate an analytical solution for a general non-Fickian desorption problem.We use Laplace transforms along with the residue theorem, following an approach similar to that in [25].
Let us denote by c 1 the concentration of drug; then the evolution of non-Fickian drug diffusion in a polymeric membrane of thickness L can be described by where We close the system with initial conditions and boundary conditions where c out is a constant that represents the drug concentration of the external medium.
The solution of ( 25)-( 28) is detailed in Appendix A and is given by , (29) where The time variation of mass of drug inside the membrane can be calculated as then from (29) we get 4. Multiscale boundary conditions.We note that the concentration c 1 defined in the previous section does not exactly satisfy the problem for the concentration c c defined in section 2, but motivated by the fact that there is a big difference of the spatial scales between the coating of the stent and the release media, the time evolution of the diffusion problem for c l and c w is much slower than that for c c [10,32].Then it is reasonable to use the expression (29) to replace the subproblem of drug diffusion in the coating of the stent ((2), with initial condition (8) and boundary conditions ( 21) and ( 22)) with a Robyn boundary condition.
In order to couple the coating with the lumen and the wall, we use the continuity of the fluxes to replace the original problem for c c by the boundary conditions where (33) for i = l, w.We consider that c out is identified with c l on Γ c,l and with c w on Γ c,w .
From ( 29) and (33) the flux of outgoing drug of the stent can be calculated as for i = l, w.

5.
Modeling the blood flow.Since blood is considered an isotropic Newtonian and laminar flow in the coronary arteries where the drug-eluting stents are typically implanted [14,9], we will assume the steady Navier-Stokes equations to model blood flow in the lumen.In fact, as discussed, e.g., in [2,28], the relevance of non-Newtonian effects in midsize arteries is negligible.For plasma filtration through the wall, assuming that both the intima and the media layers are porous media with isotropic and uniform permeability, we assume the Darcy equation.
Let us denote by p i the pressure field in each subdomain i = l, w.Then, the model is defined by the following system of equations: where μ b and μ p are the dynamic viscosities of the blood and the plasma, respectively, ρ b the density of the blood, and κ w the intrinsic permeability of the arterial wall.The previous system of equations is closed with boundary conditions where u in : R 2 → R 2 is a given function and p 0 ∈ R is a constant.
As interface conditions, we consider the continuity of the normal components of the stresses and the continuity of the normal components of the velocity, and the tangential component of the velocity in the lumen is set to zero, In order to numerically approximate the solution of problem ( 35)-( 46), we employ a finite element method in combination with a Picard's iteration to linearize the nonlinear advective term in (35).We approximate the plasma filtration problem in Ω w in terms of the pressure p w by quadratic Lagrange elements, recovering a posteriori the velocity field u w .The problem of blood flow is approximated with quadratic and linear Lagrange elements for u l and p l , respectively.

6.
A multiscale model for drug delivery.Before presenting the numerical scheme for the drug delivery problem, we briefly outline the model for drug release in the coating and drug transport in the lumen and the wall.The multiscale model is defined by and boundary conditions 7. Overview of the numerical scheme.In order to numerically approximate the solution of problem (47)-( 66), we propose a coupled implicit-explicit (IMEX) finite element method, since the problem of drug transport in the blood flow is advection dominated.We combine the IMEX method with the SUPG-method introduced in [5].
We will assume that each subdomain Ω i is a convex polygonal domain, equipped with a family of quasi-uniform triangulations T h,i made of affine simplexes K that are conforming on Γ l,w , Γ c,l , Γ p,w , and Γ c,w .
For the numerical solution of the diffusion problems in the lumen, in the wall, and in the plaque we define the following bilinear forms: The mass transport problems in Ω l , Ω w , and Ω p can be summarized as: find at each time step the functions and We note that in the IMEX method (67)-(69) we use the information from the previous iteration as boundary data.Therefore, it is possible to decouple the solution of the global discrete problem and separately solve each local subproblems on Ω l , Ω w , and Ω p .

Numerical results.
In what follows we consider different idealized stent geometries (Figure 5) for our numerical simulations.We approximated each geometry using the information available in the literature, which do not represent faithful copies of the original models.Figure 5(a) corresponds to a generic mesh stent, while Figure 5(b) is an approximation of the Promus element stent, manufactured by Boston Scientific Corporation.The Dynalink-E stent, manufactured by Abbott Vascular, is represented in Figure 5(c), and the Nobori stent, manufactured by the Terumo Corporation, is depicted in Figure 5(d).
In each case we consider a 2×14 mm stent with rectangular struts of size 0.15×0.19mm and a coating thickness of 0.014 mm.The embedment of the stent struts into the arterial wall is assumed at 50% for all the simulations.
The viscosity and the density of the blood are set, respectively, to μ b = 3.5 × 10    For the arterial wall domain we consider a typical left anterior descending coronary artery.We set the thickness of the healthy part of the intima to 0.2 mm, and the diseased part to 0.65 mm [15].Since the media is not affected significantly by the disease process, we set its thickness to 0.35 mm for both the healthy and the diseased parts [19].The thickness of the plaque is set to 0.15 mm (see Figure 6).
In this study, the release of the drug Sirolimus is considered.The drug diffusion coefficients in the media, lumen, and polymer are set to D me = 2.5 × 10 −10 m 2 /s, D l = 1.5 × 10 −9 m 2 /s, and D c = 1.2 × 10 −16 m 2 /s, respectively [23,24].The diffusion in the intima is in general approximated to be one hundred times faster than that in the media [19], and therefore we set the drug diffusion coefficient in the intima to D in = 2.5 × 10 −8 m 2 /s.The diffusion in the plaque can be approximately 85 times faster than in the media [19]; therefore, we set the diffusion coefficient of free drug in the plaque to D p = 2.94 × 10 −12 m 2 /s.
The permeability of the IEL is set to P IEL = 1 × 10 −9 m/s, as estimated in [19] from the transport of albumin through the IEL.The initial concentration of drug loaded in the stent is set to c 0 = 1000 mol/m 3 .
By M ci , we denote total mass of drug released at time t for the concentration c i , defined as In the case of the mesh stent (Figure 5(a)) we were able to generate an appropriate computational mesh for the coating and solve the original model described in section 2. In Figure 7 we present the time variation of the mass of bound drug in the media layer for the original complete model and the proposed multiscale model, normalized with respect to the maximum value of the mass for the complete model.We observe that the multiscale model is a good approximation of the complete model, with a maximum underestimation of approximately 15% for the first three days of the simulation.After this stage, multiscale and complete models produce similar qualitative results, with a maximum overestimation of approximately 10%.
The next results were obtained using the proposed multiscale model.Figure 8 shows the concentration of free drug in the arterial wall after 1-4 days, for the Dynalink-E stent (Figure 5(c)).According to this figure, we can observe that during the first days of the release there is a higher accumulation of drug in the intima because the diffusion coefficient in the intima is higher than the diffusion coefficient in   the media, and due to the presence of the plaque and the IEL acting as a barrier to diffusion from one layer to the other.In order to investigate the time variation of the mass of bound drug for the Dynalink-E stent (Figure 5(c)), we have plotted in Figure 9 the results considering the intima and the media layers and the plaque, normalized with respect to the maximum value of the mass in the intima.We observe that in each case, first the mass increases rapidly to an upper bound and then it decreases asymptotically.We also note from Figure 9 that during the first hours of drug release the amount of mass  in the intima is higher than in the other areas.Also as expected, the amount of mass in the plaque is smaller than in the other areas for the complete time domain.
Figure 10 displays the mass of free drug in the media for the different stent models, normalized with respect to the maximum value of the mass for the Dynalink-E stent (Figure 5(c)).We can see from Figure 10 that each stent has a different quantitative release profile.These variations account for the fact that the physical barrier to the flow in each case generates a different velocity field in both the lumen and the wall, influencing the mass transport.
The effect of the D v on the time variation of the mass of free drug in the media for the Nobori stent (Figure 5(d)) is analyzed in Figure 11.In this figure, the mass is normalized with respect to the maximum value of the mass when D v = 9 × 10 −19 s.The case D v = 0 corresponds to Fickian diffusion.We observe that the mass is an increasing function of D v .We note also that an increase in D v amplifies the effect of the relaxing flux term J NF .
Figure 12 shows the time variation of the mass of free drug in the media as a function of E 0 , for the Nobori stent (Figure 5(d)), normalized with respect to the maximum value of the mass when E 0 = 50000 Pa.We observe that the mass is an increasing function of E 0 .Like in the previous case, an increase in E 0 amplifies the effect of the relaxing flux term J NF .
In order to study the influence of the plaque in the transport problem, in Figure 13 we show the time variation of the mass of bound drug in the media for three different cases, considering the mesh stent (Figure 5(a)), normalized with respect to the maximum value of the mass in the absence of the plaque in the intima.The case of the solid calcified plaque is represented by excluding from the computational domain of the arterial wall Ω w the domain of the plaque Ω p .We observe that the presence of the plaque is important for the study of the drug deposition in the arterial wall.In particular, we point out from Figure 13 that there is a considerable difference in the quantitative behavior during the time domain considered.
In the previous figure we note a slight difference of the drug accumulation between the solid and soft plaques.We believe that this is motivated by the fact that we are considering the same values for the parameters associated to binding in the media, the intima, and the plaque.It is known that a lipid accumulation can alter the microviscosity of the tissue and local binding properties [19].However, to the authors' knowledge, there is not enough information in the literature on how to estimate the parameters for binding in the plaque.
Finally, to supply further evidence concerning the time variation of the mass of bound drug in the media, we have investigated the influence of the forward rate constant k f,p ; see Figure 14.To obtain these results, the Boston stent was considered, and the mass was normalized with respect to the maximum value of the mass when k f,p = 0.1.We observe that the mass is a decreasing function of k f,p .We note that the changes in the parameter k f,p can influence the deposition of drug in the media.Therefore, some effort must be made in the future in order to understand how to estimate the parameters for binding in the plaque.

Conclusions.
In this paper we presented a multiscale model for non-Fickian drug delivery in a coronary DES.Assuming that most of the drug delivery occurs in the direction normal to the coating, we approximated the diffusion problem inside the coating by a Robin boundary condition.We observed numerically that the reduced model is a good approximation of the complete model.
The main advantage of our approach is that it allows the incorporation of experimental rheological information about the polymer-solvent system into the transport problem, and it will significantly reduce the computational cost of simulations over realistic three-dimensional stent geometries.Finally, we remark that the Robin boundary condition for the outgoing flux of drug deduced in this work is not restricted to stents; in fact, it can be relevant in other applications where a solvent is slowly released from a thin polymeric coating.Substituting the previous expression for C 1 (s) into (73) gives c1 (x, s) = c out cosh (λ(s)x) s cosh (λ(s)L) + c 0 (cosh (λ(s)L) − cosh (λ(s)x)) s cosh (λ(s)L) .(74) In order to calculate the inverse Laplace transform of (74) with the aid of the residue theorem, we begin by considering the first term of (74).There is a singularity at s = 0, and its residue can be calculated as where i represents the imaginary unit.Then, from (72) we get that the roots must satisfy the equation which leads to (30) and (31).The residues at s = s jn for j = 1, 2 can be calculated as Concerning the second term of (74), we note that the singularity at s = 0 is removable.For the singularities at s = s jn for j = 1, 2, the residues can be calculated as

F
(c) = −(D∇c) represents the Fickian part of the flux, J NF (σ) = −(D v ∇σ) represents the non-Fickian part of the flux, and σ stands for the stress.The parameters D and D v represent respectively the Fickian diffusion coefficient and the non-Fickian diffusion coefficient.

Fig. 7 .
Fig. 7. Time variation of the mass of bound drug in the media layer for the complete model versus the multiscale model.

Fig. 9 .
Fig. 9. Mass of bound drug in the arterial wall layers and in the plaque.

Fig. 10 .
Fig. 10.Mass of free drug in the media for different stent geometries.

Fig. 11 .
Fig. 11.Mass of free drug in the media layer as a function of Dv.

Fig. 12 .Fig. 13 .
Fig. 12. Mass of free drug in the media layer as a function of E 0 .

Fig. 14 .
Fig. 14.Mass of bound drug in the media layer as a function of k f,p .