Analysis of Blood Stasis for Stent Thrombosis using an Advection-Diffusion Lattice Boltzmann Scheme

1. Introduction

Cardiovascular disease is the world’s leading cause of death [1], with over three million people suffering from atherosclerotic cardiovascular disease in the UK. Atherosclerosis, characterized by the narrowing or blockage of arteries due to plaque build-up, can be treated using percutaneous interventions, such as balloon angioplasty and stent implantation. Inevitably, such procedures damage the monolayer of protective endothelial cells (ECs), which line the arteries and play a vital role in minimising inflammation and thrombosis. Consequently, scar tissue may form, leading to in-stent restenosis (ISR), with a reported incidence rate during the bare metal stent era of $\sim 17-41\%$[2]. Although stent thrombosis (ST) is increasingly rare, with studies reporting incidence rates ranging from approximately $0.7\%$ to $16\%$, this infrequent yet severe complication carries a mortality rate of $30\%$ [3,4] and a considerable burden of medication.

The dynamics of thrombus formation is multi-factorial but the central pillars of a quantitative, clinically relevant understanding are physical transport (species advection-diffusion) and hydrodynamics, because thrombus forms in the presence of:

• elevated concentrations of transported chemical species and activated platelets (which have been exposed to high fluid shear);
• increased fluid residence time;
• suppressed fluid and wall shear rates.

Haemodynamic shear rate, $\dot{\gamma }$, and endothelial shear stress (ESS) history play pivotal roles in activating advecting platelets [5]. This fact, together with the need to consider relative residence time, makes it necessary to assess the history of fluid parcels, for a quantitative understanding of ST. Our fundamental premise is that the introduction of stent struts into the vessel induces local haemodynamic perturbations which couple effects $1.-3.$ above. This coupling, elaborated below, can be expressed using data from suitable single framework fluid dynamics simulations. We test this hypothesis by measuring residence times and stress due to fluid mechanics in the presence of a set of clinically relevant flow perturbations, representative of both stent design and stent deployment scenarios.

Stent struts vary in size and shape and for this work are generalised as rounded-square or circular (see Figure 2). For a correctly deployed stent, elevated wall shear stress (WSS) is observed at corners and along the top of the strut. In these regions, margination forces platelets down the shear gradient, towards the boundary [6] while proteins like the von Willebrand factor (vWF) form nets, thereby increasing the number of platelet binding sites [7]. Meanwhile, distal to the strut, flow separates into a re-circulation region. Within such regions the vessel wall experiences low WSS and platelets and vWF accumulate. Clearly, the strut and the endothelial cells immediately distal to it define a pro-thrombotic domain in which stent-induced haemodynamic perturbations couple to trigger thrombus formation [8]. We investigate the extent to which this coupling is influenced by variation in strut geometry and the effectiveness of stent deployment (stent malapposition). In passing, we remark that low WSS results in vessel inflammation and smooth muscle cell proliferation, leading to the accumulation of smooth muscle cells [9], contributing to ISR in this same region.

The application of computational fluid dynamics (CFD) to local haemodynamic behaviour facilitates an understanding of how fluid transport links to ISR and ST. A majority of the prior art has used traditional, proprietary CFD software to investigate stent strut design and stent positioning, etc. [10,11,12]. The coupling of a mass transport equation has also been used by various authors to model drug transport, both into the lumen and the artery wall, from a drug eluting stents (DES). In some of the earliest work, Balakrishnan et al. implemented a CFD and mass transfer model to predict drug deposition for single and overlapping DES [13]. Other mass transfer models have also been used with traditional CFD [14,15,16,17,18].

The use of meso-scale modelling, such as the lattice Boltzmann method (LBM) has precipitated notable advances which rely on key attributes of the method, canonically (i) a facility with complex geometries, (ii) an intrinsically parallelisable algorithm and (perhaps most important) (iii) compliance with multi-physics and multi-scale modelling approaches [19,20,21]. Exploiting these advantages, studies which are relevant to this work have applied LBM to blood flow in stenosed arteries [22,23,24,25], inferring relationships between flow dynamics, the stenosis and arterial atherogenesis. These works investigated flow around stents with particular reference to the evolution of ISR [26], to stent structure and positioning [27] and to stent selection [28]. Further afield LBM has been applied to aneurysms, formed as artery wall remodelling [29,30]

