An intricate interplay between stent drug dose and release rate dictates arterial restenosis

and a


Introduction
Cardiovascular diseases are the leading cause of mortality in the world, with the number of deaths increasing annually [1].Of these, coronary artery disease (CAD) is the most prevalent, associated with the development of atherosclerotic lesions, where multiple constituents including macrophages, lipids, and cholesterol build up over time.Consequently, blood flow becomes restricted due to the reduced lumen, necessitating treatment strategies that restore patency, alleviating the risk of adverse events that may result in death.
The treatment of CAD with endovascular stents is now common practice.High in-stent restenosis (ISR) rates were associated with the original bare metal stents (BMS), devoid of drug [2,3].Innovation led to the inclusion of anti-proliferative drugs, embedded within a polymeric coating on the surface of the stent, released locally to the injury site with the aim of preventing the excessive neointimal tissue growth that leads to ISR.Clinical evidence charted a marked improvement in ISR rates, reducing from 17 to 41% with BMS to below 10% with drug-eluting stents (DES) [4].However, with DES, a small but notable increase in stent thrombosis was observed in some studies [5], believed to be a consequence of delayed arterial healing [6].Thus, dual antiplatelet therapy (DAPT) duration was extended to provide protection against thrombosis, in some cases for several years [7].Moreover, various studies have implicated either the polymeric coating and/or the drug in delayed re-endothelialization of the stented vessel [8][9][10][11].As a result, further research explored improved designs including: thinner struts that may induce less damage to the endothelium [12,13], the use of more biocompatible coatings and alternative drugs [9,14], and the removal of the polymer coating entirely [7].These newer generations of DES typically demonstrated superiority over earlier designs with respect to target lesion revascularisation (TLR) and stent thrombosis [4].
Compared with BMS, the improvement in patient outcome following PCI with DES is profound.However, despite the aforementioned technological advancements, repeat revascularisation remains an issue, particularly so when considering more complex lesions [15,16].For example, the recent FAME 3 study failed to show non-inferiority of percutaneous coronary intervention (PCI) with a 3rd generation stent to coronary artery bypass grafting (CABG) with respect to the incidence of death, myocardial infarction, stroke, or repeat revascularisation after one year [17].Particularly, the reported repeat revascularisation rate in the PCI arm at 1-year was 5.9%, expected to increase annually at a rate of 3-5%.While there are many factors that potentially explain such stubbornly high repeat revascularisation rates, the drug delivery protocol is under-explored as a potential target for further optimisation.As such, the focus of this article is to present a novel computational modelling framework that investigates how the drug delivery strategy modulates arterial restenosis.
Pioneering work by Evans et al. [18] introduced complex autonoma (CxA) models that permitted the description of the main physical and biological elements associated with restenosis, initiating a series of models that explored an agent-based approach.Boyle et al. [19] presented a cell-centred model that encompassed multiple species within a two-dimensional (2D) setting, where the constituents followed a set of rules in response to injury and inflammation.The model was extended to three dimensions (3D), where different stent design considerations were investigated [20,21].Zahedmanesh et al. [25] adopted a similar approach, identifying relationships between stent deployment and arterial stress which influenced the restenotic response.Keshavarzian et al. [26] presented a highly coupled agent-based model (ABM), depicting the behaviour of many chemical species involved in restenosis.
Moreover, a series of 2D ABMs [27][28][29][30] were extended to 3D [31,32,44], where many significant processes involved in and contributing to restenosis were investigated.Initially, studies explored the impact of strut dimensions and deployment depth within a 2D setting [27] before examining the significance of the regeneration of the endothelial layer, and the consequence this may have on the heterogeneity of restenotic tissue [28,29].Particular emphasis was then placed on deployment depth, where results demonstrated that deeper vessel injury facilitated the migration of smooth muscle cells (SMCs) into the lumen more easily [30].These works culminated in a fully coupled model within a 3D environment for the first time, investigating the impact of stent design and re-endothelialization on restenosis [31].The results illustrated reasonable agreement with trends present in clinical data before the model was further extended to account for ECM synthesis [32,44].
In contrast to ABMs, the continuum approach describes the evolution of populations of cells and other constituents such as growth factors (GFs) through PDEs.For example, Escuer et al. [38] detailed species evolution and tissue growth following a finite element (FE) simulation of stent expansion.Damage, a consequence of stress-induced stent deployment, is presented as an artificial construct emulating inflammation that upregulates various species that control species evolution.Assessing various stent designs, a general agreement with trends present in clinical literature was demonstrated.Moreover, Boland et al. [39] employed a ghost mesh strategy to simulate restenosis, highlighting the importance of arterial stresses that lead to ISR.Such a strategy permits an easier coupling to the mechanical aspects of the problem, whilst simultaneously modelling the explicit behaviour of constituents that contribute to neointimal hyperplasia.
Various studies exist that detail stent expansion and subsequent tissue growth within the FE framework.To the best of our knowledge, Fereidoonnezhad et al. [42] presented the first kinematic growth model, where increased mass as a consequence of tissue growth following balloon angioplasty was explored.He et al. [41] implemented a similar growth model, including tissue damage through the Mullins effect within their kinematic growth model.Three different stent designs were simulated, where the authors focused on the significance of stent materials, expansion severity and the likelihood of ISR following the deployment of overlapping stents; with subsequent efforts extended to patient-specific geometries [40].Maes et al. [43] adapted an ABM within their homogenized constrained mixture model framework to assess the arterial response to balloon expansion.Gierig et al. [36] considered a highly coupled chemo-mechano-biological model of restenosis, where the change in mass of local constituents drives the growth of tissue at a macroscopic level.In light of these efforts [36,[40][41][42][43] that couple the biological and mechanical framework, the need for further experimental evidence of soft tissue mechanics is made clear by the authors, such that the predictive capacity of in silico models can be improved.
Models simulating computational fluid dynamics (CFD) following stent deployment have established a link between low wall shear stress (WSS < 0.5 Pa) and ISR outcome [34,[45][46][47].For example, Corti et al. [45] devised a multi-scale model framework applied to human femoral arteries, where WSS directly influenced the functionality of ECs, such that areas exposed to lower levels of WSS gave rise to increased renarrowing of the lumen [34,45,48].Chen et al. [49] adopted a similar framework but within coronary arteries and coupled WSS to a damage index through a linear interpolation which directly influenced the behaviour of constituents present in restenosis.In agreement with [33,34,45,48,50], greater neointimal thickness was observed in regions where WSS was low.Furthermore, Jansen et al. [51] devised a multiscale bio-chemo-mechanical model of intimal hyperplasia, where haemodynamic indices are coupled to restenosis in an idealised artery through an influx of growth factors.Although CFD is included within these models [34,45,49], the role of stent deployment on arterial healing is not considered, as is present in [19][20][21]25,38]; a feature ubiquitously linked to neointimal thickening [52,53].
In summary, there is no existing comprehensive model of ISR, with each model having its own set of assumptions and limitations.However, regardless of the model framework selected, each of the aforementioned models fail to adequately couple drug transport and tissue retention with cellular function; a key limitation that we address in this article.First, we present a brief overview of current drug transport models and their limitations before providing a short review of the existing models that consider the effect of stent-induced drug delivery on restenosis, highlighting their limitations.

Mathematical and computational modelling of DES
Over the past two decades, mathematical models that investigate stent drug elution kinetics and subsequent tissue retention have elucidated the significance of controlled release on patient outcome .From these studies, the importance of prolonged specific receptor saturation on patient outcome was highlighted, a key characteristic associated with drug efficacy [54,56,58].The current state-of-the-art models [60,61,75] employ a 2D-axisymmetric multilayer geometry inspired from those proposed in [71,72,76], but incorporating two phases of nonlinear saturable binding.These efforts primarily focused on drug transport and subsequent tissue retention, highlighting the role of stent expansion [60], vessel disease and curvature [61], and different endovascular devices (balloons versus stents) [75] on indicators associated with device safety (tissue drug content) and efficacy (specific receptor saturation).The significant limitation of these state-of-the-art models is the consideration of a static artery: they neglect the influence of specifically bound drug on cellular function, and are therefore unable to model the influence of drug on tissue growth.

Mathematical and computational modelling of arterial re-narrowing in response to stent-induced drug delivery
Given that DES are known to significantly reduce the occurrence of ISR, it is somewhat surprising that none of the aforementioned ISR models (Section 1.1) adequately incorporate state-of-the-art models of stent-based drug delivery, transport within arterial tissue (Section 1.2) and importantly, the effects of drug on SMC proliferation, the key biological process that contributes to restenosis.To the best of our knowledge, only a few models have attempted to consider the effect of drug on restenosis.However, these models adopted a significantly simpler drug transport model.Specifically, Caiazzo et al. [77] and Tahir et al. [27] were the first to incorporate the effect of drug within their agent-based modelling framework of ISR.Only free drug concentrations were considered, with binding kinetics and receptor saturation neglected, features that were subsequently shown to be of critical importance [54].In the following years, efforts shifted to a continuum modelling framework, where Rossi et al. [55] simulated the effect of drug on cell growth through a linear dependence coupled to an apoptotic term, where the drug was cytotoxic as opposed to cytostatic.Moreover, Peddle et al. [78] adopted a generic two-species model, emulating both SMC phenotypes.They assessed the role of drug as a growth inhibitor or phenotype switching blocker, where release kinetics are simulated through an exponential condition and drug transport neglected.More recently, Maes et al. [43] presented a preliminary approach to detail the inhibition of SMCs following the deployment of a drug-coated balloon (DCB).Without simulating drug transport explicitly, they examined how SMC phenotype switching may be mediated following drug release by adjusting the appropriate rate constant.This subsequently influenced the proliferative capacity of active (synthetic) SMCs where vascular remodelling through collagen deposition is also inhibited.
In this article, we address the aforementioned limitations by developing a comprehensive computational modelling framework that simulates stent deployment, drug delivery and subsequent drug transport and retention within a remodelling artery.Addressing the limitations present in [43,78] we consider a drug transport model that simulates the release of free drug from a polymer coating into the arterial wall, where it is able to bind both specifically and non-specifically, an important feature absent in [27,77].In line with current thinking (e.g.[54]), we assume that that the level of specific bound drug is associated with drugeffectiveness.The key novelty of our approach is the explicit coupling of spatiotemporal drug delivery and interference with the cellular processes that give rise to restenosis.We use the modelling framework to assess the impact of stent design, exploring a variety of different strut configurations including different thicknesses, embedment depths and inter-strut spacing.Moreover, to establish the significance of controlled release, several different drug release profiles and initial drug loadings are considered.

Model geometry
A 2D-axisymmetric geometry corresponding to an idealised segment of a coronary artery is considered for all simulations.A two-layer arterial wall hosting the intima and media is simulated (Fig. 1), with layer thicknesses obtained from the literature [38,60].The baseline model (Table 1) emulates current stent designs which consider an ultra-thin strut thickness (s t < 80 μm) coated with a polymer that has a generic thickness of 8 μm (δ p ). Struts are separated by a spacing (δ ss ) equivalent to 7×total strut thickness (T st ), accounting for the circumferential coating thickness, and are half-embedded within the arterial wall prior to stent expansion [60,73].
In the present study, fluid flow within the lumen is not considered, with justification for neglecting the influence of fluid dynamics on restenosis and drug transport provided in Sections 2.4 and 2.6, respectively.Taken together with the idealised cylindrical arterial configuration, this assumption allows further reduction in the complexity of the geometry and, consequently, computational time.Specifically, through symmetry, the model governing equations need only be solved on half of the arterial wall domain, with an axis of symmetry imposed between the middle struts (Fig. 1).Our baseline model therefore considers 2 equally spaced struts, with symmetry implying a stent totalling 4 struts.However, the number of struts is not particularly important, as highlighted in [38], where the tissue growth simulations clearly demonstrate identical patterns of tissue growth between struts, justifying our approach to simulate a 'shorter stent'.Moreover, we confirmed that symmetry was indeed observed by consideration of the full geometry, as well as configurations involving varying numbers of struts (4, 7 and 10).
Fig. 1a illustrates the arterial segment pre-expansion, whilst Fig. 1b shows the altered geometry following stent expansion emphasised by the difference in lumen radii (r l,h versus r l,ref ).In the baseline model, the stented domain covers a lesion length of 4 × δ ss .Extending this domain a further 3 × δ ss from the centre of the outermost strut completes the therapeutic domain, where indicators associated with drug efficacy and safety are computed [60,61].Further domain extensions are necessary such that boundary conditions upstream/downstream from the lesion can be justified [38,60].All domains and boundaries present in the model are denoted by Ω and Γ, respectively.

Computational model overview
The arterial healing process is immensely complex, encompassing a plethora of species which interact in response to a mechanical stimuli, aiming to heal the vessel in a timely manner.The model presented by Escuer et al. [38] is built upon here, with improvements made to various aspects.
Fig. 2 schematically summarises the model framework, hosting three distinct steps: stent expansion, species evolution, and remodelling.Initially, stent expansion deforms the artery where the new geometry (Fig. 1b) and Von Mises (σ vm ) stress map are exported and utilised within the species evolution model.High levels of stress trigger a cascade of biological events through a variable called damage which, in our model, contributes to the increase in various species that stimulate healing, specifically growth factors (GFs) and matrix-degrading metalloproteinases (MMPs).The former promotes cell growth explicitly, whilst the latter, through the degradation of extracellular matrix (ECM), induces a switch in SMC phenotype to an active state, where cells are free to migrate and proliferate.As SMCs proliferate, additional ECM is synthesised as part of the healing process, providing structural support to the vessel.Note that the location of drug interference is also highlighted, where a reduction in SMC proliferation is expected following drug release, which subsequently influences the healing process.Within the remodelling step, it is the SMCs and ECM which primarily comprise the restenotic tissue, alongside the regenerating endothelial layer.
Relative to Escuer et al. [38], amendments to the stress-strain model and key constituents (SMCs and ECM) within the species evolution model are considered.Particularly, the constitutive model of the artery is improved to account for collagen fibre dispersion within the arterial wall.As for the restenosis model, ECM synthesis is amended to prevent excess growth in low expansion scenarios as discussed in the Supplementary material (Section S3.1).In addition, SMC proliferation is altered such that growth is mediated by spatial elements and SMC density, accounting for more realistic dynamics present in vivo.

Stent expansion: mechanical model
The arterial model comprises of two (intima and media) unique layers each modelled as a hyperelastic, incompressible material.Accounting for the family of collagen fibres that reinforce the vessel wall, the strain energy function (SEF) proposed by Gasser et al. [82] and in a Longitudinal section view following stent expansion.The inset highlights the polymer coating, intima, and denuded endothelium.Ω j defines the various domains, where subscript j = p, m, i denotes the polymer coating, media, and intima, respectively.Boundaries are denoted by Γ.The healthy and denuded endothelium are defined by Γ et, h and Γ et,d , respectively.The polymer coating boundary is in contact with either the wall (Γ pw ) or lumen (Γ pl ), whilst Γ ps is the interface between the coating and the strut.The boundary between the intima and media is Γ iel whilst the media/adventitia interface is denoted Γ eel .Upstream from the stent, the boundaries for the intima and media are Γ i,up and Γ m,up , respectively.For full details of notation, the reader is referred to the main text.

A. McQueen et al.
recent study by [60] is adopted here: where I 1 = tr(C) refers to the first invariant of the Cauchy-Green deformation tensor, C, detailing the isotropic potential of noncoallagenous matrix [83].Invariants, I 4 and I 6 detail the mechanical response in the preferential directions of the fibres.Parameters k 2 and 0 < κ ≤ 1/3 are dimensionless and describe the level of collagen fibre dispersion, such that κ = 0 describes fibres that are perfectly aligned and κ = 1/3 refers to an isotropic material where fibres are randomly distributed.The remaining parameters, μ and k 1 are positive constants that have dimensions of stress.The anisotropic directions are assumed to be helically orientated at ±β degrees with respect to the longitudinal direction, such that I 4 and I 6 can be expressed as functions of the main stretches, as reported in [60], where the constitutive model is discussed in more detail.All constitutive model parameters were obtained from [60], which were calibrated against the experimental data provided in [84].Note, with all layers of tissue assumed incompressible, the volume ratio J = detF = 1, where F is the deformation gradient tensor.The metallic strut is modelled assuming a bilinear elastic-plastic relationship with a Young's modulus (E), density (ρ), Poisson's ratio (ν), yield strength (S y ), and an ultimate tensile strength (E t ) representative of cobalt chromium alloy L6505 [85].The stent coating is simulated as a phosphorylcholine (PC)-based polymer assumed to behave as a nonlinear elastic material with bilinear behaviour.Parameter values for the artery layers and stent materials are detailed in Table 2.
The strategy for stent expansion employed in [38,60] is adopted within this work.Firstly, the boundaries between the lumen and tissue, (Γ et,d and Γ et,h ), are pressurised at 100 mmHg to mimic physiological conditions.A displacement on the stent struts is subsequently prescribed to simulate stent deployment.Three clinically relevant levels of overexpansion (OE) are considered (10/20/30%): where r l,ref is the lumen radius within the stented domain, postexpansion, and r l,h is the lumen radius within the unstented domain, equivalent to that present pre-expansion.The resulting deformed arterial geometries and Von Mises stress maps, σ vm (r, z), are then used as inputs to the restenosis model.

Species evolution: restenosis model 2.4.1. Damage
As is common in the literature (e.g.[25,38]), we assume that damage induced by stent deployment is the predominant driver of the biological cascade of events that initiate the healing process.Alternative approaches in the literature include linking haemodynamic indices such as WSS to a damage index [49], a function describing endothelial dysfunction [45] or an influx of growth factors [51].Damage, d (r, z, t), is modelled as a spatiotemporal variable mapped throughout the domain in response to the stresses present following stent expansion.We impose   inferior (σ vm,inf ) and superior (σ vm,sup ) limits which correlate to minimum (0) and maximum (1) levels of damage, respectively, obtained from in vitro experimental tensile tests performed on coronary arteries [84].For intermediate levels of stress, in the absence of data to suggest otherwise, a linear interpolation is imposed between the two limits [38].The initial level of damage, d (t = 0) = d 0 , is defined by: We impose that damage subsequently decays continuously as a function of the MMP concentration: where k d , the rate of decay of damage, is linked directly to the current MMP concentration, c mmp (r, z, t).

Growth factors (GFs)
Upon injury, an array of GFs are produced, including, but not limited to platelet-derived GF (PDGF), epidermal GF (EGF), and transforming GF (TGF) which behave in harmony to stimulate healing [52].Our model encompasses each of these within a single species, c gf (r, z, t), justifiable since their specific roles in restenosis are still poorly understood: where c 0 gf is the initial concentration of GFs, assumed constant across the arterial wall.The movement of GFs within the domain is detailed through a diffusive term, with associated diffusion coefficient, D gf , assumed constant and isotropic.The rate of GF production is defined by g gf , which is upregulated by damage (d) to some upper limit (K gf ), whilst GF degradation to baseline levels (c 0 gf ) is dictated by k gf .

Smooth muscle cells (SMCs)
The heterogeneity in vascular SMC phenotype is heightened in the presence of vascular injury as the artery responds to damage by initiating the healing process.Various factors mediate this switch in phenotype, from ECM degradation to the exposure of cells to circulating GFs via a dysfunctional/damaged endothelium [25,52,87].The dormant, contractile (cSMC) phenotype is governed by Eqn.(2.6), whilst the active, synthetic (sSMC) phenotype that contributes significantly to restenosis through proliferation follows Eqn.(2.7).Their respective concentrations are denoted by c csmc (r, z, t) and c ssmc (r, z, t).
∂c ssmc ∂t where the proliferation of synthetic SMCs (P ssmc ) is governed by: The rate of differentiation from the contractile to synthetic phenotype is given by k diff csmc , whilst k diff ssmc denotes rate of differentiation back to contratile SMCs from synthetic.Following Escuer et al. [38], the functions K ecm ssmc (Eqn.(2.9)) and K ecm csmc (Eq.(2.10)) modulate the differentiation between species, directly related to the ECM concentration.A surge in synthetic cells is present following degradation in ECM (c ecm < c 0 ecm ), whilst the resurgence of the contractile phenotype occurs through the remodelling process as ECM synthesis occurs (c ecm ≥ c 0 ecm ): ) , (2.9) The migration of sSMCs is enabled through a diffusive term, with associated diffusion coefficient, D ssmc , assumed isotropic and constant.A significant improvement to the restenosis model presented in [38] relates to how sSMC proliferation is modelled.This process is driven by the growth rate, g ssmc , whilst being coupled directly to the current density of sSMCs (c ssmc ), GFs (c gf ) and logistically to the ECM population (c ecm ), such that more realistic spatial restrictions on growth are introduced.We allow sSMCs to die naturally through apoptosis, at a rate k ssmc .As for the contractile SMC population, we assume that the movement, growth and death are negligible.It is assumed that contractile SMCs exist only in the media initially at some constant density, c 0 csmc , while the initial density of synthetic SMCs is assumed to be zero in all domains.

