Investigation of hemodynamic markers for stenosis development

Correspondence James M. Buick, School of Mechanical and Design Engineering, University of Portsmouth, Anglesea Building, Portsmouth PO1 3DJ, UK. Email: james.buick@port.ac.uk Abstract Different hemodynamic markers for stenosis development are considered, assessed and compared. A recently proposed numerical approach is employed, where stenosis development is modeled based on the local hemodynamic conditions at the artery wall, determined using the lattice Boltzmann method. A range of hemodynamic markers, commonly associated in the literature with the progression of atherosclerosis, were considered. It was observed that using the oscillatory shear index, which is related to the wall shear stress oscillating from its mean direction, did not produce a realistic stenosis development. The other markers produced stenosis growth comparable with observations from the literature.


INTRODUCTION
Cardiovascular diseases, of which atherosclerosis is a major contributor, account for over 35% of all deaths in the western world. 1 It has been suggested that blood flow in the artery is a localizing factor in atherosclerosis. 2 Many studies have shown an association between regions of low wall shear stress (WSS) in an artery and the position when stenosis develops. [3][4][5][6][7][8][9][10][11][12][13] It has also been suggested that oscillations in the shear stress, measured through the oscillatory shear index (OSI) and/or the related relative residency time (RRT), are additionally important. 4,5,7,10,12,[14][15][16] From Gnasso et al.'s study, 3 Ku et al., 4 described increased residence time along the outer wall of the sinus which is caused by oscillations of fluid velocity about a mean value close to zero. This occurrence delays the convection of fluid and traps fluid elements near the outer wall for several cycles despite the absence of a clear region of stasis or of an area of permanent boundary layer separation. It has also been suggested 17 that OSI and RRT are the key hemodynamic parameters, rather than the time-averaged (TA) WSS. Other authors have suggested low near-wall velocity 7,8 and also high near-wall viscosity (due to the non-Newtonian nature of blood). 8,18 The lattice Boltzmann method (LBM) 19 is a numerical technique for performing fluid flow simulations which is based on a kinetic approach, rather than the more traditional computational fluid dynamics (CFD) approach of solving the Navier-stokes equation. Application of the LBM to stenosis growth related problems has been considered by a number of authors. Tamagawa and Matsuo 20 presented a model for blood clotting using a tracer particle which adhered to the wall following a rule based on it being in a region of low shear rate and also close to the wall. Tamagawa et al. 21 extended this work. This work was extended to also include surface tension and particle concentration in the adhesion rule. An alternative approach 22,23 also considered tracer particles with adhesion occurring based on the particles achieving a required residence time in a near-wall location. Moiseyev and Bar-Yoseph 24 considered coagulation of the blood based on an integrated accumulation of the stress and a sufficient residential time. Karimpour and Javdan 25 considered an LBM model where stenosis development was modeled based on the OSI at the artery wall. Regions of high OSI were deemed to be prone to atherosclerosis and layers of deposits were progressively put down in such regions. This resulted in a layer of deposit at each step. Stamou and Buick 26 simulated the early stages of stenosis development using a model where stenosis development was focussed in regions of low time-averaged velocity. In contrast to the model of Karimpour and Javdan,25 this model enabled the stenosis to develop in a local manner, rather than in layers, and provides a more localized and detailed description of the stenosis as it develops, albeit at a computational cost. Yan et al. 27 proposed an LBM model where cell adhesion was simulated using a Monte Carlo model. The results showed that the model was appropriate for simulating wall adhesion; however, it was focussed on an idealized micro-vessel. A further LBM approach 28 used a regular channel with similar dimensions to the carotid artery. A thrombosis process, involving platelets transported through the blood, was then introduced to provide a growth model at the artery wall.
In addition to applications specific to stenosis growth, the LBM has also been applied to a range of applications relating to the carotid artery; generally where simulation is performed either for a healthy artery, or for one with a prescribed level of stenosis. Pontrelli et al. 29 provide an overview of recent developments, particularly relating to hemodynamics applications with multi-scale aspects and demonstrates how they can be applied to problems such as atherosclerosis. Kim et al. 30 preformed simulations in arteries with varying pre-determined degrees of stenosis and investigated the reflection of the pressure pulse at the stenosis and how this changed as the level of stenosis becomes more severe. Wang et al. 31 performed LBM simulations in both a straight artery section and also one with an idealized stenosis. They considered the deformation of both sick and healthy red blood cells due to the hemodynamics and the artery geometry. The simulations were focussed at micro-channels, but could also be applied at the scale of the carotid artery, and would be particularly relevant when a high level of stensosis is present. Ge et al. 32 considered details of velocity and WSS fields for a number of pre-determined carotid artery geometries including an unstenosed artery. As well as focussing on the velocity and WSS from their patient specific arteries, they also considered aspects of parallel implementation required to provided a detailed 3D representation. Chen et al. 33 also considered the flow fields in a number of patient-derived artery geometries. In particular, they focussed on the differences between the flows in terms of regions or reversed flow as well as the time extent of both the reversed flows and secondary flows. Although they did not consider stenosis development, they considered these regions to be strongly correlated with the formation of plaques. Chen et al. went on to identify the features of the patient-specific geometries which led to the formation of the hemodynamic regimes associated with atherosclerosis.
Although there are many indications in the literature that hemodynamic factors are important and can act as markers for atherosclerosis, a recent review 34 has concluded that there is no consensus as to which of the markers in common use are best placed to identify early stages of atherosclerosis. Here we consider a numerical model which predicts early stenosis development based on a specified hemodynamic marker. Simulations of the blood flow, over a period of oscillation, are simulated using an LBM approach and the hemodynamic marker calculated. The stenosis model is then applied to increment the stenosis development, based on the marker. The hemodynamic parameter is then re-calculated from a further periods of LBM simulation. The stenosis development is modeled using the model using an approach 26 which used the TA near-wall velocity as the marker for stenosis development; here we expand this approach to consider a range of different markers which are commonly identified in the literature. In this way, a marker-specific development is realized, which enables comparison between the different markers and evaluation of how accurately the stenosis development is modeled with respect to observations and additionally to enable their influence on stenosis development to be assessed. A better understanding of the importance of different haeomodyamic markers will aid in the detection of early atherosclerosis onset.