Here, we will delineate and implement a novel, single framework meso-scale modelling approach, with:

• a “clotting algorithm", based upon a suitable asymptotic proxy for blood residence time;
• a simple and robust boundary closure algorithm for flow and advection-diffusion modalities, accurately dealing with geometrical complexity, such as endothelial cells and rounded stent struts.

We will apply this model to investigate and differentiate the patho-physiological haemodynamics of ST and ISR, using convenient (if approximate) metrics applied to rounded-square and circular stent struts, within the following workflow. Section 2 provides an overview of the D2Q9 LBM model, discusses our clotting algorithm and accounts for our computational domain, Section 3 reports our simulation results and Section 4 concludes our paper, with clinical interpretation of the data and indications of future work.

2. Methodology

Two-dimensional (2D) simulations are certainly convenient to represent our scenarios. Moreover, they are deemed physiologically adequate for the present application [18], in which flow is dominated by the transverse and stream-wise co-ordinates.

LBM is an unsteady, bottom-up Eulerian approach in which the probability density function (PDF) for the momentum of a fluid parcel, $f_i$, is evolved on a discrete lattice (see Section 2.1). LBM is conceived at the length-scale native to a minimal form of discrete kinetic theory and is designed to recover the continuum Navier-Stokes equations. This occurs in the thermodynamic limit, utilising a simple lattice evolution rule, or the kinetic equation, in which the momentum PDF streams between and interacts on lattice nodes (see [19,20,21]). Here, we implement a single component and single relaxation LBM, widely designated lattice Bhatnagar-Gross-Krook (LBGK) [31], which is one-way-coupled to an LBM advection-diffusion scheme with a source term. This advection-diffusion scheme effectively constitutes our clotting model, carried by a second lattice distribution. Flow and advection modes are a statistic of computations within a single framework. We first introduce the flow LBM model, then our clotting algorithm and finally disclose the straightforward coupling between the two.

2.1. Flow Solver Algorithm

The LBM-based descriptions of unsteady, weakly compressible hydrodynamics are as follows:

$${\displaystyle \frac{\partial }{\partial t}}\mathbf{u}+\mathbf{u}·\nabla \mathbf{u}=-{\displaystyle \frac{1}{\rho }}\nabla p+\nu {\nabla }^2\mathbf{u},{\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\mathbf{u}=0,\nu =c_s^2\left(\tau -{\displaystyle \frac{1}{2}}\right).$$

(1)

These originate from the lattice gas cellular automation of [32]. At the time of writing, LBM variants have been derived in several forms, by multiple routes and authors (see e.g., [19,33]). Solution of the momentum and continuity equations above is supported by evolution of the LBGK momentum PDF, $f_i$, using the following kinetic equation:

$$f_i(\mathbf{x}+\Delta t{\mathbf{c}}_i,t+\Delta t)=f_i(\mathbf{x},t)-{\displaystyle \frac{\Delta t}{\tau }}\left(f_i(\mathbf{x},t)-f_i^{\left(0\right)}(\rho ,\mathbf{u})\right).$$

(2)

Above i, $f_i$, $\mathbf{x}$, t, $\Delta t$, ${\mathbf{c}}_i$, $\tau$, $\rho$, $\mathbf{u}$, $\nu$ are the lattice link label, $i^{\mathrm{th}}$ component of the PDF, discrete location, discrete time, time step, meso-scopic lattice velocity, relaxation constant, macroscopic density, macroscopic velocity (as defined below) and kinematic viscosity. The equilibrium distribution function, $f_i^{\left(0\right)}(\rho ,\mathbf{u})$, is a second order polynomial approximation to a uniformly propagating Maxwell-Boltzmann distribution:

$$f_i^{\left(0\right)}(\rho ,\mathbf{u})=t_i\rho \left(1+{\displaystyle \frac{{\mathbf{c}}_i\mathbf{u}}{c_s^2}}+{\displaystyle \frac{{\mathbf{c}}_i^2{\mathbf{u}}^2}{2c_s^4}}-{\displaystyle \frac{{\mathbf{u}}^2}{2c_s^2}}\right),$$

(3)

where $t_i$ and $c_s$ are a lattice link weight and the speed of sound for the model $c_s^2={\displaystyle \frac{1}{3}}$, respectively. Macroscopic fields $\rho$, $\mathbf{u}$ and P are calculated through summations of the distribution function in the following way:

$$\sum _i{f}_i=\rho ,\sum _i{f}_i{\mathbf{c}}_i=\rho \mathbf{u},P={\displaystyle \frac{1}{3}}\rho .$$

(4)

The LBGK scheme outlined above is suitable for addressing Reynolds numbers characteristic of coronary artery flow:

$$Re={\displaystyle \frac{\left|\mathbf{u}\right|L}{\nu }}\approx 600.$$

Above, L characterises the spatial scale of the domain. The clotting algorithm component (for residence time) discussed in the next section, is one-way coupled to this scheme, i.e., the development of the residence time solution is influenced by the fluid dynamics, but not the other way around. The simple coupling is achieved as a parameterisation of the passing of the macroscopic velocity $\mathbf{u}$ to the clotting algorithm, which is then used in the advection term.

2.2. Domain Residence Time Calculation

Here we address the identification of fluid parcels within the domain having an elevated residence time, relative to a threshold to be determined. The residence time algorithm is intended, eventually, to inform a reduced order model of clotting, by identifying parcels which are sequestered within a finite region in which we hypothesise, their cargo of metabolites will not be refreshed. The time for which a fluid parcel has been in the simulation domain -its age- is defined as its residence time. When regions of high residence time exist adjacent to a boundary and, further, contain elevated concentrations of vWF and activated platelets, one can infer an increased risk of thrombus formation within that region. While we consider only relative residence time, our method will generalise straightforwardly to species advection-diffusion.

To assign a residence time surrogate, we imagine that at the instant fluid parcels enter the region of interest, an internal source producing a times species at a uniform rate is activated. The domain residence time scalar field value is then directly proportional to the time species concentration. Of course, time does not diffuse, so that computing a time species concentration, denoted $\varphi$ (from the solution of the advection-diffusion equation, with a constant source) requires us to parameterise within the regime of a very large Peclet number [34]. See Appendix A.1). LBM based advection-diffusion dynamics have also been variously derived [35,36,37]. Here, we employ a mature variant developed by Shi and Guo [38] which contains a source term, S (see also Chang et al., [39]). The advection-diffusion equation, with a homogeneous diffusion coefficient D, below:

$${\displaystyle \frac{\partial \varphi }{\partial t}}=D{\nabla }^2\varphi -\nabla ·\varphi \mathbf{u}+S,D=c_s^2\left({\displaystyle \frac{1}{{\omega }_d}}-{\displaystyle \frac{1}{2}}\right),$$

is supported by the evolution of a second PDF, $g_i$, as follows:

$$g_i(\mathbf{x}+\Delta t{\mathbf{c}}_i,t+\Delta t)=g_i(\mathbf{x},t)-{\omega }_d\left(g_i(\mathbf{x},t)-g_i^{\left(0\right)}(\varphi ,\mathbf{u})\right)+\Delta S_i,$$

(5)

where the equilibrium distribution function $g_i^{\left(0\right)}(\varphi ,\mathbf{u})$ is:

$$g_i^{\left(0\right)}(\varphi ,\mathbf{u})=t_i\varphi \left(1+{\displaystyle \frac{{\mathbf{c}}_i\mathbf{u}}{c_s^2}}+{\displaystyle \frac{{\mathbf{c}}_i^2{\mathbf{u}}^2}{2c_s^4}}-{\displaystyle \frac{{\mathbf{u}}^2}{2c_s^2}}\right),$$

(6)

and its kinetic evolution equation source term $S_i$ is