Extracellular matrix (ECM)
The ECM, a non-cellular constituent which provides structural support for the cellular species residing within the arterial wall, is a crucial part of the arterial healing response, and as a consequence is often the dominant species present in restenotic tissue [88].In our model, collagen is considered to be the primary component of ECM [38].The equation governing ECM is given by: ∂c ecm ∂t where c ecm (r, z, t) is the ECM concentration, initially present in all domains at a constant concentration, c 0 ecm .It is synthesised at a rate g ecm by sSMCs and degraded at a rate k ecm through an increase in MMP production.ECM is modelled to synthesise to some upper limit (K 1 ecm ), whilst degrading to a lower limit of (K 2 ecm ) such that a concentration of zero is not possible within our model.Another significant change to the restenosis model presented in [38] is the definition of K 1 ecm .We define an expression to detail how the upper limit varies with respect to the initial damage in the stented area: where θ is the mean damage, d (r, z, t), within the stented domain at t = 0.This relationship, motivated through consideration of existing literature, dictates that increased levels of damage should produce more ECM as part of the remodelling process [89].

Matrix metalloproteinases (MMPs)
MMPs are a large family of proteins that exist to mediate the levels of ECM.Our model, in line with [38], specifically refers to MMP-2.The synthesis of MMPs is governed by both c csmc and c ssmc and upregulated by d, to initiate an inflammatory response that drives healing.The rate of change of MMPs is described by the following equation: where c mmp (r, z, t) is the MMP concentration, initially taking the constant value c 0 mmp within the arterial wall.Movement of MMPs is described by a diffusive term, with an associated diffusion coefficient, D mmp , that is both isotropic and constant.The synthesis of c mmp by synthetic and contractile SMCs are governed by g ssmc mmp and g csmc mmp , respectively, limited by a threshold condition of K mmp .The rate of MMP degradation is mediated by a linear rate constant, k mmp [38].

