An LDV based method to quantify the error of PC-MRI derived Wall Shear Stress measurement

Wall Shear Stress (WSS) has been demonstrated to be a biomarker of the development of atherosclerosis. In vivo assessment of WSS is still challenging, but 4D Flow MRI represents a promising tool to provide 3D velocity data from which WSS can be calculated. In this study, a system based on Laser Doppler Velocimetry (LDV) was developed to validate new improvements of 4D Flow MRI acquisitions and derived WSS computing. A hydraulic circuit was manufactured to allow both 4D Flow MRI and LDV velocity measurements. WSS profiles were calculated with one 2D and one 3D method. Results indicated an excellent agreement between MRI and LDV velocity data, and thus the set-up enabled the evaluation of the improved performances of 3D with respect to the 2D-WSS computation method. To provide a concrete example of the efficacy of this method, the influence of the spatial resolution of MRI data on derived 3D-WSS profiles was investigated. This investigation showed that, with acquisition times compatible with standard clinical conditions, a refined MRI resolution does not improve WSS assessment, if the impact of noise is unreduced. This study represents a reliable basis to validate with LDV WSS calculation methods based on 4D Flow MRI.


Scientific Reports
| (2021) 11:4112 | https://doi.org/10.1038/s41598-021-83633-y www.nature.com/scientificreports/ (FEM) starting from a cubic interpolation of the velocity field on the nodes of a 3D mesh based on MRI voxels. Finally, Piatti et al. 13 developed an innovative method to compute WSS based on 2D and 3D Sobel filters. PC MRI-derived WSS calculation methods are often validated with Computational Fluid Dynamics (CFD) techniques in realistic vascular geometries. Nevertheless, CFD employment is limited in complex physiological conditions 14 . Alternatively, in Laser Doppler Velocimetry (LDV) campaigns, it is possible to reproduce the same experimental conditions of in vitro PC MRI acquisitions in terms of geometry, flow rate, pressure, fluid viscosity, and density, providing accurate velocity profiles to serve as a reference for PC MRI data. LDV enables an accurate measurement of absolute velocity components on a single point with the drawback of having to move the measurement volume to obtain velocity mapping, and it represents an intrinsic, robust velocimetry technique since velocity data are obtained through robust statistical post-processing tools. In the scientific literature, few works about LDV validation for PC MRI-derived WSS computation are reported. Bauer et al. 15 successfully investigated WSS accuracy from PC MRI data with LDV in a simplified model of an aortic aneurysm, employing in parallel CFD for direct comparison.
The present work aims to develop a robust method based on a hydraulic circuit, in which LDV served as a standard for velocity measurements, to be employed as a reference to assess the accuracy of 2D PC and 4D Flow MRI velocity acquisitions. This method is intended to be employed to evaluate any potential developments in acquisition and flow quantification protocols, particularly for the calculation of WSS. Firstly, a hydraulic test rig was designed and manufactured to be MRI-compatible and to allow accurate LDV measurements under both steady and pulsatile flow conditions. Next, velocity data obtained with the test rig were statistically analyzed to assess the accuracy and reproducibility of results between the two velocimetry techniques. The correspondence between LDV and PC MRI WSS profiles was statistically assessed. Finally, the efficacy of the method developed in this work was tested evaluating the influence of 4D Flow MRI spatial resolution on WSS profiles.