$$S_i=t_i{S}_0\left[1+\left(1-{\displaystyle \frac{{\omega }_d}{2}}\right){\displaystyle \frac{\varphi {\mathbf{c}}_i\mathbf{u}}{c_s^2}}\right].$$

(7)

Note, the diffusion coefficient is linked to the second distribution relaxation parameter ${\omega }_d$, that the velocity of the fluid $\mathbf{u}$ is a parameter of $g_i^{\left(0\right)}(\varphi ,\mathbf{u})$ and that the kinetic equation of the second distribution is, essentially, an LBGK scheme close to that of the momentum PDF. The source term is used to increase the residence time surrogate $\varphi$ at the specified, uniform rate, $S_0$.

Before targeting stented geometries, we validated the application of the advection-diffusion scheme in an unstented geometry (straight channel) in the large $Pe$ regime, i.e., where advection dominates and diffusion effects are approximately zero. To do this, we used the known analytical solution for the velocity profile to recover domain residency time in the channel, and compared this result with simulation data at a range of diffusion coefficient values. Results showed that by tuning the relaxation parameter ${\omega }_d$, we can sufficiently simulate the desired $Pe$ regime, where for the remainder of the work we use an ${\omega }_d=1.98$ [25].

2.3. Computational Domain

The computational domain is an idealised rectangular channel, consisting of $X_{LENGTH}$×$Y_{LENGTH}$ nodes, see Figure 1. The south boundary of the domain is modelled to represent a physical endothelium, consisting of endothelial cells which are modelled with a height of $1.5\mu m$ and length $50\mu m$. Figure 2 shows the two strut profiles considered in this work, both having the same maximum width $100\mu m$. The resolution of the domain is such that 1 lu = 1 $\mu$m and the physical-lattice scaling used can be seen in Table 1. For additional information on our boundary closure methods please see Appendices Appendix A.2 and Appendix A.3.

Figure 1. Schematic to illustrate our simulation domain$X_{\mathrm{LENGTH}}=1450$ lu and $Y_{\mathrm{LENGTH}}=1500$ lu. The boundary conditions are: North symmetry condition; East / West algorithm of Zhang and Kwok [40]; South boundary (Endothelial cells and Stent Strut) non-slip condition enforced by the Bouzidi boundary Condition [41]. See Table 1 for physical units of the simulation.

Figure 1. Schematic to illustrate our simulation domain$X_{\mathrm{LENGTH}}=1450$ lu and $Y_{\mathrm{LENGTH}}=1500$ lu. The boundary conditions are: North symmetry condition; East / West algorithm of Zhang and Kwok [40]; South boundary (Endothelial cells and Stent Strut) non-slip condition enforced by the Bouzidi boundary Condition [41]. See Table 1 for physical units of the simulation.

Figure 2. Stents Struts: Schematic of the two stent struts being simulated in this paper: - Stent 1: Circular stent strut - diameter 100$\mu$m. - Stent 2: Rounded stent strut - width and length 100$\mu$m.

Figure 2. Stents Struts: Schematic of the two stent struts being simulated in this paper: - Stent 1: Circular stent strut - diameter 100$\mu$m. - Stent 2: Rounded stent strut - width and length 100$\mu$m.

3. Results and Discussions

We will apply our model to study the local haemodynamic effects of variation in (i) strut geometry and (ii) position. To do this, we will simulate two common stent geometries, depicted in Figure 2 and investigate four different stent strut positions, specifically 50% embedded, 25% embedded, 0% embedded and -100% embedded (malapposed strut). Additionally, we will apply our model to an isolated strut, and then a group of three. Throughout, we report $Re\in [400,600]$ representative of haemodynamics in human coronary arteries. To address (i) and (ii) we compare results between simulations using metrics proposed to identify a thrombotic state or region. These metrics are relative blood residence time (BRT), wall shear stress (WSS), flow shear rate, re-circulation region size and location. It is useful to bound these metrics, to allow us to evaluate how local haemodynamics relates to a thrombotic state, and their relation to ISR and ST.

Blood Residence Time

- The residency time is a measure of the amount of time a fluid parcel has spent within the domain, and can be used to locate atheroprone regions within the domain [42]. In relation to stent haemorheology, regions of elevated residency time may accumulate platelets and adhesion proteins, promoting thrombus formation.