METHODS
The blood flow is simulated using the LBM 19 and the development of the stenosis is modeled following. 26 These two numerical approaches are outlined below. Following References 25,26,35, and 36, the simulations will be performed in 2D and a rigid wall approximation will be applied. This allows the main features of the flow to be observed and gives a good approximation in the stenosed region where the walls are not compliant.

The lattice Boltzmann method
The LBM implemented here is based on the Bhatnagar-Gross-Krook (BGK) LBM equation 19,37,38 where i is an index which labels the link directions, e i , of the regular grid used in the simulations and is the relaxation time. Here we use the D2Q9 model which consists of a 2D grid where each node is connected to its four nearest neighbors: and its four next-nearest diagonal neighbors: The blood is represented by a distribution function f i (x, t) defined on each of the links, including e 0 = (0, 0). These distribution functions evolve according to equation (1), where the left-hand-side represents streaming on the gird and the right-hand-side simulates collisions in terms of a relaxation of the distribution functions toward their equilibrium values, f eq i . The equilibrium distribution functions are found from the BGK collision operator in terms of the local velocity, u and density as where w 0 = 4/9, w 1 = w 2 = w 3 = w 4 = 1/9 and w 5 = w 6 = w 7 = w 8 = 1/36. The blood density, velocity, and strain rate tensor, 39 S , can be calculated at each site from the distribution function The relaxation time in Equation (1) determines the kinematic viscosity of the fluid as The LBM model described by Equations (1)-(6) is referred to as the single-relaxation-time LBM, due to the single relaxation parameter, . This is the approach applied here, in common with other LBM studies relating to the carotid artery 29,31,32,36 due to its efficient computation; however, it is worth noting that it does have some disadvantages. In particular, it can experience stability issues at high Reynolds number flows, where a multi-relaxation-time (MRT) 40,41 model can be applied. This includes a number of free parameters which can be selected to improve accuracy and stability. In-flow boundary conditions were applied at the artery base by setting the distribution functions to their equilibrium values, Equation (4), for a velocity pulse adapted from Reference 42 and has previously been applied in References 25,26,43 and 44. The velocity was reduced in a linear fashion to zero in a small boundary layer at the artery wall. An outflow extrapolation scheme 45 was applied at the top of the artery. Both the in and out flow boundaries were a significant distance from the region over which the results are presented here. An extrapolation no-slip boundary condition 46 was applied at the artery walls. This is a second-order scheme which enables the artery geometry to be simulated to sub-grid accuracy. 26,43,47 Simulation parameters were selected to achieve a Reynolds number (based on boundary layer thickness and the peak velocity) of 307 and a Womersley parameter of 4.5. 42 This was achieved with a diameter of 36 lattice units, a period of 66,987 time-steps, and = 0.5045.

