Baseline local hemodynamics as predictor of lumen remodeling at 1-year follow-up in stented superficial femoral arteries

In-stent restenosis (ISR) is the major drawback of superficial femoral artery (SFA) stenting. Abnormal hemodynamics after stent implantation seems to promote the development of ISR. Accordingly, this study aims to investigate the impact of local hemodynamics on lumen remodeling in human stented SFA lesions. Ten SFA models were reconstructed at 1-week and 1-year follow-up from computed tomography images. Patient-specific computational fluid dynamics simulations were performed to relate the local hemodynamics at 1-week, expressed in terms of time-averaged wall shear stress (TAWSS), oscillatory shear index and relative residence time, with the lumen remodeling at 1-year, quantified as the change of lumen area between 1-week and 1-year. The TAWSS was negatively associated with the lumen area change (ρ = − 0.75, p = 0.013). The surface area exposed to low TAWSS was positively correlated with the lumen area change (ρ = 0.69, p = 0.026). No significant correlations were present between the other hemodynamic descriptors and lumen area change. The low TAWSS was the best predictive marker of lumen remodeling (positive predictive value of 44.8%). Moreover, stent length and overlapping were predictor of ISR at follow-up. Despite the limited number of analyzed lesions, the overall findings suggest an association between abnormal patterns of WSS after stenting and lumen remodeling.

www.nature.com/scientificreports/ and 15 months after the endovascular procedure 16 . Furthermore, it was reported that, after 2 years, the occurrence of ISR in SFA is reduced and the lesions seem more related to the evolution of atheromatous disease rather than neointimal hyperplasia 17,18 . Besides common risk factors related to a patient's pathological characteristics, such as disease history, smoking, and diabetes mellitus, some lesion-specific factors, including a longer lesion length and smaller vessel diameter, are potential risks of ISR 19,20 . Furthermore, as shown in previous studies on different vascular regions, in particular coronary arteries, the implantation of a stent within an artery alters the local hemodynamics, resulting in stagnation, recirculating flow, and areas subjected to low wall shear stress (WSS) 21 , which may promote the ISR development 22,23 . Recently, computational fluid dynamics (CFD) tools have been used to investigate the relationship between altered hemodynamics and restenosis in femoropopliteal artery atherosclerotic lesions treated with angioplasty or stenting 24 . A predictive statistical model of the treatment outcome at 6-month, which includes information about the artery hemodynamics immediately after intervention, was successfully proposed. However, a direct, local comparison between the luminal distribution of WSS-based descriptors and the lumen remodeling was not performed.
Accordingly, the present work aims to expand the previous analysis by investigating the impact of local hemodynamics on the vessel lumen remodeling, focusing on human SFAs treated with self-expandable stents. Specifically, patient-specific CFD simulations were performed in stented SFA models reconstructed from computed tomography (CT) images acquired at the first follow-up after stent implantation (i.e. 1-week follow-up). The local hemodynamics computed at 1-week follow-up was then related to the lumen remodeling occurring at 1-year follow-up. The analysis was conducted both at the global and local levels, taking into consideration demographic and clinical information, such as age, length of stented region, stent overlapping, and failure at 2-year follow-up, in addition to the hemodynamic descriptors and the lumen remodeling at 1-year follow-up. Figure 1 shows the general workflow applied to analyze the relationship between local hemodynamics and lumen remodeling in the stented region of patient-specific SFA models. The following sections detail the available clinical dataset, the methods adopted to reconstruct the SFA models, perform the CFD simulations and analyze the results.  25 . The 1-week model is used to perform computational fluid dynamics (CFD) analysis, prescribing as inlet boundary condition the velocity waveform derived from the patient's Doppler ultrasound (US) image. Finally, the local hemodynamics is compared to lumen remodeling.