Wall Shear Stress

- In the 1970s it was hypothesized that the build-up of atheroma correlates with lower WSS [43]. Over the last couple of decades, more evidence has been provided to support this assertion [44,45,46] with low WSS regions now being characterised as $WSS<0.4-1.0$ Pa [44,47,48]. The build up of atheroma is associated with regions of low oscillatory WSS, which typically occurs in regions of flow separation and re-circulation. Higher WSS also plays a role in atheroma formation. The region of higher WSS found on the adlumninal surface of the struts results in platelet geometry change and activation. These platelets can then become trapped in the lower wall shear region, or wake behind the strut [49]. When discussing our results, we identify high WSS as $WSS>7$ Pa [49].

Shear Stress

- Shear rates in medium-to-large-sized arteries generally fall within the range 100 to 1000 $s^{-1}$[50]. However, stenosed regions experience larger shear rates. Indeed, it has been shown that at and above super-critical shear rates ${\dot{\gamma }}_{\mathrm{crit}}\approx 5000s^{-1}$, the vWF is activated in the diseased vessel which see’s a narrowing of the artery, and forms long strings that increase platelet binding [51,52]. Accordingly, we take shear rates above ${\dot{\gamma }}_{\mathrm{crit}}$ to indicate an increased risk of atheroma formation, caused by the effective local narrowing of the artery due to the presence of the strut.

To summarise, we are particularly concerned with regions which will link to atherosclerotic plaque formation, i.e., areas of low WSS ($<1$ Pa), areas of high WSS ($>7$ Pa), and regions of high shear rates ($>5000$$s^{-1}$). Once identified, we will use information on flow stream paths, derived from the Stokes stream function, in tandem with the relative residence time information, to identify elevated risk of atherosclerotic plaque formation, linking to biological responses.

3.1. Single Stent Strut Simulations

Here, we simulate a single stent strut (both circular and rounded stent strut geometries as shown in Figure 2) at a range of stent positions (100$\mu m$ malapposed, 0% embedded, 25% embedded, 50% embedded), at Re $=409,469,528$. In each case we consider relative residence time, velocity and shear rate within the domain- where haemodynamics are most of interest. Figure 3 and Figure 4 provide an example of simulation results corresponding to the circular and rounded struts in flow characterised by Re=528. To facilitate comparison, Table 2 and Table 3 summarise key metrics for both the circular and rounded struts at Re $=409,469,528$.

For both stent geometries, we observe that the more embedded the stent strut, the smaller the flow disturbance caused. This observation can be quantified by the re-circulation length given in Table 2 and Table 3. Comparing results for struts which are 50%, 25% and 0% embedded, for both the circular and rounded struts, the re-circulation size decreases as the percentage embedding increases. This suggests that for less embedded stents, there will be a larger distal region of low WSS. However the influence of stent geometry manifests when comparing circular and square stent struts when malapposed, with circular stent struts producing smaller re-circulation regions, due to their streamlined profile maintaining flow attachment for longer. In addition, for both stent geometries, malapposed struts cause a re-circulation region which is detached from the stent strut, in comparison to 0%, 25% and 50% embedded simulations, where the re-circulation region is attached to the strut. This observation indicates a potential change in thrombotic region location based on stent strut positioning.

Consider the impact of stent geometry on the local haemodynamics, and how this relates to thrombus formation. Comparing WSS over the stent struts it is evident that the rounded square strut has a larger maximum WSS than the circular strut, for equivalent Re. We further observe the maximum WSS increases with Re number and decreases as the stent becomes more embedded, as expected. Interestingly, across all simulations of circular and rounded stents struts, maximum WSS values are above those of the high WSS level ($>7$ Pa [49]). We will return to this point. Figure 3 and Figure 4 plot the shear rates in the region local to the stent. These plots clearly show, for all simulations, that the areas of low shear rates correspond to areas where re-circulation regions are present, and that areas of high shear rates appear at the top centre of the strut and the post top left rounded corner of the struts, for the circular and rounded square struts, respectively. For the circular strut, the malapposed positions resulted in shear rates above the ${\dot{\gamma }}_{\mathrm{crit}}$, whereas for the rounded square strut, the malapposed and 0% embedded positions across all simulations have shear rates above the ${\dot{\gamma }}_{\mathrm{crit}}$. This suggests that both stent position and geometry play a critical role vWF activation and elongation.

