Validation of ultrasound velocimetry and computational fluid dynamics for flow assessment in femoral artery stenotic disease

Abstract. Purpose To investigate the accuracy of high-framerate echo particle image velocimetry (ePIV) and computational fluid dynamics (CFD) for determining velocity vectors in femoral bifurcation models through comparison with optical particle image velocimetry (oPIV). Approach Separate femoral bifurcation models were built for oPIV and ePIV measurements of a non-stenosed (control) and a 75%-area stenosed common femoral artery. A flow loop was used to create triphasic pulsatile flow. In-plane velocity vectors were measured with oPIV and ePIV. Flow was simulated with CFD using boundary conditions from ePIV and additional duplex-ultrasound (DUS) measurements. Mean differences and 95%-limits of agreement (1.96*SD) of the velocity magnitudes in space and time were compared, and the similarity of vector complexity (VC) and time-averaged wall shear stress (TAWSS) was assessed. Results Similar flow features were observed between modalities with velocities up to 110 and 330  cm/s in the control and the stenosed model, respectively. Relative to oPIV, ePIV and CFD-ePIV showed negligible mean differences in velocity (<3  cm/s), with limits of agreement of ±25  cm/s (control) and ±34  cm/s (stenosed). CFD-DUS overestimated velocities with limits of agreements of 13±40 and 16.1±55  cm/s for the control and stenosed model, respectively. VC showed good agreement, whereas TAWSS showed similar trends but with higher values for ePIV, CFD-DUS, and CFD-ePIV compared to oPIV. Conclusions EPIV and CFD-ePIV can accurately measure complex flow features in the femoral bifurcation and around a stenosis. CFD-DUS showed larger deviations in velocities making it a less robust technique for hemodynamical assessment. The applied ePIV and CFD techniques enable two- and three-dimensional assessment of local hemodynamics with high spatiotemporal resolution and thereby overcome key limitations of current clinical modalities making them an attractive and cost-effective alternative for hemodynamical assessment in clinical practice.


Introduction
Peripheral arterial disease (PAD) is a major cardiovascular disease burden with a global prevalence of 10% to 25% in age groups above 65 years, 1 with the femoral artery as a key site of atherosclerotic manifistations. 2It is well known that areas with complex flow structures, such as recirculation zones, induce low and oscillatory wall shear stress (WSS), thereby forming a local risk factor for atherosclerosis. 3Flow disturbances induced by the femoral bifurcation, the triphasic flow waveform including a pronounced backflow phase, 4 and the confined anatomic space of the superficial femoral artery (SFA) 5 may be the leading causes for the high prevalence of atherosclerotic lesions in this region.The association between adverse flow patterns and the onset and progression of atherosclerosis in patients has so far not been studied in-depth, mainly due to a lack of cost-effective and validated in-vivo flow measurement techniques.
The current diagnostic workflow to examine blood flow in PAD patients involves duplex ultrasound (DUS), a combination of conventional B-mode ultrasound and pulsed-wave Doppler measurements.][9] Besides DUS, phase contrast MRI (PC-MRI) can also be used to study time-resolved threedimensional hemodynamics.However, PC-MRI has limited availability, is expensive, time consuming, and has a low spatiotemporal resolution.Moreover, it averages the phase information over multiple cardiac cycles thereby losing insights in cycle-to-cycle variations. 10,11ultiple innovative approaches, using a variety of imaging modalities, have promising and synergistic potential for blood flow quantification in patients.First of all, several ultrasound techniques, such as multi-directional Doppler imaging and tracking of blood speckle or ultrasound contrast agents (UCA), can be used to measure two-dimensional local hemodynamics. 12,13dvances in plane wave ultrasound imaging improved the acquisition process by enabling a framerate of over 1000 frames per second (fps), thereby vastly overcoming the limited temporal resolution of PC MRI.For complex flows in stenosed vessels, speckle tracking presents a key advantage: it does not suffer from a tradeoff between tracking low and high velocities, in contrast to Doppler-based imaging where the pulse repetition frequency influences its sensitivity to certain velocities due to aliasing.Tracking of UCA, referred to as ultrafast ultrasound particle image velocimetry or echoPIV (ePIV), has been validated in vitro against the gold standard, laser optical PIV (oPIV), in various anatomically realistic models, such as the left ventricle 14 and the carotid bifurcation. 15Furthermore, its capability to measure blood flow in stented femoral arteries of patients with PAD has recently been demonstrated. 16However, these in vitro and in vivo studies did not investigate significant transitional flows, i.e., the flows with mixed laminar and turbulent characteristics that could arise downstream of a stenosis.The accuracy of ePIV in stenotic lesions, which comprise large velocity gradients and many small-scale vortices, is thus yet to be demonstrated.
Alternatively, computational fluid dynamics (CFD) modeling in patient-specific geometries can be used to simulate three-dimensional velocity profiles and associated flow parameters. 17FD has shown its benefit across the field of cardiovascular medicine to study complex physiological flow and to calculate different hemodynamical parameters, such as stenotic pressure gradients 18 and WSS. 19For pressure gradients in stenosed coronary arteries, extensive clinical validation has been performed, 20 but the validation of velocity fields computed by CFD in pulsatile stenotic flows is limited; one study in particular has demonstrated the ability of a spectralelement method for simulating the complex pulsatile transitional flow in a stenosis. 21omparisons of this direct numerical simulation to Reynolds-averaged Navier Stokes (RANS) turbulence models in CFD highlighted the shortcomings of these models when assessing the transitional flow behavior.Most importantly, these RANS models provided erroneous predictions in the onset of turbulence and the relaminarization length. 21An under-resolved CFD simulation without an explicit turbulence model or a large eddy simulation model, where van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . . the computational cost is reduced by implicitly or explicitly modeling turbulence for only the small, subgrid, length scales (instead of all length scales with RANS models), is likely better suited for simulating stenotic flows in patients. 22econd, CFD results critically depend on the applied boundary conditions, which consist of the segmented geometry and the flow conditions enforced on the inlets and outlets.Most validation studies have applied optimal boundary conditions as they were based on the experimental technique that it was compared to. 23,24Thus, CFD accuracy with current clinical techniques is directly limited by the spatial resolution of the applied imaging modality.Flow information for setting the boundary conditions is typically only available from DUS.
Given the higher resolution, absence of angle-dependency and a 2D or 3D blood velocity quantification, ePIV and CFD-or a combination of both-might provide novel insights in local hemodynamics.This provides the ability to assess new parameters of flow disturbances in patients, such as WSS and vector complexity (VC).VC has recently been used to evaluate lesion severity in SFA stenoses and outperformed the conventional velocity ratio obtained with DUS. 25 As such, ePIV and CFD might improve the diagnosis of PAD and provide a more patient-specific and durable treatment.However, the implementation of ePIV and CFD in clinical practice requires validation in clinically relevant flow conditions.This study investigates the accuracy of ePIV and CFD with ePIV and DUS-derived boundary conditions for determining spatiotemporal velocity vector fields in a non-stenosed and stenosed femoral bifurcation model by comparison to oPIV as a ground-truth reference technique.