Endothelial cells (ECs)
Endothelial cell regrowth is the cornerstone of arterial healing.The presence of this layer of cells ensures natural order within the vessel, where a release of numerous chemical mediators (e.g.nitric oxide) diffuse through the wall to maintain quiescence [28,90,91].The ability to suppress tissue growth is halted following stent expansion causing endothelium denudation which stimulates the healing response through the aggregation of inflammatory cells.Similarly to Escuer et al. [38], EC regeneration is modelled as: where c ec (r, z, t) is the density of ECs, migrating as a result of diffusion, D ec , assumed isotropic and constant.Endothelial cell proliferation is driven by a constant growth rate, g ec , which approaches an upper limit, K ec , emulating a fully repaired endothelial layer.Within the stented domain, we assume that the concentration of ECs is approximately zero, mimicking endothelial denudation as a consequence of stent deployment.In the domains unaffected by stenting, the initial EC concentration (c 0 ec ) is equivalent to K ec .

Model parameters
The majority of model parameters associated with the restenosis model have previously been reported [38] and are presented in the Supplementary material (Table S1).However, due to amendments in our restenosis model, two new parameter values are introduced, g ssmc and K 1  ecm , with the former chosen such that the proliferative time scale is comparable with [38] and the latter described by equation (2.12).

Boundary conditions
Within the arterial domain, continuity of concentration and flux is assumed at the intima-media interface (Γ iel ) for all species that reside across the two domains (damage, GFs, SMCs, ECM, and MMPs) [38].Moreover, zero-flux ( − n⋅ ( − D j ∇c j ) = 0), where n is the unit normal vector, and j denotes each layer of arterial tissue, was imposed on all relevant boundaries: lumen-wall interface (Γ et,h and Γ et,d ), strut-wall interface (Γ pw ), external elastic lamina (Γ eel ), upstream (Γ j,up ) and at the axis of symmetry.In the case of Γ et,h , Γ et,d , Γ pw , and Γ eel , the application of zero-flux ensures that species are unable do leave the tissue domain, thereby contributing to tissue growth.The choice of zero-flux boundary conditions at Γ j,up is motivated by the fact that the upstream boundary is chosen to be sufficiently far from the therapeutic domain, where quantities of interest are calculated, meaning that this boundary condition has negligible impact on the results, in agreement with previous studies [60,67,71,72].At the axis of symmetry in the centre of the domain, the zero-flux condition ensures no concentration gradient across this axis.The ECs are exposed to zero-flux on all boundaries, including Γ iel to ensure that the species exists only within the intima.