Methods
Experimental set-up. The experimental set-up ( Fig. 1) was designed to keep limited dimensions, to be easy to move between facilities, and to be compatible with LDV and MRI constraints. It was based on an MRIcompatible gear pump (CardioFlow 5000, Shelley Medical Technology), a test section and a reservoir open to the www.nature.com/scientificreports/ atmospheric pressure (Fig. 1a), completed with several connecting flexible PVC pipes for ease of assembly and transport between LDV and MRI facilities. The pump provided a TTL signal to synchronize both LDV and MRI acquisitions, the period of flow rate signals was T c = 1 s for both steady and pulse acquisitions. Made of a mixture of 40% by weight of glycerol in water, the working fluid reproduced blood density ( ρ = 1098 kg/m 3 ) and viscosity ( µ = 0.00345 Pa·s ) closely 16 . The test section was a compact 25 × 25 × 345 mm squared PMMA pipe (Fig. 1b).
The measurement area was located on a cross-section placed 265 mm downstream of the test section inlet, to limit the influence of the test section inlet and outlet. Investigations were performed under laminar steady and pulse flow conditions, deploying an ultrasound flowmeter (Flowmax 44i, MIB Gmbh) to ensure the flow rate control in the measurement location, Q exp ( Q exp = 98.67 ml/s and Q exp = 16.47 + 35.64 sin (2πt + 0.50) ml/s for steady and pulse flow rate respectively). During MRI measurements, a plastic box filled with agarose gel (A8963, Applichem) surrounded the test section to reproduce human body tissues and to reduce noises (Fig. 1c).
Phase contrast MRI. Firstly, one standard 2D PC and one standard 4D Flow MRI sequence were used on a 1.5 T MRI system (Aera 1.5T, Siemens), their parameters were selected as in current clinical practice, defining then tests called Steady2D X0 , Pulse2D X0 , Steady4D X0 , and Pulse4D X0 (Table 1) 17 . In view to fit the maximum expected velocity v max , the Velocity ENCoding (VENC) was adjusted as in Eq. (1) ( v max = 28 cm/s and v max = 12 cm/s for steady and pulse tests respectively).
Next, a second experimental campaign was performed employing the 4D Flow sequence. The size of the voxels was progressively reduced from 2.22 × 2.22 × 2 to 1.86 × 1.86 × 1.86 , and then to 1.50 × 1.50 × 1.50 mm. The refinement was carried out with isometric voxels to simplify the analysis; the starting size was considered as "almost" isometric since the smallest dimension was shorter than the others by about 10%. Consequently, the Field of View (FOV) was proportionally scaled. Acquisition parameters are resumed in Table 1, defining instead tests subsequently called Steady4D X1 , Steady4D X2 , Pulse4D X1 , and Pulse4D X2 . Markers on the experimental setup enabled the placement of the measurement location at the center of the MRI scanner, drastically reducing the effects of the magnetic field inhomogeneities. The velocity signal from the static homogeneous region composed by the agarose gel was exploited to define the impact of the velocity offset from magnetic field inhomogeneities and noise. The distribution of this signal was a Gaussian with a non-zero mean and variance, representing the phase offset and the velocity noise respectively. As a first approximation, the influence of the offset was considered negligible; its ratio with the mean of the velocity signal ( v fluid ) in the fluid domain was lower than 3%. To evaluate the impact of noise on PC MRI acquisitions, Velocity to Noise Ratio (VNR) was defined as in Eq.
(2) ( σ v gel is the variance of velocity noise in the agarose gel calculated as the standard deviation of the velocity signal of pixels in the agarose gel).
(1) VENC ≈ v max + 10%v max for steady tests v max + 30%v max for pulse tests www.nature.com/scientificreports/ Laser doppler velocimetry. For LDV investigations, a low concentration ( C volume ≈ 0.05 %) of 10 µ m silver-coated hollow glass spheres (S-HGS, Dantec) was added to the mixture. They exhibited a high refractive index, increasing the LDV SNR with respect to non-coated hollow glass sphere seeding. The relative density gap with the fluid ( ρ fluid /ρ particle ≈ 0.70 ) did not induce undesirable behaviors like a lift or gravitational effects, particle accumulation, or modifications of the working fluid physical properties. The axial velocity along the x-axis (see Fig. 1b for axes reference) was measured by an LDV system (fp50 LDV system). The optical path of the laser beams passing through the upper face of the test section and the fluid was defined according to Snell's refraction law (Fig. 1b) (Eq. 3, α represents the incidence, in the air, or refraction, in the fluid, angle, and n the refractive index).
The location of a measurement point was calculated considering the refractive index mismatch at the optical interface along the z-direction, discounting the window thickness. To optimize the whole acquisition time, measurements were obtained on a square grid of 63 × 63 points equally spaced by 0.40 mm; this length was the LDV data spatial resolution. For each measurement point, the velocity was recorded for 30 s and velocity data were obtained with a home-made Matlab script (Matlab R2016b, MathWorks). For steady tests, a two-class Gaussian Mixture Model was applied to the velocity sample distribution to separate the signal Gaussian distribution, determining the mean velocity and measurement uncertainty σ s , from velocity outliers, determined during zero flow rate measurements. For pulse acquisitions, the velocity signal was sampled in periods as long as those of the pump flow rate; all periods were overlapped and finally fitted with a sinusoid, according to the Womersley theory of velocity profiles under pulse flow conditions, for which velocity at each point has the same waveform of the inlet flow rate 18 . The coefficient of determination R 2 (best fit for R 2 = 1) was employed to assess the uncertainty of the measurement. A CFD verification study was performed to preliminary control LDV measured velocity profiles on a simple test case for which CFD is deemed reliable. Numerical simulations were performed using the Finite Volume Method (Star CCM+, Siemens). Calculations were based on a Cartesian mesh of 1725000 0.5mm-size cubic elements and results were then linearly interpolated to fit the LDV grid (Matlab R2016b, MathWorks) for direct comparison. For the Steady test, the mean value of σ s was 6.63 × 10 −4 m/s, while for the Pulse test the mean value of R 2 was 0.82. The Locally Weighted Scatterplot Smoothing (LOWESS) method (Matlab R2016b, MathWorks) was applied to experimental data to find velocity profile boundaries, smooth noises, and scale results in a uniform grid. Ultimately, a broader LDV grid was generated to approach a similar resolution of MRI data to allow direct comparison, with the average of all LDV points contained in a window with the identical size of one MRI pixel.
Wall shear stress calculation. WSS is defined as (4) where µ is the blood dynamic viscosity, v the blood velocity, and n the normal direction to the arterial wall. The methods proposed by Stalder et al. 10 (2D-WSS) and Potters et al. 11 (3D-WSS) were employed to assess WSS from PC MRI and LDV data-sets in the same location of the test section, in view to demonstrate the efficacy of this experimental system to validate PC MRI post-processing algorithms for WSS computation. For 2D-WSS, from 2D PC or sliced 4D Flow MRI data, the lumen is contoured through B-splines and within this fluid domain, the blood velocity profile is interpolated with B-splines from which its derivative is calculated. In the 3D-WSS method instead, velocity derivatives are calculated on the 3D lumen segmentation from B-splines interpolation of blood velocity profiles, exploiting only a fixed number of velocity points close to the boundary.
Both methods were tested independently of the lumen segmentation technique. The lumen segmentation was obtained from the test section CAD (SolidWorks 2016, Dassault Systèmes) as an stl file overlapping the test section fluid domain in the anatomy images. Markers on the TS permitted to obtain data from different datasets at the same location. The obtained data were employed as inputs for 2D-WSS. 3D-WSS calculations need to be performed from 3D data sets. For 4D acquisitions, such 3D data sets are readily available, but it is not the case for the two dimensional LDV and 2D PC MRI measurements. Thus, for these two latter, 3D data sets were artificially built by duplicating the 2D measurements in several planes in the x-axis (see Fig. 1b for axes reference), spaced by 2 mm as in the 4D data sets. For 4D Flow MRI data resolution refinement, WSS was calculated only with the 3D-WSS method.
For the comparison between LDV and CFD, WSS was calculated with a first-order finite difference scheme (Eq. 5, µ is the fluid viscosity, s the spatial resolution, v wall and v wall+1 the velocity at the boundary and at the boundary nearest point, respectively).
(3) n air sin α air = n fluid sin α fluid Linear regression Y = aX + b was performed for velocity and WSS profiles, X = [x i , . . . , x N ] and Y = [y i , . . . , y N ] were respectively the N velocity points of LDV and PC MRI datasets. Coefficients a and b were determined for the analysis. Additionally, Root Mean Square Error (RMSE), Coefficient of Determination R 2 , and L2 norm error ǫ L2 were included in the analysis (Eq. 7) Correlation between different velocity or WSS data was assessed with the Pearson Correlation Coefficient (PCC), while consistency was assessedwith Intraclass Correlation Coefficient (ICC). Coefficients a and b, RMSE, R 2 , ǫ L2 , PCC, and ICC were evaluated for each acquired phase and subsequently averaged over the 30 samples. Their scores for the exact correspondence between diverse groups of data are reported in Table 2.