Femoral Bifurcation Model
The geometry of the femoral bifurcation was based on a literature review (Appendix A) and was constructed as follows: common femoral artery (CFA) diameter of 8.9 mm, SFA diameter of 6.2 mm, and deep femoral artery (DFA) diameter of 6.1 mm.To facilitate imaging, the bifurcation was kept in a single plane with a 37.5-deg angle between the SFA and DFA.Two models were made, one without stenosis [control model, Fig. 1(a)] and one with a 75%-area concentric stenosis of the CFA [stenosed model, Fig. 1(b)].
Since no single available material was satisfactorily suited for both oPIV and ePIV, twin models were produced, one for oPIV and one for ePIV.The oPIV models were fabricated using the lost core technique: the geometries were first printed (Ultimaker S5, Ultimaker, Geldermalsen, The Netherlands) with hollow acrylonitrile butadiene styrene (ABS) and, after acetone vapor smoothing, casted in polydimethylsiloxane (PDMS, Sylgard 184, Dow Corning, Midland, Michigan).After curing, the ABS parts were manually extracted.Residual ABS that could not be removed was dissolved by running an acetone solution through the flow lumen.The ePIV models were constructed by 3D-printing the geometries with a flexible resin (Flexible 80A; Form 3, Formlabs, Somerville, Massachusetts) with 50 μm resolution and 1 mm wall thickness.The flexible resin had a speed of sound that matched that of the blood mimicking fluid (BMF) and had low enough attenuation that allowed for UCA visualization and flow quantification.
To assess the geometric similarity of the established lumens in the optical and acoustical twin model, a cone-beam computerized tomography (CT) scan (Artis Pheno, Siemens AG, Berlin, Germany, spatial resolution 0.47 mm) was made.The geometry was segmented with a subpixel-accurate level-set approach.Subsequent analysis of the centerline and inscribedsphere based radius of the models showed a comparable SFA-DFA angle (oPIV model 37.8 deg; ePIV model 37.4 deg) and deviations in the radius of the CFA and SFA in the order of 0.25 mm (Appendix B).

Experimental Flow Setup
The models were placed in a circulating flow setup, as schematically presented in Fig. 1(c).A triphasic flow profile, derived from DUS measurements in the CFA of seven healthy volunteers (four male, three female; aged 20 to 30 years), was replicated at a frequency of 1 Hz, i.e., 60 beats per minute (Appendix A).A hydraulic piston pump (SuperPump, ViVitro labs, Victoria, Canada) was used to enforce time-varying forward flow.The downstream vasculature of both the SFA van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .and DFA was modeled with a tunable Windkessel model (air chamber) and tunable resistance.An electronically controlled valve was used to reproduce the backflow phase, during which the pump valve is closed and the built-up Windkessel pressure drives flow in retrograde direction through the opened electronic valve (electronically synchronized with the piston pump).Ultrasonic transit time flow sensors in the setup (Sonoflow CO.55, Sonotec, Halle, Germany) were used to establish target mean flow rates by tuning the piston amplitude and resistances.The amplitude of the piston was set such that a flow volume of 9.6 mL was ejected during forward flow, of which 2.3 mL returned during the backflow phase.The SFA and DFA resistances were tuned to enforce a mean 44%:56% flow split (Appendix A).A mixture of 44 w% water, 34.5 w% glycerol, and 21.5 w% urea was used as BMF, which had an optical refractive index matched to PDMS. 26Density of the BMF was measured at 1.144 kg∕m 3 with a density meter (Anton Paar GmbH, Austria) and viscosity at 4.17 mPa.s with an Ubbelohde viscometer (Xylem Analytics, Germany), which comply with the physiological range of blood.