Tissue growth: remodelling
A model for growth of biological tissue in a continuum framework [92] is applied here, adopted previously in [38,93,94].Growth is simulated as a stationary problem, where after 300 days, as proposed in [38], the transient simulation for species evolution is concluded and various species contribute to the volumetric growth of tissue where the balance of mass in the system must satisfy: where i indicates the species contributing to growth; ρ 0 i is the concentration of species i in the reference configuration; Π i defines source and sink terms of these species which encompass the processes of proliferation, production, apoptosis and degradation and; M i are mass fluxes of the species considered.
The total material density is thus the sum of the concentration of all individual species, such that ρ 0 = ∑ i ρ 0 i .These densities evolve in time as a consequence of equation (2.15), where growth implies an increase in total species concentration, whilst shrinking of the geometry corresponds to the opposing behaviour.Assuming these volumetric changes in species are locally isotropic, then F g i = ρ 0 i /ρ init i I can be defined as the growth deformation tensor, where ρ init i is the initial density of species i in the reference configuration and I is the isotropic tensor of the second order [92,93].Thus, under the small strain approximation, where growth is isotropic, then: where v i is the velocity of material points [38,93].Contributing to growth in our model are four species: ECs, ECM, and SMCs (contractile and synthetic).Thus, isotropic growth as a consequence of species evolution can be defined as: where Δc i , defines the change in species concentration.Note, smc accounts for the total SMC population, inclusive of both phenotypes.In this approach, it is the volume (V i ) of cellular components and the density (ρ i ) of non-cellular components that contribute to growth.The volume of EC cells is assumed spherical (Eq.(2.18)) [25,38] and SMCs ellipsoidal (Eq.(2.19)) [87]: ) where ECs have a characteristic radius, r ec and SMCs have a typical radius and length of r smc and l smc , respectively.

