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

Abstract

An advection-diffusion solver was applied to assess how stent strut shape and position impact the development of a pro-thrombotic region within the stented human artery. Presented here is a suitably parameterised advection-diffusion equation with a source term that is spatially uniform within a certain sub-domain of interest to compute a “time concentration”. The latter will serve as a surrogate quantity for the “age” of fluid parcels, i.e., the time the fluid parcel has spent in the sub-domain. This is a particularly useful concept in the context of coronary artery haemodynamics, where “stasis of blood” (or residence time) is recognized as the most important factor in thrombotic initiation. The novel method presented in this work has a very straightforward and convenient single lattice Boltzmann simulation framework encapsulation. A residence time surrogate is computed, presented and correlated with a range of traditional haemodynamic metrics (wall shear stress, shear rate and re-circulation region shapes) and finally, the role of these data to quantify the risk of thrombus formation is assessed.

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 (EC), which line the arteries and play a vital role in minimizing inflammation and thrombosis. Consequently, scar tissue may form, leading to in-stent restenosis (ISR), with a reported incidence rate during the era of bare metal stents of ∼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 of $30\%$ [3,4] as well as a considerable burden of medication.

The dynamics of thrombus formation is multifactorial, 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 the following:

• 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.

The haemodynamic shear rate $\dot{\gamma }$ and endothelial shear stress (ESS) history play pivotal roles in activating advecting platelets [5]. This fact, along with the need to consider the relative residence time, makes it necessary to assess the history of fluid parcels for a quantitative understanding of ST. This work’s fundamental premise is that the introduction of stent struts into the vessel induces local haemodynamic perturbations which couple the effects 1.–3. above. This coupling, elaborated below, can be expressed using data from suitable single-framework fluid dynamics simulations. This hypothesis is tested 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 generalized as a rounded-square or as circular (see Figure 1). For a correctly deployed stent, elevated wall shear stress (WSS) is observed at the corners and along the top of the strut. In these regions, margination forces platelets down the shear gradient toward the boundary [6], while proteins such as the von Willebrand factor (vWF) form nets, thus increasing the number of platelet binding sites [7]. Meanwhile, distal to the strut, the flow separates into a recirculation region. Within such regions, the vessel wall experiences low WSS and platelets and vWFs 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]. This study investigates the extent to which this coupling is influenced by variation in strut geometry and the effectiveness of stent deployment (stent malapposition). In passing, low WSS is observed to result 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 hemodynamic behavior facilitates an understanding of how fluid transport links to ISR and ST. Most 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 into the artery wall, from 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 DESs [13]. Other mass transfer models have also been used with traditional CFD [14,15,16,17,18]. It is important to note that although our work primarily focuses on assessing thrombotic risk in coronary arteries, the base model presented here can also be applied to other flow studies where the residence time is a key metric. For example, in the study of cerebral aneurysms, the blood residence time is of significant interest [19], as disturbed flow patterns, such as stagnation and recirculation, can contribute to thrombus formation within the aneurysm sac. Of course, CFD lends itself not only to cardiovascular applications, but also to other health applications when combined with biomechanics, for example, hydrocephalus [20].

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 [21,22,23]. Exploiting these advantages, studies which are relevant to this work have applied LBM to blood flow in stenosed arteries [24,25,26,27], inferring relationships between flow dynamics, the stenosis and arterial atherogenesis. These works investigated flow around stents with particular reference to the evolution of ISR [28] to stent structure and positioning [29] and to stent selection [30]. Further afield, LBM has been applied to aneurysms, formed as artery wall remodelling [31,32]. Other studies have also investigated solving and applying the advection-diffusion equation to study biomedical flows using traditional approaches, such as in the work of Lynch et al. [33], who used a stabilising finite element approach to target a wide range of Péclet numbers. Traditional Galerkin finite element methods for advection-diffusion problems face significant challenges in high Péclet flows, including non-conservation, numerical diffusion and instability. Work has been done to overcome these problems; for instance, Hughes et al. [34] stabilised Galerkin methods by adding a residual term to ensure global conservation, resulting in better stability and accuracy. While stabilised Galerkin methods address conservation and stability through added residuals, LB solvers, such as the one presented here, inherently achieve local conservation and computational efficiency in advection-dominated flows governed by the Navier–Stokes continuity and convection-diffusion equations, and without requiring such modifications. We note in passing that while the latter formulation is appropriate for our RRT, stasis of blood application, it is not currently fitted to problems which require intimate adjustment of the Navier–Stokes equations, for example, to the Allen–Cahn–Navier–Stokes form [35].