Optical Particle Image Velocimetry
Fluorescent spheres (rhodamine-coated; size, 1 to 20 μm; density, 1190 kg∕m 3 ; Dantec Dynamics A/S, Skovlunde, Denmark) were added to the BMF as tracer particles and illuminated with a thin laser sheet (sheet thickness maximum: 0.8 mm, 5-W DPSS laser, 532 nm; Cohlibri, Lightline, Germany) that was passed through the center of the bifurcated flow lumen.A highspeed camera (FASTCAM SA-Z; Photron Inc., West Wycombe, Buckinghamshire, United Kingdom) was positioned perpendicular to the laser sheet and used to capture 10 cardiac cycles of the bifurcation (55 × 55 mm 2 ) at 8000 frames per second (fps).
Data analysis was performed in MATLAB R2022b (The Mathworks, Natick, Massachusetts).The images were first filtered with a moving minimum intensity projection background subtraction (window length = 500 frames) and preprocessed with contrast-limited adaptive histogram equalization.Subsequently, the 2D velocity vector fields were derived through PIV analysis using the PIVlab toolbox, which tracked velocities in a static (Eulerian) grid of rectangular image windows. 27Three iterations of block wise normalized cross-correlation were applied with progressive grid refinement.Square interrogation windows of 64, 32, and 16 pixels with 50% overlap were used for the successive iterations, yielding a final isotropic vector resolution of 0.42 mm.Final postprocessing included a 10-ensemble temporal averaging filter, resulting in 400 vector fields per second.

Echo Particle Image Velocimetry
The ePIV measurements were performed with a Vantage 256 ultrasound system connected to a linear transducer (L11-4v; Verasonics, Kirkland, Washington).The printed geometries were encased in a reservoir and submerged in BMF.The reservoir contained an acoustical window to which the transducer was fixed to position the center of the lumen in the elevation focus, at a depth of 15 to 25 mm.As such, the alignment between the oPIV and ePIV measurements was assured.A UCA (Sonovue microbubbles; Bracco Imaging, Milan, Italy) was added to the BMF as tracer particle after which images were captured for 3 s using a three-angled plane wave acquisition scheme (4 MHz center frequency, single-cycle pulse) at 9000 fps.These settings are feasible in the clinical setting as well since the acoustical output is within safety limits.The 3-s acquisition was performed four times to capture a minimum of 10 cardiac cycles.Because of the limited aperture of the transducer (38 mm), two consecutive field of views were acquired to capture the entire region of interest.First, the CFA and proximal bifurcation were imaged.Second, the distal bifurcation, SFA, and DFA were imaged.The fields of view had an overlap of ∼1 cm, which was later used to automatically register the field of views and fuse the velocity vector fields.
The reconstructed images for the different angles were coherently compounded-yielding an effective frame rate of 3000 fps-and filtered by means of singular value decomposition.The lower and higher rank thresholds were semi-automatically found based on the spatial similarity matrix. 28Similar to oPIV, the 2D velocity vector fields were calculated through PIV analysis.
Here, a modified implementation of PIVlab was employed, as previously used for ePIV in a clinical setting. 16,29For the control model, six iterations with square interrogation windows of 32, 32, 16, 16, 8, and 8 pixels and 75% overlap were used, corresponding to a final isotropic vector grid spacing of 0.22 mm.For the stenosed model, two initial iterations with a square size of 64 pixels were added to capture the higher stenotic velocities.Correlation averaging of 10 frames was applied to increase the signal-to-noise ratio, resulting in 300 vector fields per second.Final postprocessing consisted of a spatial Gaussian filter (3 × 3 velocity vectors) and a threeensemble temporal moving average filter.

Doppler Ultrasound
Pulsed-wave Doppler measurements were performed at the CFA, SFA, and DFA using a conventional clinical US system with a linear 14L5 transducer (ACUSON S2000; Siemens AG, München, Germany).Clinical guidelines for vascular laboratory examinations 30 were followed: velocities were obtained at the center of the lumen with a sampling volume of 1 mm and a pulse repetition frequency of 13 kHz.The beam-to-velocity angle was kept below 60 deg.

Computational Fluid Dynamics
The SimVascular finite-element solver (Stanford & UC Berkeley, Palo Alto & San Francisco, California) was used for CFD simulation of the two geometries.This open-source solver has been optimized for cardiovascular flows with pre-conditioners accounting for resistive or Windkessel boundary conditions and stabilization of reversed flow in the outlets. 31The geometries were segmented from the CT-scans of the PDMS models, using a regularized level-set algorithm. 32Tetrahedral elements were used to construct a volumetric mesh with a three-element prism boundary layer to capture the high shear gradients near the wall.The vessel wall was considered rigid and the density and dynamic viscosity of BMF were used as material properties.
For both geometries, a mesh convergence study was performed for peak systolic flow rate.Convergence was defined as an overall smaller than 3%-change for the centerline velocities in the CFA and SFA for a doubling of the mesh points.Mesh convergence for the stenosed geometry, in which transitional flow was present (characterized by significant turbulent kinetic energy), was achieved for the time-averaged velocity along the centerline (Appendix C).
For the two geometries, the flow boundary conditions were derived from both DUS and ePIV, as these techniques are suitable for obtaining flow boundary conditions in patients.For DUS, the centerline velocity was converted to a flow rate over time, by using Womersley theory under the assumption of a fully developed flow profile and the vessel diameters measured on CT. 33 For ePIV, which provides a 2D spatial velocity profile along the diameter in the CFA, SFA, and DFA, the flow rate was computed as follows: where Q is the flow rate, R is the radius of the cross-section, and u is the axial velocity.With the velocity only known for polar angle θ ¼ 0 and θ ¼ π, this equation assumes that the angleaveraged velocity for a radial coordinate is equal to the average of the values measured at θ ¼ 0 and θ ¼ π.