Methods
1-year post-operative CT scan protocol were recruited between 2007 and 2012 at the Malcom Randall VA Medical Center (Gainesville, FL, USA). Within the obtained clinical dataset, 7 patients, for a total of 10 stented lesions, presented all the imaging data necessary to perform the present study. In particular, each patient presented both CT and Doppler ultrasound (DUS) images at 1-week (baseline) and 1-year post-operative follow-up. Furthermore, the dichotomous outcome of the intervention (i.e. success/failure) was also available at 2-year follow-up, even though not supported by imaging data. The failure was defined when a new procedure of revascularization was required.
Baseline patient demographics, as well as medical history, are summarized in Table 1. The patients were all male. Their age was 62.3 ± 6.6 years (mean ± standard deviation). All lesions (A to K, Table 1) were treated with the EverFlex self-expanding stent (EV3, Medtronic, Dublin, Ireland). The stent length was 113.8 ± 38.6 mm, with five cases presenting overlapping devices. The stent overlapping was observed and defined by visual inspection of CT images and when two devices were overlapped, they were counted as a unique stent (stented region length 165.0 ± 99.6 mm). According to the 2-year follow-up, 5 out of 10 lesions failed.

Human subjects/informed consent statement. This study was approved by the Institutional Review
Board at the University of Florida and conformed to the Helsinki Declaration on human research of 1975, as revised in 2000. Written informed consent was obtained from the patients. No animal studies were carried out by the authors for this article.
Three-dimensional vessel reconstruction and computational analyses. Three-dimensional (3D) SFA geometrical models were reconstructed both at baseline and 1-year follow-up (Fig. 2) using a previously validated semi-automatic methodology 25 . Briefly, the CT images were segmented by means of an active contour method, based on a level set algorithm. Both the stented and the non-stented regions were segmented through calibrated thresholds. Calcifications, metallic artifacts, and thrombi within the stent were automatically removed. The common femoral bifurcation was included in all vessel reconstructions to take into account its impact on the local hemodynamics 25 .
The SFA geometrical models were discretized using tetrahedral elements with curvature-based refinement and 5 boundary layers of prismatic elements near the wall. The element size was defined according to a previously conducted mesh independence study 25 , resulting in computational grids ranging from 2,202,629 to 8,535,680 elements (SFA models with lesion B and K, respectively). The software ICEM CFD (v. 18.2, Ansys Inc., Canonsburg, PA, USA) was employed for the discretization process. Transient CFD simulations were carried out on the baseline SFA models using Fluent (v. 18.2, Ansys Inc., Canonsburg, PA, USA).
Regarding the boundary conditions, a velocity waveform derived from a patient's DUS image was imposed at the inlet 25 . A flow-split of 67% and 33% was applied to the outlet of the SFA and the profunda femoral artery, respectively 25 . The no-slip condition was prescribed at the wall boundaries, assuming the vessel walls as rigid. Blood was modelled as a non-Newtonian fluid using the Carreau model, with a density of 1060 kg/m 3 . Details about the simulation settings are presented elsewhere 25 .
Analysis of the results. Comparison between local hemodynamics and lumen area change. The local hemodynamics of the baseline SFA models was analyzed by computing three WSS-based descriptors quantifying low and/or oscillatory WSS, previously indicated as promoting factor of ISR 26 . In detail, the time-averaged WSS (TAWSS), the oscillatory shear index (OSI) and the relative residence time (RRT), were computed as follows: www.nature.com/scientificreports/ where WSS is the WSS vector and T is the duration of the cardiac cycle. As previously conducted 27 , the WSS-based data of all lesions were combined to define objective thresholds of disturbed shear stress: the 33th percentile was identified for TAWSS, the 66th percentile for OSI and RRT. The following thresholds were found: 1.05 Pa, 0.20 and 1.68 Pa −1 for TAWSS, OSI and RRT, respectively. Then, the percentage of surface area exposed to TAWSS (OSI and RRT) lower (higher) than the thresholds was computed and identified as TAWSS33 (OSI66 and RRT66).
The 3D distributions of the WSS-based descriptors along the stented portion (i.e. vessel region of interest) were also re-organized into two-dimensional (2D) maps by means of the open-source software VMTK (Orobix, Bergamo, Italy) (Fig. 3A). The 2D maps were discretized into cells of 1 mm in the axial direction and 1° in the circumferential direction. The 2D maps were then circumferentially averaged (Fig. 3A) to obtain one-dimensional (1D) maps with axial resolution of 1 mm, being conservative with respect to the CT axial resolution (lower than 1 mm in each acquisition, range = [0.87 mm, 0.98 mm]). The circumferentially averaged WSS-based data allowed the match between the local hemodynamic results at baseline and the data of lumen remodeling at 1-year post-intervention.
The lumen remodeling was quantified by calculating the lumen area change, intended as the difference between the cross-sectional lumen area at baseline and the one at 1-year follow-up. A positive (or negative) lumen area change value indicates that the remodeling at 1-year takes place in the direction of lumen area reduction (or increase). To do that, from the stented region of the 3D reconstructed vessels subdivided into 0.2 mm axiallyspaced cross-sections were extracted using an automatic algorithm developed in Grasshopper (v. 5, Rhinoceros, Robert McNeel & Associates, Seattle, WA, USA). The lumen area was averaged every 1 mm in the axial direction (i.e. considering 5 cross-sections per millimeter) to obtain 1D maps with the same axial resolution as that of the WSS-based descriptors.
Independence of the hemodynamic results. The statistical analysis of the local association between the WSSbased descriptors and lumen remodeling requires that the data points are spatially decorrelated 28,29 . Thus, a spatial autocorrelation study was first performed to ensure the independence of the hemodynamic data points. This analysis provided a decorrelation length (L decorr , Fig. 3B), namely the minimum distance between the hemo- www.nature.com/scientificreports/ dynamic data points so that they can be considered as independent 28,29 . The analysis was based on the definition of an index, called z-score, derived from that defined by Peiffer et al. 28 . In particular, referring to each segmentspecific result along the 1D map as a cell, for each i-cell the z-score was computed as: www.nature.com/scientificreports/ where i is the cell, x i the central value of the hemodynamic variable, µ the mean and σ the standard deviation of the considered cells. This index was computed for each cell at different distances from the target cell, so that for each cell it is given Z 1,i , Z 2,i ,…, Z n,i where 1, 2, …, n defines the distance between the target cell and the inspected cell ( Fig. 3B, left panel). For instance, for n = 1, Z 1,i , µ 3 and σ 3 were evaluated only on 3 cells, for Z 2,i on 5 cells and so on. According to its definition, the z-score is greater when the value of the target cell is more distant from the mean value of the inspected cell i. Once computed for each cell, the resulting z-score was averaged based on the distance, obtaining an average Z 1 , Z 2 ,…, Z n . The decorrelation length was then defined as the first peak in the diagram 'averaged z-score-cell distance' (Fig. 3B, central panel) 28,29 . Knowing the decorrelation length, different combinations of independent data points can be created. An average of the combinations was then calculated, thus obtaining a single 1D map (Fig. 3B, right panel), namely an array of the results for each lesion and for each variable of interest, which was used in the subsequent statistical analyses. Data were presented as either mean ± standard deviation or median (interquartile range), depending on the distribution. The normality of the distributions was tested by means of Kolmogorov-Smirnov test. In case of a low number of data points (e.g. global analysis with 10 data points), non-parametric tests were employed. To measure the strength and direction of the association existing between the hemodynamic descriptors and the lumen area change, the Spearman's rank-order correlation coefficient was considered in case of monotonic distribution both at global and local, individual level. Furthermore, similarly to previous studies 30,31 , the positive predictive value (PPV) was calculated for each hemodynamic descriptor to evaluate the probability that vessel regions exposed to TAWSS33 (OSI66 and RRT66) can successfully identify corresponding regions with high lumen area change. In this analysis, the second tertile of the lumen area change distribution was used to discriminate between low/ high lumen area change.