The relative residence time (RRT), or stasis of blood, is most frequently calculated at the vessel wall, where its value is an accepted proxy for the probability of thrombotic initiation. Almost always, RRT is computed as a derived metric, obtained from the time-averaged instantaneous wall shear stress (AWSS) and oscillatory shear index (OSI), which are themselves obtained from the fluid’s velocity field. At the time of writing, the most widely used surrogate metric is assigned as follows:

$$RRT\propto {\left[\left(1-2\times OSI\right)\times AWSS\right]}^{-1}$$

(1)

Clearly, the above estimate of RRT is defined only at the wall. While this measure is useful in physiological flow calculations, at the boundary [36,37], to the authors’ knowledge, only the method advanced in this article can compute RRT at locations within the domain bulk, for example, in eddies. Furthermore, Equation (1) only provides a surrogate measure, not the residence time itself.

This work describes the implementation of a novel, single-framework meso-scale modelling approach, with the following:

• 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.

The model is then applied 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 the clotting algorithm and accounts for the computational domain, Section 3 reports the simulation results and Section 4 concludes this research, with clinical interpretation of the data and indications of future work.

2. Methodology

Two-dimensional (2D) simulations are certainly convenient to represent the scenarios in this research. 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 [21,22,23]). Implemented here is a single component and single relaxation LBM, widely designated the lattice Bhatnagar–Gross–Krook (LBGK) [38], which is one-way-coupled to an LBM advection-diffusion scheme with a source term. This advection-diffusion scheme effectively constitutes the clotting model, carried by a second lattice distribution. The flow and advection modes are a statistic of computations within a single framework. The following first introduces the flow LBM model, then the clotting algorithm, and finally, discloses 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).$$

(2)

where $c_s$ is the speed of sound for the model ($c_s^2={\displaystyle \frac{1}{3}}$) and $\tau$ is the relaxation constant for the scheme. These originate from the lattice gas cellular automation of [39]. At the time of writing, LBM variants have been derived in several forms, by multiple routes and authors (see e.g., [21,40]). The momentum and continuity equations are solved using the kinetic-scale LBGK evolution equation, which can be shown to be recovered through Chapman–Enskog analysis. For details on this analysis, readers are referred to [21,40].

$$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).$$

(3)

Above, i, $f_i$, $\mathbf{x}$, t, $\Delta t$, ${\mathbf{c}}_i$, $\rho$, $\mathbf{u}$, $\nu$ are the lattice link label, ith component of the PDF, discrete location, discrete time, time step, meso-scopic lattice velocity, 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),$$

(4)

where $t_i$ is the lattice link weight. This form of the equilibrium distribution function $f_i^{\left(0\right)}$ satisfies the conservation laws for mass and momentum and incorporates the influence of the local density $\rho$ and velocity $\mathbf{u}$ in a way that is consistent with the kinetic theory. For further justification and derivation of this form, readers are referred to [21,40]. The 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 .$$

(5)

The LBGK scheme outlined above is suitable for addressing the 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

Central to this work is 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 it is hypothesised that 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 this work considers only the relative residence time, the method will generalise straightforwardly to species advection-diffusion.

To assign a residence time surrogate, it is imagined that as the instant fluid parcels enter the region of interest, an internal source producing a time 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 [41] (see Appendix A.1). LBM-based advection-diffusion dynamics have also been variously derived [42,43,44]. Employed here is a mature variant developed by Shi and Guo [45], which contains a source term, S (see also Chang et al. [46]). 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_0,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,$$

(6)

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),$$

(7)

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{{\mathbf{c}}_i\mathbf{u}}{c_s^2}}\right].$$

(8)

Note that 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$.

To provide further clarity of the full computational algorithm, please see Figure 2 below, which shows the coupling between the population and species solver. Note that the evolution of the species distribution function, $g_i$, is one-way coupled to the evolution of the distribution function for the overall mass and momentum transport, $f_i$. One way to introduce two-way coupling would be to make the fluid kinematic viscosity (collision parameter $\tau$) depend upon the species concentration. In this work, however, the species concentration, g, is a proxy for the fluid residence time (stasis of blood).

Before targeting stented geometries, the model is validated for 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, the known analytical solution for the velocity profile is used to recover the domain residency time in the channel, and this result is compared with simulation data at a range of diffusion coefficient values. The results show that by tuning the relaxation parameter ${\omega }_d$, the model sufficiently simulated the desired $Pe$ regime, where for the remainder of the work, a value of ${\omega }_d=1.98$ is used [27].

2.3. Computational Domain