Analysis of oPIV, ePIV, and CFD Data
For all modalities, 10 cardiac cycles were aligned at peak systole and ensemble-averaged over time to compare the characteristics of the mean transitional flow.The cycle-averaged data were used to extract various flow parameters including spatial and temporal velocity profiles, VC (see below), and WSS.The spatial and temporal velocity profiles were obtained at the CFA, SFA, and DFA of the control model, and at mid stenosis and at the post-stenotic jet of the stenosed model.For all temporal velocity profiles, the absolute peak differences between the three modalities relative to oPIV were calculated.
To statistically assess the agreement of the three techniques relative to oPIV, velocity vector data from ePIV and CFD were resampled using linear interpolation (i) in space at the vector locations of the oPIV grid and (ii) in time to 100 sample points of the averaged cardiac cycle.The mean and standard deviation (SD) of the difference in velocity magnitudes for all spatiotemporal velocity vectors, referenced to oPIV values, were calculated to assess the agreement in velocity in space and time between the modalities.This was repeated for all velocity vectors during peak systole to assess the agreement in space between the modalities at this single timepoint.The limits of agreement were defined as mean AE1.96 SD, covering 95% of the area under a normal curve. 34wo flow-related parameters were derived from the velocity vector fields and compared between the modalities.First, VC was obtained as a metric to describe the distribution in vector direction in which completely unidirectional parabolic flow will have a value of 0 and recirculating complex flow with a large spread in vector direction will have a value toward 1. 35,36 VC was computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 1 3 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 7 ; 7 3 6 x ¼ 1 n E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 7 ; 6 8 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 7 ; 6 5 4 where θ i is the beam-to-flow angle of the vector v x;i ; v y;i at position i, x and y are the averaged beam-to-flow angles represented on the unit circle, and n is the total number of vectors in the region of interest.VC was derived in a region around the stenosis, i.e., 1 cm upstream until 1 cm downstream of the stenosis.The same area was selected in the control model.Second, timeaveraged WSS (TAWSS) was determined for the inner and outer walls of the CFA and SFA for the oPIV and ePIV data by interpolating the wall-parallel velocity component along a line perpendicular to the wall.Subsequently, a weighted cubic spline was fitted to the velocity component parallel to the wall and used to calculate TAWSS. 37For both CFD simulations, TAWSS was provided by the solver based on the shear rate tensor.

Control Model
Qualitatively, similar flow features were observed in the control model for ePIV, CFD-DUS, and CFD-ePIV compared to the oPIV reference [Fig.2(a) and Video 1].For all approaches, a flow separation at the outer SFA and DFA walls near the bifurcation was observed during forward flow with high velocities at the inner walls of the SFA and DFA.During backflow, the oblique flow originating from the DFA interacted with the SFA flow and created undulating flow structures.These flow structures were slightly less defined for ePIV, yet exaggerated and with a phase-offset for CFD-DUS.Temporal and spatial velocity profiles of the CFA showed moderate to good agreement with an absolute peak velocity difference of 0.19 cm∕s between oPIV and ePIV and 21.6 cm∕s between oPIV and CFD-ePIV [Figs.2(b) and Appendix D].CFD-DUS showed higher velocities with an absolute peak difference of 53.7 cm∕s compared to oPIV and a larger diastolic phase offset.Near the wall, both CFD approaches followed the oPIV curve with a sharp transition toward lower velocities of ∼10 cm∕s, whereas the ePIV curve appeared smoother and velocities remained relatively high around 50 cm∕s.For the SFA, absolute peak velocity differences were 2.9, 46.0, and 5.0 cm∕s for ePIV, CFD-DUS, and CFD-ePIV, respectively, with respect to oPIV.The corresponding spatial velocity profiles presented higher velocities at the inner wall of the SFA due to the bifurcation.For the DFA, larger deviations were found between oPIV and ePIV, and between oPIV and CFD-ePIV, with absolute peak differences of 11.0 and 30.7 cm∕s, respectively.The difference between oPIV and CFD-DUS was slightly smaller compared to the other locations with a peak difference of 34.0 cm∕s.All modalities showed a similar trend in the spatial velocity profiles with a steep gradient at the inner wall of the DFA.Mean differences of the entire CFA-SFA region averaged over all time points, referenced to oPIV, are minimal for ePIV and CFD-ePIV with similar limits of agreement (Table 1).The spatial plot of the SD [Fig.4(a)] demonstrates the highest spread for ePIV and CFD-ePIV relative to oPIV in areas with high velocities during backflow.CFD-DUS showed an overall larger disagreement with oPIV.showed better agreement with oPIV, with an absolute peak difference of 23.2 cm∕s.During diastole, CFD-DUS presented a phase off-set, whereas the other approaches were well-aligned.No large cycle-to-cycle variations were found at the centerline velocity of the jet, whereas closer to the wall, at locations with strong spatial gradients, this variation increased.The spatial profile captured by ePIV was more parabolic compared to the flattened profile of the other techniques but estimated similar low near-wall velocities at the level of the jet.For the stenosis, however, ePIV measured increased near-wall velocities with a larger cycle-to-cycle variation.Overall, mean differences showed minimal bias in ePIV and CFD-ePIV velocity magnitudes, with a moderate positive bias for CFD-DUS [Figs.4(b) and Table 1].For ePIV, the largest differences and spread were found specifically at the inlet of the stenosis, whereas the entire stenosed region showed small positive biases for CFD-ePIV.The largest differences for CFD-DUS were found    in a similar region inside and around the stenosis.Here, the variation over time was elevated as well with increased limits of agreement compared to the control model and compared to the other modalities.