Drug transport model
The polymeric coating is presumed to be durable, where the release of drug is a diffusion driven process: where c p (r, z, t) is the free drug concentration in the polymer coating, with an effective constant diffusion coefficient of D p .At t = 0, drug is assumed to reside only within the polymer coating at some initial constant concentration, c 0 p = M 0 /V p , where M 0 is the initial mass of drug per strut and V p is the volume of the coating per strut.
The concentration of free drug within the arterial wall, c w (r, z, t), is described by a reaction-diffusion equation: where c w is initially zero, D w is the drug diffusivity tensor in the wall (Eq. (2.22)), noting that diffusion in the radial (D r ) and axial (D z ) directions differ: Note that in the present model the influence of advection on drug transport is neglected.This assumption is valid in cases where the arterial Peclet number (Pe) is less than 1, in which case diffusion dominates over advection.
Free drug depletes due to binding reactions within the arterial wall, with the concentration of specific bound drug and non-specific bound drug denoted by b s (r, z, t) and b ns (r, z, t), respectively.The rate of change of concentration of specific and non-specific bound drug are governed by: where k x on and k x off (x = s, ns) denote the binding on and off rates whilst B x c (x = s, ns) refers to the maximum local density of binding sites present, which in the baseline model is assumed constant.
We consider our model drug to be sirolimus, the most widely studied drug in the context of drug transport and retention following delivery from DES and for which a full set of model parameters is available in the literature [54,61,71].For completeness, these are displayed in the Supplementary material (Table S2).Further motivation for considering sirolimus lies in the fact that previous models which implicated saturation of specific receptors with sirolimus in DES efficacy ignored the cellular response to the drug [54,59].
A zero-flux condition is imposed at the interface between the coating and metallic strut (Γ ps ) for free drug within the coating (c p ).Moreover, upstream from the stent (Γ i,up , Γ m,up ) and at the axis of symmetry in the middle of the domain, zero-flux is imposed for free drug within the wall (c w ).Continuity of both concentration and flux is assumed at the interface between the coating and wall (Γ pw ), and between arterial layers (Γ iel ).Following [75], where it was shown that drug is rapidly washed away by the blood, leading to very low drug concentrations in the lumen, an infinite sink condition is imposed at the lumen coating interface (Γ pl ) and at the lumen-wall interface (Γ et,h , Γ et,d ).Finally, an infinite sink is also considered at the boundary between the media and adventitia (Γ eel ) [71].

Coupling drug release and transport to restenosis
Since DES efficacy has been associated with saturation of specific binding sites, assumed to exist within SMCs, our baseline model assumes that sSMC proliferation is reduced in proportion to the fraction of specific binding sites that are bound by drug.Thus, the sSMC proliferation term (P ssmc ) is amended to: where the cytostatic effect of the drug is considered such that cell proliferation is reduced as a consequence of specifically bound drug (b s ) approaching its upper limit (B s c ); denoted by the under-brace termed 'Drug Interference'.This occurs as free drug (c w ) infiltrates the wall and binds either specifically or non-specifically (b ns ) as a result of equations (2.23) and (2.24), respectively.Cell proliferation is then completely ceased when all specific receptors are bound with drug (b s = B s c ).This is in agreement with static drug transport models [54,56], where the saturation of specific receptors has been directly linked with the efficacy of DES.To be consistent with these models, the maximum binding site density (B s c ) is assumed constant in our baseline model.However, to assess the validity of this assumption, an alternative approach is considered (Section 2.7.1) that probes a dynamic, variable binding site density that varies in proportion to the sSMC density.
Following the recent findings of McQueen et al. [95], two alternative approaches to modelling the effect of drug are also considered.In each approach cell proliferation is reduced through a Michaelis-Menten term.The reader is referred to Section S1 in the supplementary detail for full details on the model equations of these alternative drug models.

Impact of binding site density
Sirolimus is known to exert its effect through binding to the FKBP12 receptor located within SMCs.This impacts the functionality of FKBPrapamycin-associated protein (FRAP), inhibiting the biological cascade of events stimulated as a result of vascular injury.Current models which incorporate drug transport assume a static arterial wall, in other words the SMC density, and therefore the density of FKBP12 binding sites, remains constant [54,59].Whether this assumption remains valid in a healing (dynamic) arterial wall, has yet to be fully resolved.Recently, we [95] assessed the impact of a variable binding site density within a standard in vitro cell proliferation assay, where SMCs were exposed to constant doses of drug.We now extend these efforts to consider the in vivo situation, enabling us to investigate the influence of spatiotemporal delivery of drug on tissue growth.Here, the heterogeneous nature of damage following stent deployment leads to spatiotemporal variation in sSMC proliferation.
In an animal study involving pigs, [54] noted a minimal change in receptor expression of FKBP12 following injury, apparently consistent with the assumption of a constant density of specific binding sites.Other efforts which explored this behaviour noted an increase in receptor expression shortly after vascular injury, which decreased in time as the neointimal tissue became more mature [96,97].To explore this issue further, we now assume that the specific receptor density varies with the sSMC population in a spatiotemporal manner.The following expression replaces the previous constant binding site density (B s c ) in equation (2.25): where R c is a non-dimensional quantity that captures the receptor density per cell.To overcome a singularity in equation (2.25), it is assumed a small density (10 − 3 ) of sSMCs exist initially.The value was chosen to correspond to a density an order of magnitude below that present in the drug-free simulations after 5 seconds.Simulating the behaviour over 50-300 days, analysis of this assumption displayed a negligible difference in the overall model response when values an order of magnitude above and below were considered.
In the absence of reliable estimates for R c , we assume: where X ssmc corresponds to one of three values associated with the spatially-averaged sSMC concentration across the therapeutic domain: (i) maximum in time (ii) time-averaged over 50 days or (iii) timeaveraged over 21 days.Detailed in Supplementary material, Fig. S4 is used to obtain the different X ssmc values for (i)-(iii), with the corresponding values for R c presented in Table S4.

Numerical methods
The commercially available software COMSOL Multiphysics 5.6a (COMSOL AB, Burlington, MA, USA) was used to create the finite element (FE) mesh and to numerically solve the computational model in three consecutive steps: (i) a stationary mechanical analysis emulating stent deployment to quantify tissue damage, (ii) a time-dependent highly nonlinear coupled PDE system depicting the evolution of species present in ISR within the deformed geometry from (i), and finally (iii) a further stationary mechanical step to depict how transient species evolution corresponds to tissue growth.A mesh independence study was considered, investigating the impact of mesh size on % area restenosis, and the ECM concentration.A mesh was selected such that perturbations in these indices were below 1%.This corresponded to an 'extremely fine' mesh in the absence of drug, where the domain size varied across the different strut configurations.Further refinement in mesh quality (addition of boundary layers) was necessary to ensure convergence in the coupled drug transport model.The domains of the model were discretized with a combination of triangular and quadrilateral elements, resulting in a mixed mesh ranging from 65,000 to 150,000 elements in size.A direct linear solver (MUMPS) was used to solve the stationary problems with a relative error tolerance of 10 − 3 .The implicit backward differentiation formula (BDF) method was used for the time discretization of the transient study step, consisting of the restenosis and drug transport models, with variable order of accuracy varying from one to five in order to obtain better stability and variable time step size.Further stability was ensured by enabling the nonlinear controller.The relative and absolute tolerances were set to 10 − 3 and 10 − 4 , respectively.Using 16 cores of an AMD Ryzen TM 9 5950X CPU @ 3.40 GHz processor, the computational time for each of the simulations carried out varied between 1.5 and 23 hours, as a consequence of perturbations in the geometrical configuration and model set-up.