Results
Preliminary verification of LDV with CFD. Results of the LDV verification with CFD are displayed for the pulse test in Fig. 2. Velocity profiles show a good agreement between LDV and CFD, as well as correlation plots with regression lines. WSS profiles display significant discrepancies between CFD and LDV, as CFD did not reproduce the exact same experimental conditions. Specifically, results for values of the boundary coordinate l between 25 and 75 mm might suggest that an imperfect alignment of the flow at the experimental test section inlet occurred.
PC MRI velocity profiles. Results for velocity data centerlines and correlation plots between LDV and PC MRI are reported in Fig. 3 for Pulse4D X0 , Pulse4D X1 , and Pulse4D X2 tests (see Table 1 for test labels), while the quantitative analysis is summarized in Table 3. The effect of the resolution refinement on the magnitude images is displayed in Fig. 2. Even though the refinement improved the image resolution, for the fine resolution, overall image quality decreased. There was globally a good agreement between LDV and PC MRI results. Both 2D PC and 4D Flow MRI reproduced LDV velocity profiles well, and the same feature was found with the refinement of the spatial resolution in 4D acquisitions. 2D PC MRI results were noisier than those of 4D Flow MRI. Noise impacted more data close to the boundaries and at low flow rate phases, as evidenced by the velocity oscillations. Concerning the quantitative analysis, the VNR was higher for 4D Flow than 2D PC MRI, and for steady than  Fig. 4 and in Table 4. Profiles of 3D-WSS calculated from 4D Flow MRI data displayed generally a satisfactory agreement with those obtained from LDV velocity data. The refinement of the voxel size introduced oscillations in 3D-WSS profiles, particularly        www.nature.com/scientificreports/ for t ≥ 0.5T c . Concerning WSS quantitative analysis, overall results for 2D-WSS were inconsistent, particularly for 4D Flow acquisitions, for which figures of mean WSS errors increased importantly. Also, for Steady tests, correlation and consistency were weak and moderate respectively, while for Pulse tests, there was no correlation and consistency was poor. Results of 3D-WSS instead recorded better performances. Figures were in the same order of magnitude for 2D PC and 4D Flow MRI, but the results of pulse acquisitions were slightly worse. For Steady tests, correlation was strong and consistency excellent while, for Pulse tests, consistency was moderate and consistency was good for 4D Flow and moderate for 2D PC MRI. For data resolution refinement, results registered improved performances for steady tests, and scores worsened with the decrease of the voxel size. However, for steady tests, consistency and correlation remained excellent and strong respectively.