Vector-Derived Parameters
VC showed good agreement between all techniques for both models, except for the diastolic phase offset and corresponding time mismatch in peaks during flow reversal (Fig. 5).In the control model, mean peak systolic VC was <0.01 for all techniques.In the stenosed model, the mean peak systolic VC was elevated compared to the control model and equaled 0.17 AE 0.04, 0.16 AE 0.02, 0.13 AE 0.04, and 0.12 AE 0.02 for oPIV, ePIV, CFD-DUS, and CFD-ePIV, respectively.When only the post-stenotic recirculation was analyzed, excluding the jet, VC increased to 0.69 AE 0.20, 0.75 AE 0.15, 0.34 AE 0.13, and 0.41 AE 0.15.
Overall higher TAWSS values were found with ePIV, CFD-DUS, and CFD-ePIV compared to oPIV in the control model (Fig. 6).Lowest WSS was located at the CFA wall where the DFA splits off, whereas the highest WSS was found at the inner SFA wall.Mean TAWSS for the segment given in Fig. 5(e) was 1.43, 2.50, and 1.80 Pa for ePIV, CFD-DUS, and CFD-ePIV, respectively, compared to 0.97 Pa for oPIV.For the stenosed model, highest TAWSS values were found within the stenosis for all techniques.At this location (around z ¼ −1 cm), both CFD approaches measured two-to threefold higher TAWSS values compared to oPIV.For ePIV, similar values to oPIV were found for the left CFA vessel wall, yet a twofold increase in TAWSS was present at the right CFA vessel wall.This asymmetry was not seen in the other techniques.