Summary of investigated cases
A variety of different cases are considered, detailed in Table 3.Firstly, Study A considers the impact of stent geometry and expansion ratios in the absence of drug delivery.The baseline model considers three different levels of expansion (10, 20, 30%) with a fixed strut thickness (ST), level of stent strut embedment (SE) and strut spacing (SS).Then, ST, SS and SE are varied in turn for the highest expansion ratios.In practice, since we keep the domain size fixed, varying SS corresponds to different stent designs, where the strut number varies within the lesion.As such, two different cases are proposed: SS small and SS large .The former comprises of 7 struts within a fixed lesion of 5mm, whilst the latter hosts just 4 struts, where SS is significantly increased as a result; an approach similar to that proposed in [98].This is repeated in Study B, where a variety of initial coating masses and polymer diffusion coefficients are considered (Table 4).Perturbations in the latter influence the percentage of drug released (RP), highlighted in Fig. 3.The polymer diffusion coefficient (D p ) is estimated such that we observe: (RP1) very slow release in which much less than a third of the drug has been released within 30 days; (RP2) a sustained release such that approximately 80% of drug is released across 30 days [99]; and (RP3) a very quick release (≈10 days).The initial drug masses are quoted per strut.The largest drug mass (DM3) was selected to portray a drug loading consistent with other models in the literature that are based upon commercial stents [75].Since this is at the higher end of drug loading on commercial stents, we proposed three arbitrary smaller masses (DM0-2) in an attempt to assess the impact of lower initial drug loadings on restenosis.This idea has emerged recently in various clinical trials involving either sirolimus [100,101] or other drug compounds such as paclitaxel [102,103], where DES drug loadings are comparable to the lower values considered in our analysis.However, when analysing perturbations in SS, drug masses are reported with respect to the entire stent (Supplementary material Table S3).
Study C then explores the influence of a variable binding site density

Table 3
Summary of the different cases considered.MM refers to Michaelis-Menten kinetics, which are described with either free or bound drug.If drug is considered, then in all cases, variations in the initial drug mass (DM) and drug release profile (RP) are simulated (Table 4).for an % over-expansion of 20% and each of the DM and RP values in Table 4.
Finally, in Study D, two different approaches to modelling the influence of drug on SMC growth are considered, each based upon Michaelis-Menten kinetics, but differing by whether free drug (Eq.(S.1)) or bound drug (Eq.(S.2)) dictates the effect on cell proliferation.

Analysis of results
To facilitate comparison between the different cases considered, a variety of different performance indicators are considered.Of most significance is % area restenosis (% AS).Assuming that the smallest radius within the stented region defines the severity of restenosis, we propose that the maximum level of restenosis present can be computed, assuming a cylindrical geometry, by comparing the radius pre and posthealing: where r l,g is the minimum radius of the lumen within the stented region following remodelling and r l,ref is the radius of the artery immediately after the stenting.The level of tissue growth is inherently linked to sSMC proliferation and ECM synthesis, which are coupled to multiple other processes involved in the healing response (Fig. 4).
In the absence of drug, tissue growth over a 300-day period is considered, similar to [38].However, when the performance of DES is analysed, the period of interest is reduced to 50 days which significantly reduces computational time, whilst still capturing the overwhelming majority of tissue growth.Focusing particularly on drug efficacy and safety, two time-dependent quantities are analysed across this period of interest.Firstly, drug content (DC, μg drug/g tissue), associated with safety, is given by: where V w denotes the volume of the arterial wall, ρ t is the density of wet arterial tissue, and M w is the molecular weight of drug (Table S2 Supplementary material).Simplifying the drug transport model to assume the same transport properties and continuity between layers means that this quantity, and others, do not need to be computed on a layer-by-layer basis.Moreover, in various studies that explore drug release kinetics, the saturation of specific binding sites has been associated with stent efficacy [54,59].In our baseline drug model, we assume cell proliferation is mediated by specific receptor saturation (RS s ), which is defined by: where subscript T denotes the different binding site considerations, either constant (T = c) or variable (T = v).

Study A: Impact of expansion and geometry in the absence of drug
Fig. 5 displays the results of Study A where the impact of overexpansion and stent geometry are assessed over a period of 300 days in the absence of drug delivery.Across all cases in Study A, approximately 90% of growth occurs within the first 50 days as shown for the baseline configuration in Fig. 5a, providing rationale for reducing simulation time within subsequent drug simulations (Study B-D).This is also in agreement with [38], where the majority of growth is present within the first 2 months.Moreover, despite significant changes in the model from [38], trends present in restenotic species behaviour are broadly similar.Table S6 (Supplementary material) presents the quantitative values of % area restenosis at 50 and 300 days for each case in Study A.
Our model is able to qualitatively capture trends present within the clinical and modelling literature.Firstly, within a given geometry, three different clinically relevant expansion profiles were considered.In agreement with physiological trends, as the severity of expansion increases, the level of restenosis observed increases (Fig. 5a) [27,52,88,104].This is simply a result of increased vessel stress which corresponds to higher levels of damage initially throughout the arterial wall.Similarly, increasing strut thickness results in an increase in % area restenosis, a trend observed clinically [105] and across various models of ISR (Fig. 6b) [25,27].Our model indicates near-negligible differences in restenosis across three (30/50/70%) pre-embedment profiles (Fig. 5c) as observed in [38].This result may be explained by the simplified stent deployment adopted in our model, in that malapposition or deep wall tears are not possible.Finally, a comparison between different stent designs within a fixed domain demonstrates that an increased number of struts (reduced inter-strut spacing) produces greater levels of restenosis.Similar trends were found in [98], where high (SS small ) and low stress (SS large ) stent model configurations within an idealised cylindrical geometry were compared to histological data.