The computational domain is an idealised two-dimensional channel, with a fixed number of nodes in the ordinal directions. Specifically, the number of nodes along the x(y)-directions is denoted $X_L$ ($Y_L$) (see Figure 3). 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\mathsf{\mu }$m and length $50\mathsf{\mu }$m. Figure 1 shows the two strut profiles considered in this work, both having the same maximum width $100\mathsf{\mu }$m. The resolution of the domain is such that 1 lu = 1 $\mathsf{\mu }$m and the physical-lattice scaling used can be seen in

Table 1. Scales: Each column shows the physical scale range, simulated scales and the simulated scale in lattice units.

Linking Physical to LB Scales: 1 $\mathsf{\mu }$m = 1 lu
	Physical Size Range	Simulated Size	Lattice Units
Stent Height	60 $\mathsf{\mu }$m–160 $\mathsf{\mu }$m	100 $\mathsf{\mu }$m	100 lu
EC Height	0.1 $\mathsf{\mu }$m–10 $\mathsf{\mu }$m	1.5 $\mathsf{\mu }$m	1.5 lu
EC Length	50 $\mathsf{\mu }$m–70 $\mathsf{\mu }$m	50 $\mathsf{\mu }$m	50 lu
Artery Diameter	3000 $\mathsf{\mu }$m	3000 $\mathsf{\mu }$m	3000 lu

. For additional information on the boundary closure methods, please see Appendix A.2 and Appendix A.3.

3. Results and Discussion

Here, the model is applied to study the local haemodynamic effects of variation in (i) strut geometry and (ii) position. To do this, two common stent geometries are simulated, depicted in Figure 1, and four different stent strut positions are investigated, specifically, 50% embedded, 25% embedded, 0% embedded and −100% embedded (malapposed strut). Additionally, the model is applied to an isolated strut, and then a group of three. Throughout, this work reports $Re\in [400,600]$ representative of haemodynamics in human coronary arteries. To address (i) and (ii), the results are compared between simulations using metrics proposed to identify a thrombotic state or region. These metrics are the 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 its 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 [49]. 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 hypothesised that the build-up of atheroma correlated with lower WSS [50]. Over the last couple of decades, more evidence has been provided to support this assertion [51,52,53], with low WSS regions now being characterised as WSS < 0.4–1.0 Pa [51,54,55]. The build up of atheroma is associated with regions of low oscillatory WSS, which typically occur 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 the wake behind the strut [56]. When discussing the results, high WSS is identified as WSS > 7 Pa [56].

• Shear Stress

-

Shear rates in medium-to-large-sized arteries generally fall within the range 100 to 1000 s⁻¹ [57]. 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 5000{\mathrm{s}}^{-1}$, the vWF is activated in the diseased vessel, which sees a narrowing of the artery, and forms long strings that increase platelet binding [7,58]. Accordingly, shear rates above ${\dot{\gamma }}_{\mathrm{crit}}$ are taken 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, this work is 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 ${\mathrm{s}}^{-1}$). Once identified, this information is used 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, a single stent strut is simulated (both circular and rounded stent strut geometries, as shown in Figure 1) at a range of stent positions (100 $\mathsf{\mu }\mathrm{m}$ malapposed, 0% embedded, 25% embedded, 50% embedded), at Re $=409,469,528$. In each case, the relative residence time is considered, and the velocity and shear rate within the domain where the haemodynamics are of most interest. Figure 4 and Figure 5 provide an example of simulation results corresponding to the circular and rounded struts in flow characterised by Re = 528. To facilitate comparison,

Table 2. Circular stent strut: Max WSS found at the top centre of strut. Relative location is position relative to the origin at the bottom centre of the strut.

Reynolds Number	Embedded %	WSS	Shear Rate	Residence Time
		Max (Pa)	Max (${\mathbf{s}}^{-\mathbf{1}}$)	Max (lu)	Relative Location ($\mathsf{\mu }$m)	Re-Circulation Length ($\mathsf{\mu }$m)
Re = 409	−100	11.4	5406	328	(345, 72)	45
	0	9.91	4455	337	(49, 40)	131
	25	8.9	3894	229	(50, 30)	70
	50	7.41	3156	114	(−33, 44)	24
Re = 469	−100	12.4	5881	385	(354, 60)	117
	0	10.9	4921	414	(48, 36)	140
	25	9.77	4304	270	(48, 11)	74
	50	8.15	3478	120	(50, 1)	25
Re = 528	−100	13.3	6319	441	(366, 51)	165
	0	11.9	5357	486	(47, 32)	160
	25	10.6	4690	299	(48, 11)	77
	50	8.9	3781	120	(50, 1)	26

and

Table 3. Rounded square stent strut: Max WSS found at the post top-left rounded corner. Relative location is position relative to an origin taken at the bottom centre of the strut.