Discussion
The accuracy of blood flow assessment was investigated in models of the femoral bifurcation by clinically applicable ePIV measurements and CFD simulations.Using oPIV measurements as the gold standard, ePIV and CFD-ePIV demonstrated good overall agreement in vector velocities in a stenotic flow with limits of agreement of 0.1 to 1.5 AE 34 cm∕s.The more clinically available CFD-DUS approach was less accurate with limits of agreement of 16 AE 55 cm∕s in the stenosed model.For the vector-derived parameters, good agreement was found for VC between all four techniques in regions with laminar flows as well as regions with complex flow phenomena, such as a recirculation.TAWSS showed reasonable agreement for all modalities in the control model.Larger deviations of TAWSS were found in the stenosed area, whereas good agreement was found outside of this region for ePIV.The accuracy of the oPIV-reference for TAWSS assessment was debatable, however.
All three approaches demonstrated largely similar flow dynamics in the control and stenosed model.The complex transitional dynamics of the post-stenotic jet were well captured by the ePIV and CFD techniques, including its breakdown, recirculating flows near the jet and high spatial velocity gradients at the border of the jet.This generates confidence for using these techniques to assess complex flow velocities in clinical practice in general and more specifically in patients with PAD.
Several disagreements with oPIV were present.For ePIV, spatial velocity profiles in the post-stenotic jet and in near-wall regions had a smaller gradient compared to oPIV.In addition, velocities measured within the stenosis were higher than obtained with oPIV.The differences in the spatial velocity profiles can be attributed to smoothing by intrinsic factors that include spatial averaging by the UCA point-spread function, which is larger than the speckle size, and by the PIV correlation kernel dimensions. 38The latter could be improved by using anisotropic kernels aligned with the spatial velocity gradients.The Gaussian spatial averaging in the post-processing step likely also contributed to these differences, but this step is required to suppress noise and is therefore part of the intrinsic limitation of the technique, so far. 16The increased velocities within the stenosis for ePIV possibly originated from a difference in compliance between the optical and acoustical models.The thin-walled ePIV resin model exhibited less pulsatile wall motion (0.5% to 1% as observed by pixel displacement) than the oPIV PDMS block model (2% to 4% as observed by pixel displacement), which led to a larger stenotic diameter in the oPIV phantom during peak systole relative to the ePIV phantom.Under the same peak flow rate, this would be associated with higher velocities in the ePIV phantom.The difference in phantom compliance For CFD-ePIV, steep gradients in velocity profiles were sufficiently sampled, but still similar limits of agreement as ePIV were found in the order of AE25 cm∕s (control) and AE34 cm∕s (stenosed) relative to oPIV.Judging from the evolution of the temporal velocity profiles [e.g., Figs.2(b), 2(c) and Figs.3(b), 2(c)], an important contributor to this spread was the slight offset in timing of the backflow and subsequent forward diastolic flow.Here, the difference in compliant behavior was likely responsible too.Second, peak forward flow was overestimated in the control model by CFD-ePIV, possibly due to a positive bias of the flow rate estimation by ePIV, as higher near-wall velocities were registered by ePIV.This has an important contribution to the overall flow rate estimation as velocities are multiplied by the radius (Eq.1).For CFD-DUS, an even larger positive bias in the velocity magnitude was present, as well as important differences in the phase and velocity profile of backflow and diastolic flow.These differences can be explained by spectral broadening 39 and finite transit time. 6Clinical DUS examinations that we followed, especially for the small sampling volume of 1 mm, have been associated with peak velocity overestimation of 10% to 50%, depending on the device used. 6The variability in errors between devices, as well as the high interobserver variability of 20% for the SFA, 8,9 makes it hard to systematically correct for these errors in clinical scenarios.The phase difference in backflow resulted from deviations from a (theoretical) Womersley profile during this phase, due to reversed flow through the bifurcation creating strongly asymmetric spatial flow profiles.
VC was very high during flow reversal phases, which complicates interpretation of this parameter for vessels that exhibit backflow.Instead, a focus on peak systole could be relevant to quantify areas with disturbed flow, as VC was consistently higher at this time point for the stenosed area relative to the control model.TAWSS demonstrated similar profiles along the CFA and SFA wall, with an overall increase in magnitude for all three modalities compared to oPIV in the control model and inside the stenosis.This positive bias in TAWSS for ePIV might in part be caused by the higher near-wall velocities in ePIV, as discussed previously.Another potential factor is the finer vector grid spacing of ePIV, which could also explain the asymmetry in TAWSS values for the stenotic model in ePIV.A smaller distance between the nearest wall velocity for the right wall relative to the left wall could have led to increased TAWSS at the right wall.The CFD approaches showed overall increased TAWSS values.For the PIV techniques, a suboptimal resolution for the near-wall velocities 40 and the chosen curve fitting method are known to impact TAWSS estimates. 41These factors limited oPIV assessment of TAWSS, making it hard to judge to what extent ePIV and CFD approaches were inaccurate and to what extent the oPIV reference was truthful.This could be improved by using smaller camera windows and using interfacial PIV, which provides more accurate, higher TAWSS values relative to the standard oPIV applied in this study. 19hile ePIV was able to measure TAWSS in the physiological order of magnitude, translation into the clinic remains challenging as a precise detection of the vessel wall and corresponding segmentation of the flow lumen is difficult.Furthermore, current vector resolution limits the accurate assessment of near-wall velocities.For CFD approaches, the clinical application for estimating TAWSS can be more straightforward.
This study extends previous experiments that have validated ePIV and finite-element CFD in various geometries and modeling strategies.These include ventricular, 14 aortic, 42,43 coronary, 23 and carotid 15,22 artery phantoms.For ePIV, the differences in velocity magnitude of about 9% to 17% have been reported, 14,15 which was of similar order in our study (4% to 11%).For CFD, we used flow boundary conditions from ePIV or DUS, instead of the benchmark technique (oPIV or PC-MRI).Here, we wanted to study the technical validity of the CFD methodology, where our approach relies on ultrasound-based methods that can be applied with low cost (DUS) or with more advanced ultrasound imaging (ePIV).Furthermore, our study supports the validity of an under-resolved direct numerical simulation strategy for estimating cycle-averaged stenotic flows, as previously reported for a carotid artery stenosis. 22uture validation studies can focus on a comparison of the ePIV and CFD techniques to other clinically available vector flow imaging modalities.These include PC-MRI, ultrasound blood speckle tracking techniques, vector Doppler techniques, and 3D-ultrasound. 13,44Depending on patient anatomy, flow conditions, and clinical purpose, all techniques have different strengths and weaknesses, which will require different and possibly synergistic approaches where necessary.PC-MRI can provide time-resolved 3D velocity fields and is less susceptible to calcifications but is relatively limited in spatial and temporal resolution, prone to stent artifacts, and requires interleaved sampling, which filters details on cycle-to-cycle variations of complex flow behavior and may be challenged by variable heart rates. 45Native blood speckle tracking techniques are cost-effective, but may be less attractive due to limited image contrast for deeper-lying vessels. 46CFD approaches can provide high-resolution 3D velocity fields and WSS at low additional cost, but its accuracy will depend on input and output flow boundary data and the accuracy of CFD computed WSS is uncertain.
The in vitro validation of vector flow imaging modalities is an important step toward its potential implementation in clinical practice.The 2D or 3D visualization and quantification of blood flow characteristics could provide in-depth insights in local hemodynamics superior to DUS examination as used in current clinical decision making.In addition, vector derived parameters, such as VC and WSS, could be used to quantify complex flow phenomena and identify locations potentially at risk for plaque formation, in-stent restenosis 47 or stentgraft thrombosis. 48An extra potential benefit of CFD models is the ability to predict the effect of different treatment strategies on local hemodynamics to select a treatment that creates an optimal hemodynamic environment, with the aim to reduce reintervention rates. 49,50The feasibility of ePIV and CFD in patients with atherosclerosis has already been shown in small patient populations. 16,47,48,51,52However, longitudinal studies with larger cohorts are necessary to show the potential benefit of these novel modalities in clinical practice and to study the relation between vector derived parameters and clinical outcome.
This study has several limitations.The most important limitation was the use of different models for the oPIV and ePIV measurements due to incompatible transmission constraints.Care was taken to minimize differences in anatomic parameters (AE0.25 mm) and flow rate (AE0.1 mL∕s), and furthermore a difference in compliant behavior (2% to 3% more radial distension) was present.It can be estimated that in the CFA, a 5.6% difference in radius would lead to an 11.5% decrease in velocity in the oPIV phantom relative to ePIV, with an additional 6% decrease during systole due to the increased pulsatility of the oPIV model.In the stenosed CFA, the same radial difference could cause a 23% difference in velocity and could explain the elevated velocities in the ePIV model, with similar implications for the TAWSS.Furthermore, although the US and laser plane were carefully positioned, slight misalignment could have been present and have influenced the outcomes.Also, the PIV analysis differed slightly between both modalities.Equal kernel sizes were used for oPIV and ePIV, yet ePIV had a finer grid because of the 75% overlap (oPIV 50% overlap).A 75% overlap was not available in the PIVlab toolbox that was used for the oPIV analysis.Differences in grid spacing were offset in the statistical assessment between the modalities, as all grids were interpolated on the oPIV grid.Furthermore, the effects of non-arterial tissue and imaging artefacts by calcifications, stents or other pathologies were not investigated, which can constitute major limitations for ePIV imaging in patients 16 or for CTA segmentation 53 as input for the CFD approaches.Last, although all ultrasound hardware, pulse characteristics, and post-processing steps were set according to clinical use, tissue attenuation was not taken into account.