Stenosis growth model
The stenosis growth model was first applied in Reference 26 where a full description is given. Here we briefly present the key steps: 1. The LBM is used to simulate a single period of motion in the artery. During the period, the hemodynamic properties are monitored at the artery walls. 2. The site where the stenosis growth is to be modeled it determined based on the values of the markers. In Reference 26, this was selected as the site (mesh point) with the minimum TA velocity. Here we consider a range of markers, defined in Section 3. 3. At this site, the geometry of the boundary is augmented into the fluid. This is performed by moving the boundary point (point where the link cuts the boundary) for the selected site. This involves three parameters: (a) The augmentation direction. This is determined by finding the angle between each possible link and the artery wall. The link closest the wall normal is selected and the boundary point on this link is augmented along the normal.

4.
A further period is simulated with the next wall geometry and the process is repeated.
In Reference 26, it was shown that the model was independent of these augmentation parameters if the width was selected as 1 and the length was less than or equal to 0.5 lattice units when the TA velocity was used as a marker. Preliminary simulations suggested the same range applied to the other markers employed here and a value of 0.3 lattice units was selected for the augmentation length to allow for a safety margin.
In this model, the stenosis develops every period and simulates the sequential development, based on the selected marker, rather than the actual rate of development. Thus the motion of the wall in step 3b does not need to consider any transfer of the wall motion into the fluid, since this would not happen every period of oscillation.

RESULTS AND DISCUSSION
In Reference 26, the TA near-wall velocity was used to determine the wall position, as the end of each period, where stenosis development is modeled. This is defined as where u (1) is the near-wall velocity measured 1 lattice unit from the wall. 26 Here we consider how this, and a range of other local hemodynamic markers affect stenosis development. The additional selected criteria are: where u (1) t is the tangential component of u (1) ; RRT ∝ SI, defined as the fraction of the period during which u (1) is less than 1% of the peak pulse velocity; and The WSS is found as the magnitude of w where 36 S is found from Equation (5) where n is the normal to the boundary, determined as n ⋅ b = 0, where b is a local vector along the artery boundary. 26 The tangential component of u (1) is also found using b. Summation over repeated Greek indices is assumed in Equation (14) and the definition of TAD II . Although a non-Newtonian model is applied here, TAD II gives an indication of the variation in viscosity in a non-Newtonian fluid such as blood, with high viscosity corresponding to low TAD II . It also gives an alternative measure of the shear stress at the wall. The constant of proportionality for RRT was found by normalizing with respect to the value three diameters upstream of the bifurcation. 48,49

Stenosis development
The development of the stenosis is shown in Figure 1. Figure 2 shows the development in layers in terms of the normalized parameter T * s representing the number of sites converted from wet to dry, normalized by the length of each simulation: T s = 747. All the simulations, except for the one based in the OSI, show a similar final geometry for the stenosed artery, which is consistent with that reported in the literature. [50][51][52][53] The simulation based on OSI is presented for a shorter time, due to the stenosis developing around the limits of the buffer region. Despite the general similarity, there are a number of differences: both in the shape of the final stenosis; and also its development. For example, compared to the others, Figures 1(E) and 2(E) show the stenosis development extending further into the artery on the ECA and extending less far along the ICA. Figures 1(G) and 2(G) show that when TAD II is used as a marker the formation occurs fully on the ICA before switching to the ECA. This is in contrast to the other cases where the development in the ECA occurs predominantly toward the start of the process. To further investigate the differences between simulations, Figure 3 shows ΔT * s defined as T * sX − T * sY , where T * sX and T * sY are the values of T * s when the hemodynamic marker X = u (1) and Y ∈ {u (1) t , TAWSS, RFI, SI, RRT, TAD II , O:WS}. OSI was not considered here since the behavior in Figures 1(H) and 2(H) is very different to the other cases. Figure 3(A,B) indicates that the markers u (1) , |u (1) t | and TAWSS are all enabling the stenosis to form in a very similar manner. The other parameters (RFI, SI, RRT, TAD II , and O:WS) cause the stenosis to develop in significantly different manners, compared to u (1) . Figure 4 shows ΔT * s for these parameters where X = O:WS and Y ∈ {RFI, SI, RRT, and TAD II }. Considerable similarity is seen in Figure 4(C) and to a lesser extent in Figure 4(A). This indicates that as a marker for stenosis development O:WS and RRT act is a very similar manner, with RFI also producing a similar outcome. The results in Figures 1 to 4 indicating that using low time-averaged near-wall velocity, time-averaged tangential near-wall velocity, and TAWSS to determine how and where the stenosis develops all give a virtually identical outcome, and so can be used interchangeably as markers for stenosis formation. The results obtained using RRT and the ratio O:WS also showed considerable similarity to each other and so these two markers can also be applied interchangeably. These two sets of markers along with RFI, SI, and TAD II , although predicting different progression for the stenosis, all predict similar final geometries. OSI, on the other hand, predicts a significantly different process which is not consistent with observations in the literature. [50][51][52][53] This suggests that although the OSI can be considered important in determining where the stenosis will develop, it is not on its own a marker for stenosis development; and must be combined with other hemodynamic properties, such as TAWSS in a combined parameter such as O:WS. Since the OSI parameter does not predict a realistic stenosis, it will not be considered further.