Reynolds Number	Embedded %	WSS	Shear Rate	Residence Time
		Max (Pa)	Max ( ${\mathbf{s}}^{-\mathbf{1}}$)	Max (lu)	Relative Location ($\mathsf{\mu }$m)	Re-Circulation Length ($\mathsf{\mu }$m)
Re = 409	−100	12.3	6245	362	(328, 79)	170
	0	10.4	5053	360	(51, 49)	154
	25	9.03	4300	267	(51, 33)	87
	50	7.35	3457	151	(51, 14)	41
Re = 469	−100	13.3	6753	428	(333, 67)	231
	0	11.5	5582	438	(51, 45)	161
	25	9.97	4750	316	(51, 27)	94
	50	8.09	3810	167	(51, 9)	43
Re = 528	−100	14.3	7218	495	(337, 57)	304
	0	12.4	6077	510	(51, 42)	169
	25	10.9	5174	354	(51, 24)	99
	50	8.83	4143	172	(51, 7)	44

summarise the key metrics for both the circular and rounded struts at Re $=409,469,528$.

For both stent geometries, it was observed 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 the 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 the thrombotic region location based on stent strut positioning.

The impact of the stent geometry on the local haemodynamics and how this relates to thrombus formation are considered. 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. It was further observed that the maximum WSS increases with Re number and decreases as the stent becomes more embedded, as expected. Interestingly, across all simulations of the circular and rounded stents struts, the maximum WSS values are above those of the high WSS level (>7 Pa [56]). Figure 4 and Figure 5 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 the 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 in vWF activation and elongation.

Finally, the relative residence time is examined across the simulations. The BRT field variable is best evaluated within the context of flow along with the streamlines and shear rate plots. It is noted 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, one can 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 the 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 was assessed within the steady flow regime and is defined, for present purposes, by the intersection of the relative residence time with the other metrics of wall shear stress, shear rate and distal reattachment distance, tentatively to identify risk areas which promote thrombus formation. It is observed, 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 the max WSS being above the critical limit stated earlier (7 Pa) and reported in [56]. 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 ${\mathrm{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 a lower magnitude of WSS, lower shear rates, smaller re-circulation regions and smaller relative residence times.

3.2. Multiple Stent Struts Simulations

The 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, on how effectively the stent has been implanted.

Here, steady flow over three adjacent equi-spaced struts of both circular and rounded geometries is simulated at a range of stent positions (100 $\mathsf{\mu }\mathrm{m}$ malapposed, 0% embedded, 25% embedded and 50% embedded) at $Re=(464,533)$ to assess the flow dynamics between the stent struts. In each simulation, BRT and the Stokes stream function overlaid on a shear rate field are considered, across a cropped section of the domain where the local haemodynamics are of most interest. Figure 6 and Figure 7 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 arise for the inter-strut re-circulation zones interactions. From Figure 6, for the 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 the 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 6 and Figure 7 demonstrate a similar trend to the isolated struts, with less-embedded struts showing higher relative residence time. Observed is a more complex BRT field for multiple struts, although some established conclusions are supported by these results as follows: (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. Conclusions

The core development in this work is the introduction and application of 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, it is postulated that the concentration of a diffusion limited time species scalar field serves 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.

Like all in silico models, our algorithm for the simulation of the advection-diffusion process suffers limitations when one attempts to apply it in a wider context to that reported here—thrombus risk. Specifically, when applied in any large Péclet number regime (as the diffusion coefficient approaches zero), instabilities arise, manifesting initially as short-wavelength, high-frequency oscillations. Further, ubiquitous numerical diffusion effects will be present due to the compromise forced by the computational resource, the diffusion coefficient value and, again, numerical stability. To address this, we have, here, proposed a way to control diffusion effects, by subtracting an open channel solution to “cancel” diffusion effects associated with the wall’s no-slip condition. Finally, a more stable (but expensive) diffusion algorithm might enhance the stability at higher Péclet number flows; but this development would not impact the framework reported herein for a novel utilisation of advection-diffusion dynamics to identify thrombosis risk.

The data that are computed to find an assessment of risk associated with stent deployment scenarios and strut profiles are the 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. This work has considered a finite, but clinically representative, set of geometries, from which consistent observations emerge.

The 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 the von Willebrand factor (vWF) to elongate and hence encourage platelet activation. Importantly, flow streamline distributions imply that fluid parcels containing activated vWFs would then advect these species, resulting in the accumulation of key proteins and platelets in low-shear-rate, long-residence-time re-circulation zones. The 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. Highlighted also is the importance of the flow interaction between adjacent struts, modulated by the 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 the stent strut shape, relative position and malapposition all couple strongly to the local hemodynamic metrics, including the 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 correlate the individual modalities and to establish the role of a number of data science tools which require only such data as are presented here.