Conclusion
This study has shown that ePIV and CFD-ePIV can accurately estimate complex stenotic transitional flows.CFD-DUS mostly overestimated the velocity magnitudes, which can be attributed to DUS measurement inaccuracies and the limited information that it provides for a flow rate reconstruction.Nevertheless, similar post-stenotic and bifurcation related flow patterns were observed between all modalities with comparable VC values and trends in TAWSS.Although the applied ePIV and CFD techniques require further development for a routine application in clinical care, they are attractive from a cost-effectiveness perspective and provide the ability to assess hemodynamic insights with high temporal and spatial resolution.
6 Appendix A: Averaged Anatomy and Flow Curve

Anatomy
The diameters of the three femoral arteries were based on a literature study, yielding the characteristics as shown in Table 2.The weighted mean and weighted standard deviation of the reported van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .cohorts were computed.The mean diameters of the three arteries were used to design the femoral bifurcation phantom models.For the SFA, a diameter of 6.2 mm instead of 6.3 mm was used.

Flow Curve
Pulsed-wave Doppler ultrasound (DUS) measurements were performed at the CFA, SFA, and DFA in seven healthy participants using a conventional clinical US system with a linear 9L3 transducer (iU22; Philips, Amsterdam, The Netherlands).Clinical guidelines for vascular laboratory examinations were followed: velocities were obtained at the center of the lumen with a sampling volume of 1.5 mm and a pulse repetition frequency of 2.1 to 4.2 kHz.The beam-tovelocity angle was kept below 60 deg.
The measurements of the DFA were of mediocre quality (broad spectrum with unclear peak velocity), so only the measurements of the CFA and SFA were used.Under the assumption of fully developed Womersley flow, the flow rate was computed similar as described for the DUS measurements in the paper methods using Eq. ( 1).Because of significant differences in interval between peak forward and backward flow, a median value of this interval in the seven participants was computed.The timescale of the flow rate curves was rescaled so that all had an interval equal to this median value.Subsequently, the ensemble average curve was computed.
The CFA curve served as reference for the inflow boundary condition (Fig. 7).The DFA flow rate was considered equal to the CFA minus SFA flow curve.Together with the SFA curve (Fig. 7), they were used to compute the mean flow split and the capacitance split for the outflow Windkessels.

Appendix B: Phantom Geometries
From the cone-beam CT scan, the lumen of the oPIV en ePIV models was segmented with the VMTK level-set algorithm 32 based on the contrast between the wall and the lumen filled with air.
With the inscribed-sphere VMTK algorithm, the radius along the centerline of the models was assessed.The radius along this centerline for the models is shown in Fig. 8.

Appendix C: CFD Mesh Convergence for Stenosed Bifurcation
For the stenosed model, the convergence of a steady inflow, equal to the peak systolic flow rate (54.2 mL∕s), was assessed for three meshes.The mesh convergence for velocity, pressure, VC, and wall shear stress was assessed.Mesh size was set such that mesh 2 included roughly twice as much points as mesh 1, and vice versa for mesh 3 relative to mesh 2. For all meshes, a regional mesh refinement (twice as many points per volume) was applied within a spherical volume that included the stenosis and its downstream region extending up to the bifurcation.The approximate mesh characteristics were as follows.For the applied steady inflow rate, velocities in the stenotic jet of about 310 cm∕s were expected.This is associated with a Reynolds number of 3780 in the region downstream of the stenosis: where v is the spatial mean velocity, ρ is the fluid density, L is the local diameter, and μ the viscosity, all in mm, kg, and s.Under Kolmogorov homogeneous isotropic turbulence theory, the associated length scale of the smallest flow structures would be E Q -T A R G E T ; t e m p : i n t r a l n k -; s e c 8 ; 1 1 4 ; 6 6 8 η ¼ Thus, for mesh 3, the edge length in the area of transitional flow is one order of magnitude higher than the estimated Kolmogorov length scale under isotropic turbulence.
For the points along the common femoral-SFA centerline (see Appendix B), the temporal transitional flow statistics were assessed by computing the temporally averaged solution and the variance for the pressure and velocity magnitude variables (Fig. 9).For the mean velocity and pressure field, mesh 2 was considered converged relative to mesh 3 (mean absolute difference along centerline in pressure of 0.2% and in velocity of 2.2%, with a maximum point velocity difference in post-stenotic area of 8.5%).Nonetheless mesh 3 was used for further (pulsatile) simulations.
From the pressure variance and turbulent kinetic energy (Fig. 9), it can be observed that the mesh was not converged for the fluctuating components of the turbulence, with the turbulent kinetic energy decreasing for more refined meshes, where a mean spatial absolute difference along the centerline in turbulent kinetic energy of 21.2% was present between mesh 2 and mesh 3.
For VC, the same region of interest as in the main method that included the stenosis and the post-stenotic jet was used for analysis.For the three meshes, the value of VC was 0.121, 0.117, and 0.117 for mesh 1, mesh 2, and mesh 3, respectively.
For wall shear stress, the left wall that spanned the CFA, the stenosis, and the superficial artery was used for analysis.The average wall shear stress along this wall was 8.68, 9.20, and 9.87 Pa for mesh 1, mesh 2, and mesh 3, respectively.Ashkan Ghanbarzadeh-Dagheyan received his PhD in 2018, focusing on methods for breast cancer detection using microwave and ultrasound imaging, from Northeastern University, USA.His is currently a postdoctoral researcher at the University of Twente.His main focus is the application of biomedical imaging in disease detection and characterization, especially in cancer and cardiovascular maladies.He currently works on improving methods for ultrasound velocimetry and applying photoacoustic imaging in cancer detection on a cellular and tissue level.