Geometrical development
To further assess the effect of the different markers, Figure 5 shows the area converted from wet to dry on the ECA and ICA as a function of T * s . On both walls, the curves for u (1) , |u (1) t |, TAWSS are consistent with the first two being almost indistinguishable. The two curves for RRT and O:WS are also very consistent with each other. The curve for RFI is similar to the RRT/O:WS curves, while TAD II and SI curves are somewhat different to the others. The similarity between the development predicted by the u (1) , |u (1) t | and TAWSS markers show many similarities, such as the initial Another useful way to characterize the stenosis development is through the diameter of the artery. This is shown for a number of sections, shown in Figure 6(A), where the black geometry is for a healthy artery (T * s = 0) and the red geometry is a stenosed artery (T * s = 1) when u (1) is used as the marker. The change in diameter along the sections is shown in Figure 6(B-G), for the different markers. The curves further highlight the strong similarity of u (1) , |u (1) t | and TAWSS. RRT and O:WS are again shown to be similar; RFI is somewhat similar to RRT and O:WS in some regions, while in others difference can be observed. TAD II and SI show a different development. Development on the ECA affects the development on section S 1 , S 2 (lower stenosis development)and S 4 (upper stenosis development), while sections S 5 , S 6 are only effected by stenosis development in the upper region of the ICA. Section S 3 shows changes relating to the middle region of the ECA and the low region of the ICA. Despite the overall similarity on the predictions from the u (1) , |u (1) t | and TAWSS markers, the section also highlight some differences. For example, on sections S 1 and S 2 (Figure 6(B,C)), the TAWSS marker generally predicts a slightly earlier and more substantial development in the lower region of the ECA development, while on section S 4 (Figure 6(E)), the opposite is true in the upper region, particularly for T * s > 0.45. The main differences between the prediction of RRT and O:WS and on section S 3 (Figure 6(D)) where the RRT marker predicts a more gradual reduction in diameter up to T * s = 0.3, and on section S 5 (Figure 6(F)), corresponding to the upper ICA stenosed region, where the RRT marker predicts a larger reduction in the diameter up to T * s = 0.5. Differences between the RFI marker and the RRT/O:WS markers are also most evident in this region with RFI predicting a faster and larger development in this time frame. Figure 6 also displays the rather different predictions which arise from both the TAD II and the SI markers. For example on sections S 1 , S 2 and to a lesser extend S 4 , the SI marker predicts a smaller diameter (and so more stenosis development) on the ECA for the majority of the process, while the TAD II marker generally predicts less development. The opposite is observed on S 5 and S 6 . This is consistent with the observations from Figure 5.
The effect of the different marker on the wall hemodynamics is investigated over the time of stenosis development in Figures 7 to 9 for the ECA. The hemodynamic properties shown are u (1) in Figure 7, RFI in Figure 8, and OSI in  Figures 7(A,B), the same behavior is observed up to T * s ≃ 0.3; after this the behavior differs: in particular the low near-wall velocity region is maintained in Figure 7(B) for much longer and is still fairly prominent at T * s = 1. This is also evident when comparing Figures 8(A,B) and 9(A,B) where the regions of significant RFI are similar for T * s < 0.3; however for T * s > 0.3, there are significant differences with considerably more reverse flow when when O:WS is used as a marker-that is despite the marker acting to remove areas of low OSI (through enabling the stenosis to develop in such areas), albeit in conjunction with the TAWSS through O:WS. In each case, differences are also observed when RFI, SI, and TAD II are used as markers. Figure 7(C) shows similar behavior to Figure 7(B) in that the low velocity region is maintained longer than in Figure 7(A), but the averaged velocity alters more uniformly. Figure 7(D) shows slightly different behavior with the low averaged velocity region ending move abruptly at each y-position. Figure 7(E) also shows a different behavior, with the initial profile extending to around T * s = 0.6 before the low velocity region expands downstream. This is due to the stenosis only forming initially on the ICA, as seen in Figures 1(G), 2(G), and 5. Interestingly in Figure 1(G), the stenosis expands further downstream (y > 50mm) when TAD II is used as a marker, compared to Figure 1(A) when u (1) is used. This can be observed in Figure 7(E) where the low average velocity region also reaches beyond y = 50 mm, for T * s > 0.9. Here the low averaged velocity region appears to pierce through the higher velocity region (y> 50 mm). The observations for the RFI, SI, and TAD II markers are consistent with Figures 8 and 9, which show regions of reversible and oscillatory flow associated with low average-velocity region. Similar behavior was also observed for the ICA.