Discussions
In this work, an easy-to-use and accurate LDV-based system to validate potential improvements in PC MRI acquisitions and post-processing was designed and manufactured. The reliability of the system was preliminarily assured through the successful assessment of LDV measurement accuracy and errors. The experimental campaign was carried out proficiently, proving the efficacy of the employment of this sort of tool in research on clinical applications of PC MRI. Despite it being time-consuming, LDV enabled accurate measurements of fluid velocity profiles, providing information to test preliminary PC MRI performances and thus flow quantification, as it did for WSS calculation. Even though LDV and PC MRI campaigns were performed independently at different times, on the same test rig but in separated facilities, a good level of reproducibility of the results with the two different techniques was achieved for both velocity and WSS results. Accurate LDV tools represent, therefore, a robust alternative to CFD to simply assess 4D Flow MRI quality. Although CFD is often employed as a reference for PC MRI, its role of ground truth is not widely accepted in the scientific community. For velocity profiles, results generally showed a strong correlation and an excellent consistency between PC MRI and LDV, representing a strong basement for further investigations on post-processing algorithms, as evidenced by statistical parameter scores. Firstly, a method was employed to degrade the LDV data resolution to directly compare LDV and MRI velocity profiles, inducing a small distance shift d shift between LDV and MRI grids, but without a noticeable effect on the analysis ( d shiftmax = 0.2632 mm). No significant differences were noticed between standard 2D and 4D acquisitions, except for flow rate errors that were greater for 4D Flow MRI results since data resolution was lower. As largely discussed, validation of WSS calculation techniques with this experimental set-up represents a crucial clinical aspect. To demonstrate that, the system was proved with two algorithms to compute WSS from PC MRI data. WSS calculation from PC-MRI-based velocity data in this specific application was positively accomplished. Even if in clinical practice PC MRI data processing includes segmentation of anatomical images to obtain the domain in which blood flow quantification is performed, the analysis was simplified by testing the methods independently of the segmentation technique. The improved accuracy of 3D-WSS and its greater performances and versatility compared to 2D-WSS have been properly assessed thanks to this LDV based tool. For 3D-WSS, there were no significant differences between 2D PC and 4D Flow MRI results, while 4D Flow 2D-WSS outcomes were more inaccurate than 2D PC 2D-WSS. This high sensitivity of the 2D-WSS method to the velocity data resolution might be due to the employment of a squared geometry, representing a sensitive case for B-spline interpolation for the corners of lumen contour.
Finally, a concrete application of this approach to solve a relevant clinical challenge in 4D Flow MRI acquisitions and derived WSS calculation was investigated to additionally prove the benefit of relying on LDV reference data to optimize 4D Flow MRI. The influence of voxel size on acquired velocity data and computed 3D-WSS was successfully assessed. For velocity profiles, the agreement between LDV and 4D Flow in terms of flow rate errors, consistency, and correlation globally improved with the data resolution refinement. For the finest data resolution compared to the medium one, this agreement worsened though, which might be attributed to a higher velocity noise consequently to the voxel size reduction. This would also explain why the agreement between LDV and 4D flow MRI worsens when refining the resolution for 3D-WSS profiles, which are more sensitive to noise than velocity profiles. The impact of noise is stronger for low-velocity data, like those close to the boundaries with which WSS is calculated. This might have resulted in a reduction of the accuracy of the WSS assessment caused by a significant increase in noise while refining the spatial resolution.
It can be revealed from this experimental campaign that, as a general rule in clinical practice, for an equivalent acquisition time there is a trade-off to be found between data spatial resolution and VNR, as noise has a stronger impact on WSS calculation with respect to the velocity measurement. For potential clinical applications of current WSS computation techniques, an MRI spatial resolution improved with respect to that of standard protocols might not always imply a more accurate assessment of WSS, if VNR is highly degraded by an increase of noise in the acquired data. In this study, the effects of the temporal resolution were not investigated, but the system could additionally be employed to accomplish this task.
Limitations of the present study included the employment of an unphysiological scenario, in terms of test section shape and flow rate conditions, and the background phase removal. The selection of a squared test section did not represent a realistic anatomical geometry. Nevertheless, it represented a simple and optimal set up for accurate LDV measurements, since the optical interface crossed by LDV laser beams is constituted of a flat-screen, preventing then optical aberrations, which could instead arise from a curved surface. Therefore this squared test section was well suited to investigate the WSS biomarker. The flow rate curves were selected to reproduce the thoracic aorta fluid dynamic conditions in terms of Reynolds ( Re mean ≈ 1100 ) and Womersley ( α ≈ 20 ) numbers, compatibly with the experimental set-up constraints ( Re mean = 1273 for Steady tests; Re mean = 190 and α = 18 for Pulse tests). These unphysiological flow rate conditions led to lower magnitudes of WSS signals, which were thus harder to assess. Nonetheless, the proposed methodology is valid for the higher WSS magnitude found in physiological flow rate conditions. Further experimental campaigns with a more www.nature.com/scientificreports/ realistic scenario of flow rate and vessel geometry could be performed, however, with this experimental device. Concerning the unphysiological flow rate conditions, for steady tests, the approximate average fluid dynamic conditions of the thoracic aorta were achieved. For pulse tests instead, through the Womersley number, only the pulsatile features of blood flow of this specific cardiovascular location were replicated. In any case, even if during pulse tests velocity magnitudes were far from physiological values, this represented a challenging scenario for the assessment of WSS, since velocity gradients were lower than physiological values. Further experimental campaigns with a more realistic scenario of flow rate and vessel geometry will be performed. Lastly, preliminary PC MRI processing procedures to remove the effects of magnetic field inhomogeneities and other artefacts are mandatory in clinical investigations, but in this specific application, they were not implemented. Their effects were evaluated considering velocity data in the agarose gel surrounding the test section, and they were considered negligible in a first approximation. This might be explained by the fact that in this straightforward application the region of interest was located in the center of the scanner where magnetic fields are more homogeneous, and the FOV was properly centered in this region. Notwithstanding no correction was applied, results showed a good agreement with other experimental references (LDV, flowmeter), and WSS assessment was properly conducted. In this first preliminary phase, this aspect was not considered crucial, but it will be considered for any other further investigations involving this experimental set-up.
This work represents a strong foundation for future validation of novel developments for blood flow quantification from PC MRI acquisitions. Particularly, the method based on laser Doppler velocimetry presented here might provide reliable information about the accuracy, reproducibility, and applicability of the techniques employed to calculate in vivo WSS from 4D Flow MRI data. Based on the findings of this work, future developments of the present study are represented by the design of a novel hydraulic rig under more realistic conditions in terms of physiological and pathological anatomies, flow rates, and pressure. It could be employed as an independent standard to validate 4D Flow MRI acquisition sequences and post-processing, particularly regarding quantitative hemodynamic biomarkers for blood flow quantification, such as the wall shear stress, pulse wave velocity, or pressure mapping.