Finally we examine the relative residence time across the simulations. The BRT field variable is best evaluated within the context of flow along with streamlines and shear rate plots. We note that the regions with large relative BRTs are also regions where flow re-circulation occurs, at low velocities and low shear rates. In addition, by following the streamlines, we imagine how platelets and vWF proteins convect to these high residence time regions, generating a pro-thrombotic region and state. Furthermore, Table 2 and Table 3 show larger residence times when the stent is either malapposed or 0% embedded, when compared with the 25% and 50% embedded, which accords with a hydrodynamics of larger re-circulation regions, when the stent is either malapposed or 0% embedded.

Concluding the findings on single strut scenarios, the pro-thrombotic region has been assessed within the steady flow regime and is defined, for present purposes, by the intersection of relative residence time with other metrics of wall shear stress, shear rate and distal reattachment distance, tentatively to identify risk areas which promote thrombus formation. We observe, for both circular and rounded stent struts, that a re-circulation region forms behind the strut which is generally more prominent with less embedding of the strut. All struts experienced high WSS with max WSS being above the critical limit stated earlier (7 Pa) and reported in [49]. It was found that the square stent strut produced higher shear rates, with both the malapposed and 0% embedded strut producing shears above the stated critical value (5000 $s^{-1}$), and the malapposed circular struts also showing shear rates above the critical value. Following the streamlines from these areas of high WSS and high shear rates (where the platelets are hypothesised to be activated), the platelets will aggregate in the circulation region (region of low WSS) which will increase the risk of thrombus formation. It was also observed that the more streamlined circular strut reduced the risk of a pro-thrombotic region, by having: lower magnitude of WSS, lower shear rates, smaller re-circulation regions and smaller relative residence times.

3.2. Multiple Stent Struts Simulations

Flow dynamics in the stent region can be more complicated because of interactions between adjacent struts. The spacing between struts varies, depending upon the brand and the axial co-ordinate (i.e the position where the stent cross-section is examined) and of course how effectively the stent has been implanted.

Here we simulate steady flow over three adjacent equi-spaced struts of both circular and rounded geometries at a range of stent positions (100$\mu m$ malapposed, 0% embedded, 25% embedded and 50% embedded) at $Re=(464,533)$, to assess the flow dynamics between stent struts. In each simulation, we consider BRT and the Stokes-Stream function over-laid on a shear rate field, across a cropped section of the domain where local haemodynamics are of most interest. Figure 5 and Figure 6 provide an example of simulation results corresponding to the circular and rounded struts in flow characterised by Re=533.

A clear difference between the multiple stent strut simulations and the isolated strut simulations are inter-strut re-circulation zones interactions. From Figure 5, for circular stent struts, at the Re numbers tested, the re-circulation behind the first strut does not reach the second strut, when embedded by at least 25%. For the rounded square stent struts, the re-circulation behind the first strut does not reach the second stent strut when the strut is embedded by more than 25%. This implies that for rounded square struts there are larger regions which experience both low shear rates and, presumably, oscillatory flow. Clearly, the spacing between stent struts would also impact on whether the detached flow of one strut re-attaches before the next strut.

The residence time plots in Figure 5 and Figure 6 demonstrate a similar trend to the isolated struts, with less embedded struts showing higher relative residence time. We see a more complex BRT field for multiple struts, although some established conclusions are supported by these results (i) areas of larger BRT align with regions of low shear stress; (ii) circular strut simulations show lower BRT when compared to rounded square simulations; (iii) areas of larger BRT are found proximal to the stent struts, again, intersecting with low shear rate regions.

4. Conclusion