Fig. 1
Fig. 1 Experimental setup.(a) Control and (b) stenosed geometry used for the flow models.(c) Circulating flow loop.R, resistance.
van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .

3. 2
Stenosed Model Similar flow patterns were observed in all techniques in the stenosed model [Fig.3(a) and Video 2].All approaches showed the post-stenotic jet extending up to the bifurcation point during forward flow.High spatial velocity gradients were captured at the jet border, especially at the left side where strong reversed flow velocities were present [inset Fig. 3(a)].Post-stenotic recirculation zones were captured by all modalities.During flow reversal, all techniques showed similar dynamics including a reversed jet into the CFA with associated recirculation zones.Higher stenotic velocities, with an absolute peak difference of 57.7 and 57.0 cm∕s, respectively, were captured by ePIV and CFD-DUS compared to oPIV [Fig.3(b) and Appendix D].CFD-ePIV van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .

Fig. 2
Fig. 2 (a) In-plane velocity vector fields at peak systole in the control model.Black line and square indicate the locations where (b) velocity profiles were extracted.Cycle-averaged temporal (upper row) and spatial (lower row) velocity profiles extracted at (from left to right) the CFA, SFA, and DFA.Shading represents the standard deviation over 10 cardiac cycles.Video 1 shows animations of (a) over time and Appendix D provides spatially averaged velocity profiles of (b).

Fig. 3
Fig. 3 (a) In-plane velocity vector fields at peak systole in the stenosed model.(b) Similar spatiotemporal overview as Fig. 2, but for a location in the stenosis and downstream in the post-stenotic jet, as indicated by the black markers.Video 2 shows an animation of (a) over time and Appendix D provides spatially averaged velocity profiles of (b).
van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .

Fig. 4
Fig. 4 Mean differences (top) and standard deviation (bottom) in velocity magnitudes over the entire cardiac cycle between ePIV, CFD-DUS, and CFD-ePIV relative to oPIV in (a) the control model and (b) the stenosed model.

Fig. 5
Fig. 5 VC over time extracted in the CFA of the control (top) and stenosed (bottom) model.Periods of flow reversal are indicated in grey.Shading represents the standard deviation over ten cardiac cycles.

Fig. 6
Fig.6In-plane TAWSS at the CFA and SFA walls for the control (top) and stenosed (bottom) model obtained with all modalities.A comparison of the CFA TAWSS between z ¼ −2 and z ¼ 1.5 cm between all modalities is given on the right-hand side.

Fig. 7
Fig. 7 Flow rate of the (a) CFA and (b) SFA.

Fig. 8
Fig. 8 Segmented stenosed geometry with the computed centerline for the CFA and SFA.The vessel radius of the two models (oPIV; PDMS and ePIV; resin) along this centerline for control and stenosed model are shown in the plots.

Fig. 9
Fig. 9 Mesh convergence study for points along the CFA and SFA centerline in the stenosed model.

Fig. 10
Fig. 10 Temporal (left) and spatial (right) velocity profiles of the CFA, SFA, and DFA in the control model.Shading represents the standard deviation in space probed over five points at a 0.5 mm distance of each other.

Fig. 11 3 )
Fig. 11 Temporal (left) and spatial (right) velocity profiles at mid stenosis and the post-stenotic jet in the stenosed model.Shading represents the standard deviation in space probed over five at a 0.5 mm distance of each other.

Hadi
Mirgolbabaee received his BSc degree in mechanical engineering in 2017 from Babol Noshirvani University of Technology, Babol, Iran.He received his MSc degree in mechanical engineering in 2019 from the University of Twente, Enschede, The Netherlands.He is currently a PhD student at multi-modality medical imaging (M3i) and physics of fluid (PoF) groups at van de Velde et al.: Validation of ultrasound velocimetry and computational fluid. . .

Table 1
Limits of agreement (mean ±1.96 SD) for all vector points in the CFA and SFA between the three modalities relative to oPIV.
The mean magnitude for all velocity vectors in the oPIV model was 25.6 and 35.6 cm∕s for the control and stenosed model, respectively.

Table 2
Literature review on femoral artery diameters.
a Based on reported normal range.IVUS, intravascular ultrasound.Note: Values in Italics show the weighted mean and weighted standard deviation.