Study B: Impact of stent geometry following drug release
Considering the performance of DES now, only two expansion profiles (20/30%) where substantial restenosis is present are considered.Our moderate drug release profile, which emulates release from a Xience TM (Abbott Vascular, Santa Clara, CA, USA) stent (RP2), elutes approximately 80% of drug within 30 days [99], and close to 100% across the full 50 days.Although other articles which focus on drug retention tend to look at shorter time frames of 28-30 days [54,58,59], we believe the extension to 50 days is valid to ensure that the majority of tissue remodelling has occurred, and the subsequent effect of drug on this process can be appropriately demonstrated.
Fig. 6 presents the results of Study B1 which explores the influence of varying RP and DM across different expansion profiles for the baseline model.Firstly, our model demonstrates that, for a given expansion level (20%), increasing DM within the stated range leads to lower levels of restenosis for each RP considered.However, the non-monotonic nature of the curves reveal an intricate relationship between DM and RP.Specifically, it is not the case that simply delivering drug more quickly and/or increasing the drug loading will result in reduced levels of restenosis.For example, within a given expansion profile (20%), when a larger drug mass is simulated (DM3), a slower release profile (RP1) presents more optimal than the quicker one (RP3).However, when the initial drug mass is lowered (DM1), the quicker release profile outperforms the slower one.For either DM1 or DM3, the moderate release profile (RP2) performs best, albeit indifferent from RP1 in the latter.Therefore, our results demonstrate that for a sufficiently high initial drug loading, a sustained release will present a more continuous saturation of receptors at the cost of elevated drug content levels over a prolonged period.However, for lower drug loadings, enough drug has to be distributed throughout the wall quickly within the active phase of remodelling such that sufficient saturation of receptors is apparent.Full details on the receptor saturation profiles and tissue drug content levels are illustrated in Figs.S2 and S3, respectively, within the Supplementary material.
Considering the impact of the initial drug loading, we note that as DM is increased, the influence on restenosis becomes progressively weaker, a consequence of the model formulation where increased receptor saturation correlates directly to a reduction in sSMC proliferation (Fig. 6b).In other words, increasing the drug mass further would have a negligible effect at the possible expense of elevated DC levels.At lower drug masses, differences between expansion profiles are more prevalent, with restenosis levels decreasing as the initial drug loading is increased.Moreover, for a quicker release profile (RP3), across the same drug mass range, a less optimal response is observed.
Figs. 7 and 8 consider a variety of different stent configurations in the presence of drug delivery for a single expansion ratio.Fig. 7a and b consider the effect of varying strut thickness.As the thickness is increased, we note the same trends across the three drug release profiles, where the moderate release presents as most optimal (Fig. 7a).As the drug mass is increased, the inhibition of tissue growth is more profound (Fig. 7b).However, similarly to the baseline case, for each ST considered, it appears that % area restenosis may be reaching a plateau, whereby increasing DM further will have little impact on restenosis.This is likely due to the inter-strut spacing condition (7 × T st ), where each configuration has a different distance between two adjacent struts.Another possibility may be that a constant binding site density is assumed, which may not fully capture the dynamic response present across the different configurations.
Fig. 7c and d chart how perturbations in the drug release profile and drug mass vary the level of restenosis observed when the initial strut embedment profile is altered.The curves follow similar trends as the baseline study, A1.However, within our model that includes drug delivery, deeper stent deployment results in lower levels of restenosis.At lower drug masses, the difference in restenosis between configurations is heightened.As expected, with increased embedment we note an improvement in stent performance, where increased concentrations of drug are present in the tissue at early times.
Next we consider how inter-strut spacing influences stent performance through different designs.To emphasise the differences between the results, analysis is conducted based on the 30% expansion profile.As is apparent in other configurations, the moderate drug release profile (RP2) inhibits restenosis most profoundly (Fig. 8a).Between the designs, a more significant inhibition of tissue growth is present in SS large where spacing is higher, assuming the total drug mass is constant between the two configurations (Fig. 8b).This suggests, similar to increased stent expansion or ST, that more drug mass is required to overcome greater levels of restenosis in designs that present higher levels of stress within the wall.
Across the various configurations considered it is clear the release profile and initial drug mass are significant factors when optimising the design of DES, where precise control of release kinetics is essential to Highlights the differences between cases at higher drug masses.In all studies, when varying RP the mass is fixed at DM1 and when varying DM, the release profile is fixed to RP2.In all cases, ND corresponds to the drug-free simulation.

Study C: Significance of binding site density
Assuming that receptor saturation influences stent performance explicitly through the inhibition of sSMCs, then the definition of the binding site density is critical.However, as suggested in [96,97], an upregulation of FKBP12, the receptor associated with sirolimus binding, is present at early times, decreasing as the artery heals and the neointima matures.Exploring this further, an amendment to the binding site density is considered (Eq.(2.26)) such that the saturable limit evolves in space and time as a function of the synthetic cell population.Table S4 (Supplementary material) summarises the different approaches and values considered for R c , the non-dimensional parameter used to detail receptor expression per cell.
In Fig. 9a-c we explore three different values for R c as discussed in Section 2.7.1.These figures plot the temporal behaviour of B s v across a range of strut drug masses averaged over the therapeutic domain, where the density is notably greatest behind the stent struts, where damage is more profound.
Trivially, as the R c value is increased we note an increase in receptor density.Adopting three different estimates for the computation of R c allows us to observe a variety of behaviours in the absence of reliable estimates.At lower values (R 1 c ), we note for all curves that the receptor density expression does not surpass the steady state value estimated in [54] (Fig. 9a).Thus, it is likely the value is an underestimate of true receptor expression present in the wall.Both R 2 c (Fig. 9b) and R 3 c (Fig. 9c) In all studies, when varying RP the mass is fixed at DM1 and when varying DM, the release profile is fixed to RP2.In all cases, ND corresponds to the drug-free simulation.
Refer to Table S3 (Supplementary material) for details on the initial drug mass.present a higher expression of receptors when compared directly to the constant assumption in the absence of drug (ND), decreasing as cells differentiate back to their contractile phenotype such that expression is low at later times, as suggested in [96,97].However, as depicted in Fig. 9a-c, the binding site density is drug-dependent, where at higher drug masses, the number of binding sites decreases in response to reduced sSMC proliferation.Only R 3 c presents a binding site density which peaks above B s c in all cases.Although receptor expression (Fig. 9a-c) between the three studies (C1-C3) is different, at higher drug masses (DM2 and DM3) there is a negligible difference in restenosis level observed as charted in Fig. 9d.In other words, for these initial drug masses and this moderate release profile (RP2), receptor saturation is sufficiently high such that proliferation and restenosis are low.At lower masses, we note higher levels of restenosis when the R c value is higher.This is a result of receptor saturation being less profound, a consequence of increased receptor density.
In Fig. 10, we compare how receptor saturation varies between the constant (B s c ) and variable (B s v for R 2 c ) binding site density approaches where the number of saturated receptors increases as the initial drug mass increases.
When the binding site density is presumed constant, the receptor saturation curve rises steeply initially, followed by a slower rate of increase until peak saturation is achieved, after which saturation levels gently decline.The higher the initial drug mass, the steeper the rate of increase, and the higher the value of peak saturation observed (Fig. 10a).However, when a variable approach (B s v ) is considered, a more peculiar receptor saturation plot is apparent (Fig. 10b).This is particularly evident at lower drug masses, where an initial spike is present, rapidly declining over the first few days before steadily increasing again.The peak is a result of the initial condition prescribed for the sSMC population to avoid numerical error, where the low number of receptors are quickly saturated as a result of the rapid binding kinetics (k s on ).The fall from this spike is a consequence of phenotype conversion, where cells de-differentiate from their contractile phenotype, increasing the receptor population more quickly than drug can be eluted for the considered release profile (RP2).The linear increase which follows is a consequence of the rate of change of bound drug and sSMC density being comparable.However, it is evident for the case DM1 that the curve eventually begins to plateau due to a depletion in bound drug and decrease in sSMC density.These behaviours are less apparent at higher values of DM because of the larger masses of drug entering the artery.Similar trends are present for the receptor saturation curves following perturbations in RP (Fig. 11a).However, with the binding site density now variable, as the value of R c increases, the quicker release profile (RP3) overtakes the moderate selection (RP2) as the optimal choice for inhibiting restenosis (Fig. 11b).This is noted only for R 3  c , implying that when the binding site density is increased, drug needs to be released more quickly to achieve adequate saturation of specific receptors during the active phase of remodelling for a fixed mass of DM1.This is emphasised in Fig. 11a, where the crossover between RP2 and RP3 occurs later as R c is increased.Similar to the assumption of a constant binding site density, as the initial drug mass is increased, the  quicker release profile is sub-optimal, and a sustained release presents superior.Taken together, these results suggest that the definition of the binding site density has an impact on the optimal controlled release strategy.
Regardless of whether the binding site density is assumed constant or variable, the correct trends are present when drug mass is varied.However, it is notable that, for the same initial drug mass and release profile, receptor saturation is higher and peaks more quickly for the variable binding site density approach across 50 days, compared with the constant assumption.In both interpretations of the binding site density, as the model suggests (Eq.(2.25)), receptor saturation dictates drug efficacy when coupled to the cell population directly.Results depict that with increased receptor saturation, a more profound inhibition of restenosis is present, illustrating the presence of a maximum drug dose, above which differences in restenosis are negligible.