This paper introduces and applies a coupled LBM-advection-diffusion and Newtonian flow solver, with a straightforward and amenable boundary closure, all encapsulated within a single framework, to investigate artery scale stent haemodynamics, with particular emphasis on local, pro-thrombotic regions in the immediate vicinity of the struts. Specifically, we postulate the concentration of a diffusion limited time species scalar field as a surrogate, at high Peclet number, for the relative domain residence time of blood. This quantity is computed alongside a set of accepted mechanical metrics (see below) that might easily provide datasets from which to evaluate the potential risk of atheroma formation.

The data we compute to find an assessment of risk associated with stent deployment scenarios and strut profiles are shear rate, wall shear stress, relative residence time, and the characteristic length of the strut-distal re- circulation region. All are determined with appropriate resolution and address realistic physiological ranges, relative scales and geometries, within the steady regime. Simulations were performed both on single struts and triplets of struts, to determine how strut shape and positioning interact with the metrics identified above. We have considered a finite, but we believe clinically representative, set of geometries, from which consistent observations emerge.

Our general observations indicate that the less embedded a strut, the larger the risk of a region prone to thrombosis, based on identifying regions with large relative residence time, low shear rates and large re-circulation zones. High WSS over the top of the struts will promote von Willebrand factor (vWF) to elongate and hence encourage platelet activation. Importantly, flow streamline distributions imply that fluid parcels containing activated vWF would then advect these species, resulting in the accumulation of key proteins and platelets in low shear rate, long residence time re-circulation zones. Simulation results further show that the rounded square stent struts created regions with larger WSS, longer re-circulation characteristic length and higher relative residence time compared with the circular stent struts, which quantifies the significance of strut streamlining. Multiple stent strut simulations evidenced a similar trend, with an increased residence time associated with the rounded square stent strut. We also expose the importance of the flow interaction between adjacent struts, modulated by Reynolds number, stent position and stent spacing. All these factors impact upon whether the flow re-attaches before reaching the next stent strut, which appears to be of high importance, as it could lead to large regions of low and oscillatory shear stress.

Broadly, it is clear that stent strut shape, relative position and malapposition all couple strongly to local hemodynamic metrics, including residence time. From the data presented here, no single metric seems to be insignificant, highlighting the need for convenient single framework methodologies. This study confirms that this approach can report all the relevant fields efficiently and accurately. Further work is required to objectively to correlate the individual modalities- the role of a number of data science tools which require only such data as we present here.

Appendix A.1. Control of Diffusion

Eliminating diffusion entirely is not possible in LBM or any other numerical scheme. Simulations are therefore parameterized to emphasize advection over diffusion, by simulating at large Peclet numbers:

$$P{e}_L={\displaystyle \frac{L\left|\mathbf{u}\right|}{D}}.$$

We do not aim to facilitate precise solutions of advection-diffusion dynamics for $P{e}_L\to \infty$. Rather, to obtain a measure of relative residence time which supports quantitative comparison of clinically relevant scenarios in the form of values of $\varphi$ which stratify different patho-physiologies and characterise domain locations. This comparison might be facilitated by a normalisation, relative to the corresponding physiological case. Even so, a natural step is to minimise diffusion by parameterising as follows:

$${\omega }_d\to 2,D\to 0,P{e}_L\to \infty .$$

However, in LBGK simulation, the relaxation parameter ${\omega }_d$ cannot approach too closely its upper bound without exciting instability [19]. We defer further consideration of lattice $P{e}_L$ and $Re$ parameterisation using the lattice $\left|\mathbf{u}\right|$ and L. For results reported, we set ${\omega }_d=$1.98.

Appendix A.2. Treatment of Irregular Boundaries

One benefit of the LBM is its ability to represent complex geometries using straight-forward boundary closure rules. To model the correct boundary curvature characteristic of strut and endothelial cell geometries, an interpolated bounce-back scheme is implemented, following Bouzidi et al. [41]. This scheme determines the post-collision and post-propagation $\left(\psi \psi \right)$ distribution function $f_{-i}^{\psi \psi }(x_B,t+\Delta t)$ by interpolating between the post-collision pre-propagation $\left(\psi \right)$ distribution functions, at the fluid site immediately adjacent to the boundary $f_{-i}^{\psi }(x_B,t)$ and at the penultimate lattice site $f_i^{\psi \psi }(x_A,t)$. The formulation of $f_{-i}^{\psi }(x_B,t+\Delta t)$ depends upon a metric $q={\displaystyle \frac{\left|BW\right|}{\left|BC\right|}}$ which is the normalised distance to the boundary wall:

$$\begin{array}{cc}\hfill f_{-i}^{\psi \psi }(x_B,t+\Delta t)& =\left\{\begin{array}{cc}2q{f}_i^{\psi }(x_B,t)-(1-2q)f_i^{\psi }(x_A,t)\hfill & \mathrm{if}q<{\displaystyle \frac{1}{2}},\hfill \\ {\displaystyle \frac{1}{2q}}{f}_i^{\psi }(x_B,t)+{\displaystyle \frac{2q-1}{2q}}{f}_i^{\psi }(x_B,t)\hfill & \mathrm{if}q\ge {\displaystyle \frac{1}{2}},\hfill \end{array}\right.\hfill \end{array}$$

(A1)

To further exemplify this result and declare our notation we provide the one-dimensional schematic of Figure A1, below.

Figure A1. Bouzidi Boundary Conditions: Schematic depicting the Bouzidi boundary condition interpolation rules based on the two cases $q<0.5$ & $q\ge 0.5$, in one dimension. Note, the labelling convention used in this figure relates to that used in Equation (A1).

Figure A1. Bouzidi Boundary Conditions: Schematic depicting the Bouzidi boundary condition interpolation rules based on the two cases $q<0.5$ & $q\ge 0.5$, in one dimension. Note, the labelling convention used in this figure relates to that used in Equation (A1).

For closure of the advection-diffusion scheme, we also apply the interpolated bounce-back scheme on the $g_i$ distribution function. The remaining boundary conditions implemented for both the lattice Boltzmann and advection-diffusion solver are discussed in the next section.

Appendix A.3. Flow Boundary Conditions

Table 1 outlines: (i) the range of physical scales; (ii) the scales simulated; (iii) the related lattice units (lu). To ensure appropriate resolution, we use a conversion of 1 $\mu m$ = 1 lu. Figure 1 is a schematic of the simulation domain (top) and then an enlargement of the local flow dynamics around a strut. Figure 2 declares the two strut profiles considered in this work, both having the same maximum width $100\mu m$. Further, the ECs have a height of $1.5\mu m$ and length $50\mu m$.

A mix of boundary conditions have been implemented for the species and concentration populations. In Figure 1, we utilise the symmetry of an idealised artery channel for computational efficiency, a symmetry condition being implemented at the north boundary, for both the species and concentration populations. At the east and west boundaries for the concentration we implement periodic boundary conditions. For the momentum distribution $f_i$, we implement Zhang-Kwok pressure boundary conditions [40], which impose a pressure step across an otherwise periodic velocity boundary. The pressure in LBM is proportional to the density as $P=c_s^2\rho$. Equation (A2) describes how the pressure gradient is calculated by multiplying exiting distribution functions by constant $\beta$. The density (pressure) difference between the inlet and outlet will drive fluid in the streaming direction, by maintaining a higher (lower) density at the distal (proximal) arterial boundary.

$$f_{\mathrm{inlet}}^{\mathrm{in}}=f_{\mathrm{outlet}}^{\mathrm{out}}{\displaystyle \frac{{\rho }_0+\frac{\beta }{c_s^2}}{{\overline{\rho }}_{\mathrm{outlet}}}},f_{\mathrm{outlet}}^{\mathrm{in}}=f_{\mathrm{inlet}}^{\mathrm{out}}{\displaystyle \frac{{\rho }_0-\frac{\beta L}{c_s^2}}{{\overline{\rho }}_{\mathrm{inlet}}}}.$$

(A2)

We assign $\beta$ to ensure our simulations fall within the physiological Re number range for artery flow. See Section 3. At the south boundary, to represent endothelial cells and struts, we implement the no-slip interpolative bounce back boundary condition of Appendix A.2 for both species and concentration.