Levels of investigation
Besides the hemodynamic descriptors, the clinical information was taken into consideration. Logistic regression was performed to relate independent continuous variables to the dependent dichotomous (success-failure) variable (Tjur's pseudo R 2 was considered). Furthermore, a two-sided Fisher's exact test was performed between the dichotomous variables stent overlapping and failure at 2-year follow-up.
All the statistical analyses were conducted using SPSS Statistics (v. 25, IBM, Armonk, NY, USA) and Graph-Pad Prism (v. 8.3.1, GraphPad Software, San Diego, CA, USA). A two-tailed, p-value < 0.05 was considered to be statistically significant. Table 2 summarizes the lesion-specific decorrelation lengths obtained from the spatial correlation analysis (range = [2 mm, 6 mm]). For each lesion under consideration, the computed decorrelation length was different according to the hemodynamic descriptors and the most conservative value (i.e. the largest L decorr ) was chosen for the subsequent statistical analyses.

Spatial correlation analysis.
Global analysis. In addition to the decorrelation lengths, Table 2 reports the median values of the hemodynamic descriptors and the percentage surface areas exposed to disturbed shear computed from the decorrelated dataset of each lesion. These lesion-representative values were used in the global analysis of the results. Table 2. Results of the spatial correlation analysis in terms of decorrelation length (L decorr ) and the number of decorrelated data points that are considered for the subsequent statistical analysis (N points ). Median values (interquartile range, IQR) of time-averaged wall shear stress (TAWSS), oscillatory shear index (OSI), relative residence time (RRT), and lumen area change (ΔA), and percentage surface areas exposed to TAWSS33, OSI66 and RRT66 for each lesion under investigation. www.nature.com/scientificreports/ A statistically near-significant negative correlation was found between the median values of TAWSS and lumen area change (ρ = − 0.62, p = 0.06), suggesting that lower values of TAWSS are associated with a higher lumen area reduction. This association was statistically significant when considering the mean values of TAWSS (ρ = − 0.75, p = 0.013). Furthermore, a high positive correlation was present between the percentage surface area exposed to TAWSS33 and the lumen area change (ρ = 0.69, p = 0.026). This result indicates that larger surface areas exposed to low TAWSS are associated with a higher lumen narrowing. On the contrary, no significant associations emerged between the median values of OSI and lumen area change (ρ = − 0.40, p = 0.25) as well as between the percentage surface area exposed to OSI66 and lumen area change (ρ = − 0.15, p = 0.69). Similarly, the median values of RRT and the percentage surface area exposed to RRT66 were not significantly correlated to the lumen area change (ρ = 0.50, p = 0.143 and ρ = 0.44, p = 0.21, respectively).
Introducing the patient-specific characteristics listed in Table 1, a logistic regression (both simple and multiple) was performed to explore the relationship between the continuous variables length of the stented region (measured from the CT scans), age and hemodynamic descriptors, and the dichotomous systemic failure at 2-year post-intervention. From the simple logistic regressions, a complete separation (also called perfect prediction, i.e. the outcome variable separates completely the predictor variable) emerged between the length of the stented region with a cut-off value of 120 mm and the failure at 2-year post-intervention (Fig. 4A). The age was found not to predict the dichotomous failure (p = 0.26, R 2 = 0.15) (Fig. 4B). Regarding the hemodynamic descriptors, the TAWSS presented a nearly statistically significant association with failure (p = 0.06, R 2 = 0.28). The association between neither the OSI nor the RRT and the failure at 2-year post-intervention was significant (p = 0.11, R 2 = 0.22 for OSI; p = 0.98, R 2 < 0.0001 for RRT) (plots not shown). From the multiple logistic regression, in which all the three hemodynamic descriptors were introduced as predictors with respect to the failure at 2-year follow-up, the analysis did not show a significant difference between the success and failure groups (Fig. 4C), with a Tjur's R 2 equal to 0.49 (p = 0.033). Furthermore, the selected model resulted correct, according to the Hosmer-Lemeshow test (with a value of 6.83 and p = 0.55). The association between the stent overlapping and the failure was statistically significant (p = 0.008), indicating that the stent overlapping is a predictor of the treatment outcome (Fig. 4D).
Local analysis. Regarding the local, individual level of analysis, no significant correlation between the hemodynamic descriptors and the lumen area change was found for most of the cases (Table 3). In two cases only, the lumen area change presented a moderate correlation with TAWSS (case A: ρ = − 0.38, p = 0.04; case H: ρ = − 0.64, p = 0.02). The two lesions presenting the lowest and the highest correlation coefficient between TAWSS and lumen area change (i.e. lesions D and H, respectively) were also taken as examples to show the contour maps of TAWSS along the stented region at baseline and the corresponding vessel geometry at 1-year follow-up as well as the distributions of the variables of interest (Fig. 5A). The two cases were characterized by statistically different distributions of TAWSS (p < 0.0001), RRT (p = 0.0005) and lumen area change (p < 0.0001), according to the Mann-Whitney rank test. Conversely, no significant differences were found for the distributions of OSI of the two lesions. Lesion H, which presented lower values of TAWSS, showed a pronounced lumen narrowing as compared to lesion D.
Regarding the local, collective level of analysis, the Spearman's correlation was performed to examine the relationship between the hemodynamic descriptors. A negligible negative correlation was found between TAWSS and OSI (ρ = − 0.21). On the contrary, as expected, significant high correlations were present between TAWSS and RRT (ρ = − 0.69, p < 0.0001) and between OSI and RRT (ρ = 0.77, p < 0.0001). Spearman's correlation and linear regression analyses between the hemodynamic descriptors and the lumen area change revealed that the strongest association was present between TAWSS and lumen area change, as observed in the global analysis ("Global analysis" section). In particular, a negative correlation was found (ρ = − 0.28, p < 0.0001). Moreover, the regression line (Fig. 5B, left panel) had a negative slope, indicating that for higher values of TAWSS, the lumen narrowing reduces. A weak association emerged between OSI and lumen area change (ρ = − 0.15, p = 0.005), whereas no significant relationship was found between RRT and lumen area change (ρ = 0.07, p = 0.184). Finally, to evaluate the predictive power of the hemodynamic descriptors, the PPV was quantified. The PPV was 44.8%, 28.6% and 35.9% for the TAWSS, OSI and RRT, respectively.

Discussion
ISR in SFAs is a complex multi-factorial process widely affecting the outcome of the stenting procedure 32 . Among the different factors promoting ISR, this study investigated the impact of the altered hemodynamics after stent implantation on vessel lumen remodeling. Specifically, a computational approach, based on the 3D reconstruction of patient-specific SFA geometrical models and on CFD simulations 25 , was adopted to relate the local hemodynamics at baseline, expressed in terms of WSS-based descriptors, with the lumen remodeling at 1-year follow-up. The findings of this study suggest that in human SFA lesions an association between the TAWSS at baseline and the lumen remodeling at 1-year follow-up exists. Specifically, (i) lower baseline values of TAWSS were associated with higher lumen area reduction, and (ii) larger surface areas exposed to lower baseline TAWSS were associated with higher lumen area reduction. On the contrary, no statistically significant correlations were present between the other analyzed hemodynamic descriptors (i.e. OSI and RRT) and lumen area change.
Although computational approaches have been previously applied to investigate the influence of local hemodynamics on ISR in other vascular regions, in particular in coronary arteries 23,33,34 , their application to lowerlimb peripheral arteries has been limited. Until now, a few studies have focused on the influence of curvature and tortuosity on atherosclerosis 35,36 , or on the effect of wall motion 37  www.nature.com/scientificreports/ and relate the WSS-based descriptors to the occurrence of restenosis at 6-month after angioplasty or stenting procedure. A logistic regression model including a combination of different hemodynamic descriptors and the treatment method was built showing an accuracy of 80% in predicting the 6-month outcome. However, none of the WSS-based descriptors was able to explain the occurrence of restenosis alone. The imaging modality used for the vessel reconstruction (i.e. 2D angiography) and the lack of patient-specific boundary conditions may have affected the hemodynamic results in that study. More importantly, the analyzed arteries were only classified in  www.nature.com/scientificreports/ two groups, namely non-restenosed and restenosed arteries, based on the 6-month outcome but a local quantification of the lumen remodeling was not conducted. On the contrary, the present study provided a local quantitative comparison between the hemodynamic descriptors and the lumen area change focusing exclusively on SFA segments treated with self-expanding stents. The use of CT images for the 3D vessel reconstruction allowed obtaining more accurate SFA geometries, being at the same time less invasive. Furthermore, patient-specific boundary conditions, derived from DUS images, were here applied to the CFD models. Before investigating the link between local hemodynamics and lumen remodeling, a preliminary spatial correlation analysis was carried out to eliminate the data points that could potentially be auto-correlated or mutually influenced 29,40 . Indeed, several statistical approaches are commonly used under the assumption that the measured outcomes are independent of each other. Nonetheless, it is very liable that, especially in spatial data, some or all outcome measures exhibit spatial autocorrelation 29,40 . A previously described methodology 28 was here adapted to circumferentially averaged 1D data in order to create an ensemble of decorrelated and independent data points. The spatial correlation analysis highlighted that the minimum distance at which two data points could be considered as decorrelated is lesion-specific and should be evaluated case-by-case. The main explanation for this could be that the hemodynamic descriptors are directly related to the geometrical shape of the specific arterial model.
To explore if a correlation exists between each hemodynamic descriptor and the lumen area change, the statistical analyses were conducted both at a global and local level. The global findings, obtained considering www.nature.com/scientificreports/ one representative value for each lesion, were partially confirmed by the local analyses. In particular, also at the local, collective level, a negative correlation was found between the TAWSS and the lumen area change. The lower correlation coefficient as compared to that found in the global analysis may be explained by the higher degree of uncertainty and the inaccuracy that is introduced when a deeper, local investigation level is reached, as observed in previous works 41,42 . At this level of analysis, the total error affecting the statistical results was given by the sum of multiple factors, including (i) the 3D vessel reconstruction error, associated to the limited CT resolution, (ii) the uncertainty related to the inlet boundary condition, due to its derivation from DUS images, (iii) the computational error and (iv) the approximations related to the post-processing of the results, in particular the circumferential averaging of the results, which resulted in 1D maps that did not consider how the hemodynamic descriptors and the lumen remodeling varied circumferentially along each segment of interest. Despite the previous considerations, the local level of analysis enabled the establishing of the PPV for the lumen remodeling for each hemodynamic descriptor. While the PPV was low for the high OSI and RRT (28.6% and 35.9%, respectively), a PPV of 44.8% was found for the low TAWSS. This predictive value for the low TAWSS seems promising, indicating that the low TAWSS is a good predictor for local lumen narrowing. In fact, although a direct comparison of this index with other computational studies on SFAs and ISR is not possible, previous computational hemodynamics studies focusing on plaque progression in coronary arteries reported a lower 31,43 or similar PPVs 30 .
Additionally, this study provided insights into the relationship between the hemodynamic descriptors, some demographic and clinical variables (i.e. patient's age, length of the stented region, and stent overlapping) and the success or failure of the endovascular procedure at 2-year follow-up. No significant relationship was found between patients' age and treatment failure at 2-year post-intervention, in accordance with a previous work 44 . Conversely, the length of the stented region and the presence of stent overlapping were predictors of the treatment failure, coherently with previous findings in other vascular segments 45 .
Some limitations might affect the findings of this study. The statistical power of the used models was limited by the dimension of the clinical dataset. In particular, at the global level only 10 points per analysis were considered, possibly affecting the statistical power of the used models and the significance of the correlations. However, a local, collective level of analysis was also taken into account. At this level, the analyses counted more than 350 independent points, and the resulting findings were aligned with those of the global analyses, thus supporting and extending that level of investigation. Moreover, similarly to a previous study 24 , the SFA models did not include the stent geometry because of the limited resolution of the clinical images used for the 3D reconstruction, which did not allow detecting the stent struts. Although the effect of stent struts, stent overlapping and malapposition on hemodynamics is not taken into account, the relationship between hemodynamics and lumen remodeling emerged here. Finally, the local hemodynamics of the stented SFAs was investigated under the assumption of rigid walls and straight leg configuration, assuming that the leg movement did not markedly affect the hemodynamic results of the SFA segment. Future patient-specific studies may include the effect of lower limb movement by adopting the methodology recently implemented in an idealized moving-boundary CFD model of femoropopliteal artery 39 .

Conclusions
In this study, the impact of local hemodynamics on lumen remodeling was evaluated in human SFA lesions treated with self-expanding stents. The overall findings showed an association between abnormal patterns of WSS at baseline and lumen narrowing at 1-year follow-up. Specifically, lower TAWSS values and larger areas exposed to low TAWSS were associated with higher lumen area reduction. Conversely, OSI and RRT did not show any significant correlation with the lumen area change. Furthermore, the low TAWSS was a better predictive marker (i.e. higher PPV) of lumen remodeling as compared to OSI and RRT. In addition to the hemodynamic analysis, clinical variables such as the length of stented region and the presence of stent overlapping in those lesions treated with multiple stents were recognized as predictors of the treatment failure at 2-year follow-up. Despite these promising results, these associations need to be confirmed on a wider clinical dataset, with the final goal of defining stronger clinical predictors for ISR in SFAs. In this perspective, the presented workflow might serve as a clinical tool to predict the outcomes of a stenting procedure in SFA based on the hemodynamics analysis.