Study D: Different drug models
Finally we investigate the influence of different models of drug effectiveness [95].This is discussed in depth within the Supplementary material (Section S3.4).In summary, although the free drug Michaelis-Menten model performed best in vitro, it is unlikely that it can be implemented within a model depicting in vivo without adequate amendments to account for drug retention, a feature that is present in the bound drug Michaelis-Menten model.Thus, the receptor saturation and Michaelis-Menten bound drug models present themselves as possible options to detail the effect of drug on sSMC proliferation, where trends present are broadly similar across the two.

Limitations
A number of assumptions have necessarily been made in devising the mathematical model presented here, as justified in the preceding text.The key limitations relate to (i) the complexity of the arterial wall (idealized 2D-axisymmetric with adventitia excluded) (ii) Stent deployment (prescribed displacement) (iii) Neglecting fluid dynamics in the lumen, and (iv) including only the key biological species, based upon current understanding in the literature.In the Supplementary material (Section S4) we provide an in-depth discussion of these points.
We note that, to our knowledge, there are no published experimental studies that would enable a complete validation of the model results presented here.We hope this work will inspire the acquisition of appropriate experimental data to facilitate model validation.Notwithstanding, the demonstration of the ability of our model to qualitatively capture trends present within the clinical and modelling literature provides confidence in our approach.

Conclusion
In this paper we have presented, to our knowledge, the first mathematical model of restenosis that is explicitly coupled with state-of-theart drug kinetics.The model is capable of simulating tissue growth following stent deployment and, importantly, reveals crucial insights into how the drug delivery protocol modulates restenosis.In particular, our results have revealed an intricate interplay between initial drug loading, drug release rate and restenosis, indicating that it is not sufficient to simply ramp-up the drug dose or prolong the time course of drug release.We have demonstrated that the level of over-expansion and stent design features influence the level of restenosis observed.The model results are in broad agreement with trends reported in the experimental and clinical literature.Taken together, our results highlight the necessity of coupling drug delivery with restenosis in the pursuit of optimal DES design.

Fig. 1 .
Fig. 1.(a) Longitudinal section view of the arterial wall prior to stent expansion.(b) Longitudinal section view following stent expansion.The inset highlights the polymer coating, intima, and denuded endothelium.Ω j defines the various domains, where subscript j = p, m, i denotes the polymer coating, media, and intima, respectively.Boundaries are denoted by Γ.The healthy and denuded endothelium are defined by Γ et,

Fig. 2 .
Fig. 2. Flow chart depicting the various stages of the model framework.Step 1 -Stent deployment triggers the restenosis model by coupling stress to damage.Step 2 -Restenosis ensues, where damage upregulates species (GFs and MMPs) which interact with other constituents to drive tissue growth.Step 3 -Remodelling, where an increase in constituents gives rise to tissue growth.

Fig. 4 .
Fig. 4. Deformed configuration depicts % overexpansion, highlighting difference in lumen radii pre and post-expansion.The healed geometry highlights volume of tissue growth, where severity in restenosis can be computed by the difference present between r l,ref and r l,g .

Fig. 3 .
Fig. 3. Stent release profile for various polymer diffusion coefficients for a sirolimus coated stent.

Fig. 6 .
Fig. 6.Significance of drug release profile (RP) and initial drug mass per strut (DM) on restenosis across two different expansion profiles (20/30%) for a period of 50 days.(a) Study B1 -RP: Solid lines refer to DM3 whilst dashed lines to DM1.(b) Study B1 -DM: Solid lines correspond to RP2, whilst dashed lines refer to RP3.Inset: Highlights the differences between cases where only drug is considered.In all cases, ND corresponds to the no drug simulation.

Fig. 7 .
Fig. 7. Significance of drug release profile (RP) and initial drug mass per strut (DM) on restenosis across different stent configuration regimes: variation in strut thickness (ST) or variation in strut embedment (SE) for a period of 50 days.All simulations are conducted within the 20% expansion profile.(a) Study B2 -RP.(b) Study B2 -DM.(c) Study B3 -RP.(d) Study B3 -DM.Inset:Highlights the differences between cases at higher drug masses.In all studies, when varying RP the mass is fixed at DM1 and when varying DM, the release profile is fixed to RP2.In all cases, ND corresponds to the drug-free simulation.

Fig. 8 .
Fig.8.Significance of drug release profile (RP) and initial drug mass per strut (DM) on restenosis when inter-strut spacing (SS) is varied.All simulations are conducted within the 30% expansion profile for a period of 50 days.(a) Study B4 -RP.Inset: Highlights the differences between release profiles.(b) Study B4 -DM.Inset: Highlights the differences between cases at higher drug masses.In all studies, when varying RP the mass is fixed at DM1 and when varying DM, the release profile is fixed to RP2.In all cases, ND corresponds to the drug-free simulation.Refer to TableS3(Supplementary material) for details on the initial drug mass.

Fig. 9 . 2 c( 1 . 72 ×
Fig. 9. Analysis of different receptor per cell (R c ) values.(a)-(c): How the receptor density (B s (t)) varies as a function of time, compared directly to the constant assumption (B s c -dashed black line) for each R c value: (a) = R 1 c (1.435 × 10 − 3 ), (b) = R 2 c (1.72 × 10 − 3 ), (c) = R 3 c (2.163 × 10 − 3 ).(d) Compares % area restenosis values between the three approaches across the 50 day time period.Inset: Focused on the different drug masses.All simulations were performed on 20% expansion profile and used RP2 release profile of the baseline stent configuration.In all cases, ND corresponds to the no drug simulation.

Fig. 10 .
Fig. 10.How receptor saturation varies between the two binding site density approaches: (a) B s c (solid lines) and (b) B s v = R 2 c c ssmc (dashed/dotted lines).Simulations are considered across various drug masses where RP2 is considered within the 20% expansion ratio of the baseline stent configuration.

Fig. 11 .
Fig. 11.Analysis of receptor saturation curves for different drug release profiles (RP) for each R c value when drug mass is fixed at DM1.(a) Receptor saturation curves for RP2 (orange) and RP3 (blue) across the different R c values.(b) % Area restenosis values across all R c values and drug release profiles (RP) across a simulation time of 50 days.

Table 4
Study B-D: Different parameter values for the effective diffusion coefficient of the polymer coating (D p ), which govern the different drug release profiles (RP) and the initial drug mass present in each strut (DM).