CONCLUSION
In terms of a marker for the stenosis development, the near-wall time averaged velocity u (1) , its tangential magnitude |u (1) t |, and the TAWSS all produce very similar stenosis development and the resulting hemodynamic wall parameters are affected in a very similar manner. Markers RRT and O:WS also produce stenosis development and the associated hemodynamic wall parameters which are similar to each other.
RFI, SI, and TAD II markers produce noticeable differences between each other and between the two sets of parameters described above (u (1) , |u (1) t |, TAWSS, and RRT, O:WS). In each case, the stenosis develops in a slightly different way and the near wall hemodynamic behavior is also different. Despite the differences in formation, the form of the stenosis on both the ECA and ICA is similar for each of the markers and consistent with observations in the literature. [50][51][52][53] This suggests that each of these hemodynamic properties are good markers for stenosis development.
In contrast, when OSI is used as a marker, the simulations presented here show the stenosis being formed in a significantly different manner which is not consistent with literature observations. Low near wall velocity, as measured through u (1) produced a realistic stenosis formation as did O:WS. Although differences were observed for these two markers, both produced results consistent with literature and so it is not possible to determine which marker is superior; however, the similarity suggests that considering OSI in conjunction with the TAWSS at least in terms of the O:WS parameter considered here, does not produce a significant advantage over using low TAWSS on its own, or indeed low near wall velocity.
Details of imaged stenosed geometries in the literature are varied and focus on particular stages of development in different patients. Since stenosis development is patient specific, there is no pathway for linking these images or knowing how a stenosis will develop. Using a stenosis growth model enables this evolution to be observed. Additionally it gives a method for evaluating and comparing the many different hemodynamic markers for stenosis growth which are discussed in the literature and which can be difficult to compare between studies due to the different artery models/images which are used. Here different markers have been compared in a systematic manner to determine the effect they have on stenosis development and how they change as the stenosis develops. This showed that some of the markers, such as u (1) , and TAWSS, produced very similar results, suggesting they could be used interchangeably, It was also shown that WSS and OSI, the two most commonly used markers, produced significantly different results and that using OSI as a marker did not produce a realistic outcome. While the final outcome is affected by both the marker and the choice of stenosis model, the results from this work suggest WSS and OSI cannot be considered as interchangeable markers and that a different marker, such as TAWSS or O:WS would be a better choice than OSI. A deeper understanding of the role of these hemodynamic markers is important when interpretation the results of patient-specific simulations in terms of the potential for the onset of atherosclerosis.
Further work is required to investigate wither alternative stenosis growth models identify the same differentiation between markers, and also to expand the simulations to 3D to investigate the effect of features such as swirl which can not be captured in a 2D model. The approach to stenosis modeling applied here could be applied to facilitate a development to 3D.

PEER REVIEW INFORMATION
Engineering Reports thanks Valentina Mazzi and other anonymous reviewers for their contribution to the peer review of this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon request.

CONFLICT OF INTEREST
The authors have no conflict of interest relevant to this article.