Time-Resolved Wall Shear Rate Mapping Using High-Frame-Rate Ultrasound Imaging

In atherosclerosis, low wall shear stress (WSS) is known to favor plaque development, while high WSS increases plaque rupture risk. To improve plaque diagnostics, WSS monitoring is crucial. Here, we propose wall shear imaging (WASHI), a noninvasive contrast-free framework that leverages high-frame-rate ultrasound (HiFRUS) to map the wall shear rate (WSR) that relates to WSS by the blood viscosity coefficient. Our method measures WSR as the tangential flow velocity gradient along the arterial wall from the flow vector field derived using a multi-angle vector Doppler technique. To improve the WSR estimation performance, WASHI semiautomatically tracks the wall position throughout the cardiac cycle. WASHI was first evaluated with an in vitro linear WSR gradient model; the estimated WSR was consistent with theoretical values (an average error of 4.6% $\pm ~12.4$ %). The framework was then tested on healthy and diseased carotid bifurcation models. In both scenarios, key spatiotemporal dynamics of WSR were noted: 1) oscillating shear patterns were present in the carotid bulb and downstream to the internal carotid artery (ICA) where retrograde flow occurs; and 2) high WSR was observed particularly in the diseased model where the measured WSR peaked at 810 $\text{s}^{-{1}}$ due to flow jetting. We also showed that WASHI could consistently track arterial wall motion to map its WSR. Overall, WASHI enables high temporal resolution mapping of WSR that could facilitate investigations on causal effects between WSS and atherosclerosis.

directly related to carotid plaque development as well as its progression and possible rupture [1], [2]. This quantity is essentially a hemodynamic drag force exerted by viscous blood flowing tangentially to the vessel wall. WSS of near-zero levels is now accepted as a biomechanical factor that favors plaque development and its progression toward rupture [3], [4]. In contrast, high and focally elevated WSS levels may pose an increased risk of plaque rupture and intraplaque hemorrhage [5], thereby increasing the chance of stroke.
In practice, WSS can be estimated by measuring the wall shear rate (WSR) since WSS and WSR are related through a blood viscosity coefficient, which generally remains constant for Newtonian flow in medium and large arteries [6]. By definition, WSR is the tangential flow velocity gradient at the arterial wall, and it can be determined by measuring the intraluminal flow velocity profile through noninvasive means. Accordingly, WSR is a more accessible measurand than WSS. A few techniques for WSR mapping have been reported, such as computational fluid dynamics (CFD) simulations [7] and magnetic resonance imaging (MRI) [8], [9]. Nevertheless, CFD is not favorable to clinical workflow because it requires a 3-D image volume scan (e.g., using MRI) to construct a simulation model and needs patient-specific specification of boundary conditions that remains clinically challenging to acquire [10]. For MRI, its long scan time and the massive size of scanners have made this modality difficult to be used in bedside settings. In addition, the inherently limited spatial resolution of MRI has consistently led to WSR magnitude underestimations [11], [12].
To acquire vascular WSR information on a more routine basis, ultrasound is perhaps a more suitable noninvasive modality since it can be used in point-of-care scenarios [13]. Hitherto, a few ultrasound techniques have been proposed for WSR estimation. However, the feasibility of using ultrasound to robustly measure WSR over the entire field of view on a time-resolved basis has not been demonstrated. For methods based on single-gate [14], [15] and multigate Doppler principles [16], [17], [18], [19], they inherently lack the provision of image-level WSR information due to the use of a single ultrasound beam for measurements [20]. Conversely, methods based on the color flow imaging paradigm [21], [22] suffer from limited temporal resolution (i.e., inadequate frame rate, which is about 20 frames/s at best [23]) that is crucial to resolve how WSR evolves over various pulse cycle phases. On the other hand, contrast-enhanced image velocimetry does offer WSR maps at high temporal resolution [24], [25], but it after all requires microbubble injection that disfavors routine clinical use. More recently, the derivation of WSR has been attempted directly from vector flow images acquired from clinical scanners [26], [27]. Yet, their efficacy is inherently susceptible to near-wall flow estimation biases due to clutter filtering [28] and inaccurate wall tracking.
In this article, we present a new noninvasive WSR mapping technique that makes use of high-frame-rate ultrasound (HiFRUS) principles to track the spatiotemporal evolution of WSR in an artery segment on a time-resolved basis. Our new method, which we refer to as wall shear imaging (WASHI), is designed to robustly derive WSR values at millisecond time resolution over different arterial wall positions in the imaging view. The rest of this article will explain WASHI's theoretical principles (see Section II) and technical implementation details (see Section III). Its performance in achieving spatiotemporal profiling of WSR dynamics will also be presented in the context of phantom experiments and an in vivo pilot trial (see Sections IV and V).

II. THEORY AND TECHNICAL PRINCIPLES A. Algorithmic Overview
WASHI generally works by making use of HiFRUS to perform, at every slow-time sampling instant, WSR estimation at all detected vessel wall positions. This WSR estimation process involves three specific steps as shown in Fig. 1: 1) tracking of vessel wall positions in the imaging view; 2) derivation of the imaging view's intraluminal flow velocity map; and 3) calculation of the tangential flow velocity gradient at detected vessel wall positions as per the definition of WSR. Vessel wall positions are tracked in the first step using a semiautomatic, cross correlation-based detection algorithm on each HiFRUS image frame (see Section II-B). In the second step, because flow velocity is essentially a vector quantity, the intraluminal flow velocity map is derived by performing flow vector computation at detected blood vessel pixels using multi-angle Doppler estimation principles (see Section II-C). To obtain a consistent intraluminal velocity map, a regularization routine is implemented to smooth the flow vector estimates. In the final step, WSR values are computed by performing a 2-D spatial derivative on the intraluminal velocity map and extracting the corresponding tangential velocity gradient values at detected wall positions (see Section II-D).

B. Wall Tracking
As vessel walls are in motion due to the pulsatile nature of blood flow in arteries, the wall positions vary over a cardiac cycle. The wall positions, as an input parameter to the WSR estimator, should be dynamically tracked to accurately derive the WSR. In this work, we used a semiautomatic wall tracking technique to dynamically update the wall positions in the sequential image frames. It works by first manually identifying the vessel wall positions (depicted as gray square markers in Step 1 of Fig. 1) in the reference frame, which was essentially the first frame of the B-mode image sequence.
In the second step of the wall tracking procedure, the displacements of the wall positions in the subsequent B-mode image frames are obtained through a speckle tracking algorithm. Specifically, for each labeled wall pixel, we defined a search template that is centered upon the labeled pixel. New positions of the template over the cardiac cycle are then identified by comparing the template to its adjacent neighbors in subsequent beamformed B-mode image frames using 2-D cross correlation to identify pixels that best match the template. The interframe displacement is calculated to update positions of the corresponding markers, thereby allowing wall positions to be dynamically tracked over all image frames.
Based on the tracked wall displacements, the curvature of the vessel wall is determined to facilitate WSR computations (to be described further in Section II-D). In particular, the tangential angle of a wall segment θ is approximated as the slope between two adjacent position markers on the tracked vessel wall and is updated as the wall moves. This step served to account for the tortuous nature of blood vessels and their local curvature that resulted in varying wall tangential angles at different vessel positions. Note that the deduction of vessel wall curvature is necessary here because shear stress is a directional quantity tangential to the tissue-fluid interface.

C. Intraluminal Flow Velocity Estimation
To derive intraluminal flow velocities, the multi-angle leastsquares vector Doppler technique [29], [30] is used, where the Doppler signal is acquired from multiple beam-flow angles. In brief, a set of M-steered plane wave transmissions are sent to perform pulse-echo sensing over the entire imaging view. For each of these transmissions, parallel beamforming is used to generate multiple image frames with N receive steering angles. This results in MN image frames, each with a unique transmit-receive (Tx-Rx) angle pair, for a set of transmissions. By replicating this procedure for subsequent repeated sets of steered plane wave transmissions to acquire the Doppler signals, MN sets of image frames are generated. For each set of these image frames, slow-time processing is performed on a per-pixel basis, where the Doppler signal of each pixel is processed with clutter suppression, followed by the mean frequency estimation. Finally, for each pixel of the velocity map, the mean frequency estimates from all angle pairs collectively form an MN × 1 measurand vector u, and the axial and lateral velocities (v z and v x ) are computed by multiplying u with the pseudoinverse of an angle matrix A as where the superscript "T " denotes matrix transpose and A is an MN × 2 matrix containing the trigonometric terms related to the Tx-Rx angle pairs as described in our previous works [29], [30].
Recognizing that flow vector estimation bias may emerge near the vessel wall [29], we devised a velocity regularization framework to improve slow-flow estimation at regions near the vessel wall. This framework is important as inaccurate velocity estimates are known to significantly impact the WSR estimation performance [19]. To illustrate such an issue, the velocity profile along an axial line for an estimated velocity map [see Fig. 2(a)] is plotted in Fig. 2(b). We observed that the flow profiles at the near-wall regions are noisy and deviate from the no-slip boundary condition (i.e., zero velocity at the wall). Our regularization framework utilizes the wall positions information obtained from Section II-B (Step 1 of Fig. 1) to first identify tissue regions with no flow and imposes velocities at these regions to be zero, as shown in the second row of Fig. 2. Next, flow estimates near the wall regions (areas shaded in cyan on the third row of Fig. 2) are discarded. Velocities at these regions are then replaced with values interpolated from a spatially smoothed velocity map using a penalized least-squares regression algorithm [31], [32] that preserves the overall flow profile, thereby yielding the velocity map in the bottom row of Fig. 2.
In principle, our regularization procedure functions to smooth the velocity vector field that may contain the outliers while recovering the missing data. The smoothing process is performed using discrete cosine transform-based regression [31] that penalizes outliers by assigning a lower weight to outlying values in order to reduce their influence. In the case of the near-wall velocities, discarded velocities [shaded area in Fig. 2(f)] at the near-wall region become zero-weighted to exclude their influence on fitting. The degree of smoothness, termed the smoothing parameter, is automatically determined from the minimization of the generalized cross-validation (GCV) score that seeks to quantify the smoothness in splines [33]. Thus, the overall regularization algorithm can be automated and is independent of operator tuning. The regression procedure is iteratively performed; for each iteration, the algorithm determines new weights based on the residuals, and the weights are updated until the residuals remain unchanged, while the GCV converges to its minimum. These steps essentially perform interpolation at the near-wall region while preserving the local smoothness of the overall velocity profile. For each image frame, the regularization framework is applied to all velocity maps with the updated wall positions at different time points to generate a time series of regularized velocity maps for WSR estimation.

D. Derivation of WSR
WSR is derived as the velocity gradient tangential to the vessel wall. Instead of directly examining flow normal to the wall, we have taken the approach to calculate the velocity gradient from the respective contributions of the axial and lateral velocities before realigning them to the wall curvature. This procedure is shown in Fig. 3. First, velocity gradients are derived by convolving the regularized (axial and lateral) velocity maps with directional (horizontal and vertical) 2-D Sobel kernels to produce four spatial gradient maps. The four maps in combination represent the velocity gradients of the axial and lateral velocities when interrogated in both the horizontal and vertical directions. By combining the four spatial gradient maps, the 2-D stress tensor matrix shown in Fig. 3(c) is derived. Finally, to obtain the local WSRs, the stress tensors are aligned to the wall curvature by rotating the tensors with their corresponding wall angles. It should be noted that these operations are only performed on the known wall coordinates deduced in Section II-B. Details of each step are discussed as follows.
1) Computation of Velocity Gradient and Stress Tensor: The 2-D Sobel kernel is a spatial derivative filter [34] that performs gradient computation in one dimension while averaging samples in the orthogonal dimension. In this work, the size of the Sobel kernel is expanded by convolving the Sobel operator with an orthogonal smoothing kernel, as shown in Fig. 3(a). The smoothing kernel essentially functions as a center-weighted spatial averaging filter. The 2-D Sobel kernels are oriented in the vertical (K z ) and horizontal (K x ) directions [see Fig. 3(b)], formulated based on the orientation of the gradient kernel. Using these two forms, the velocity gradient maps G can be computed by convolving the velocity maps with the Sobel kernels that are mathematically defined as follows: where G xx and G xz are the horizontal and vertical gradient map of the lateral velocity v x , respectively, G zx and G zz are the horizontal and vertical gradient maps of the axial velocity v z , respectively, and x and z are the lateral and axial pixel sizes, respectively. Because the convolution operation is only performed along the vessel wall, the kernel is centered at the blood-wall boundary and a factor of 2 is introduced to compensate for the zero velocities in the tissue region such that continuity is maintained at the velocity field when computing the differential operations. Together, these gradients form the stress tensor τ as where the diagonal elements σ xx and σ zz correspond to the normal stresses (i.e., pressure gradients), whereas the offdiagonal elements τ xz and τ zx represent the shear stresses. There are two advantages to using a 2-D Sobel operator to derive WSR. First, the implementation of the 2-D Sobel operator is unaffected by the blood vessel's geometry to directly derive velocity gradients from the velocity maps; adjustments in response to the vessel curvature can be more efficiently performed in a later step. Second, attributed to its averaging property, the velocity gradient estimation is less susceptible to fluctuations present in the velocity map to produce a consistent gradient estimate. The size tradeoff of the Sobel kernel should be noted: a large kernel, although reduces measurement variance, would introduce an underestimation bias since the spatial gradient decreases toward the lumen center. Note that, in lieu of the Sobel operator, more specialized spatial derivative operators such as the Scharr operator may be used if preferred.
2) Derivation of Tangential Wall Shear Component: Local WSR is computed by first realigning the stress tensor to the local wall curvature approximated in Section II-B before extracting its shear component. Fig. 3(c) shows an illustration of the free body diagrams for the stress tensors τ and τ (with their stress components), respectively, before and after realignment. The shear stress components (τ xz and τ zx ) of τ in (3) are defined with reference to the image frame instead of being tangential to the vessel wall. Realignment of the stress tensor to its local tangential angle θ (estimated as described in Section II-B) is achieved by applying a rotation matrix Q using the following transformation: Upon applying the rotational transform, the rotated shear stress component (τ zx ) of τ , lower left element, is now tangential to the wall. By repeating this procedure for all the derived stress tensors along the vessel wall to their corresponding local tangential angle θ via the rotational transformation in (4), the local WSR map can be acquired.

E. Visualization
To present the derived WSR, vessel wall pixels are colorencoded to their corresponding WSR magnitude. The visualization of WSR is further complemented by flow vectors seeded using our vector projectile imaging (VPI) technique reported previously [29] for users to synchronously perceive the flow dynamics as well as appreciate the influence of flow on the wall shear. The final WASHI image frame is comprised of three stacked layers as shown in Fig. 4(a): 1) B-mode image displayed in the background to provide anatomical information; 2) a color-encoded WSR map rendered at the foreground layer superimposed onto coordinates of the vessel wall based on information acquired in Section II-B; and 3) flow projectiles rendered within the vessel lumen to depict the flow direction.
The WSR map is rendered with a bicolor hue with the colormap's brightness proportional to the absolute value of the derived WSR. A bicolor hue representation is selected to imply the flow direction at the near-wall region, as shown in Fig. 4(b). When blood flows from left to right, both the upper and lower walls are mapped with a cyan hue. In contrast, the vessel wall is mapped with a magenta hue for a right-to-left flow. In the event of counter-directional flow occurring within the same vessel segment, such as the flow recirculation in the bifurcation site in Fig. 4(a), the upper and lower walls are mapped with different hues. The use of such bicolor hues serves to highlight flow irregularities. This WSR visualization strategy, when coupled with dynamic visualization of flow projectiles, allows users to more intuitively perceive the WSR pattern of interest.

III. IMPLEMENTATION OF WASHI
A. Imaging Hardware and Data Acquisition WASHI was implemented on a research purpose, channeldomain imaging research platform (SonixTouch and L14-5 linear array; Analogic Ultrasound, Peabody, MA, USA) with a multichannel pre-beamformed data acquisition tool [35] and software beamforming [36]. In this work, a set of plane wave pulses with three steering angles (i.e., M = 3 at −10 • , 0 • , and 10 • ) was used; for each transmit event, a three-cycle, 5-MHz pulse was transmitted to insonify the region of interest. The plane wave transmission was repeated at a pulse repetition frequency (PRF) of 10 kHz, which yielded a 3.3-kHz effective PRF for each transmit angle. The peak mechanical index was 0.414 with the spatial-peak temporal average intensity of 214.6 mW/cm 2 . Image acquisition was performed for 2 s over a 6-cm imaging depth to sufficiently cover one full cardiac cycle. The acquired RF data were streamed offline to a backend computing workstation that leverages the graphics processing unit (GPU) technology [36].

B. Image Formation for Wall Tracking and Flow Vector Estimation
Image frames were generated from the acquired radio frequency (RF) data by evoking our GPU codec library [36] programmed in MATLAB (2016a, Mathworks Inc., Natick, MA, USA). For each pulse-echo event, a 64-tap finite-impulseresponse (FIR) bandpass filter (3-7 MHz) was first applied onto the RF data on a per-channel basis to suppress out-ofband noise. Next, the RF data were converted to its analytic form using an FIR-based Hilbert transformer. The analytic RF data were subsequently beamformed using the delay-andsum algorithm to form low-resolution images (LRIs) with a pixel size of 0.1 × 0.1 mm. LRIs from the same set of plane wave excitations were coherently compounded to form a high-resolution image (HRI). These HRIs were used for wall tracking (described in Section II-B) as well as to serve as the B-mode background layer for the visualization of WASHI.
For the purpose of flow vector estimation, the same analytic RF data were beamformed again but with three receive steering angles (i.e., N = 3 at −10 • , 0 • , and 10 • ) at a pixel size of 0.1 × 0.1 mm. By implementing this beamforming procedure for all RF data, MN = 9 sets of image frames were generated, each with a unique Tx-Rx angle pair. This configuration has been shown to provide consistent flow vector estimation performance [29] which is critical for WSR estimation.

C. Semiautomatic Wall Tracking
To facilitate automatic wall tracking, the algorithm was initialized with the wall position of the first HRI frame by manually placing markers on the B-mode image using MATLAB's built-in region-of-interest selection function. The markers were placed on the brightest line of the wall, which depicts the interface between vessel wall and lumen. Since the physical image resolution was near 1 mm for our given imaging parameters, each marker was roughly distanced 1 mm apart to finely trace the curvature of the vessel wall. Positions of these wall markers were updated over time by tracking their new positions using a cross correlation algorithm [37], repeated for every eight HRI frames (matching the sliding step size of flow vector processing).
To improve the wall tracking resolution, each HRI was first upsampled by a factor of 5 using spline interpolation in the axial and lateral directions, resulting in a pixel size of 0.02 × 0.02 mm. Next, a 2 × 2 mm speckle template (centered at the marker) was defined for each wall marker from the first HRI frame. New positions of each marker were determined by cross-correlating the speckle template with surrounding pixels in the subsequent frames within a search area of 4 × 4 mm. Locations that best matched the template (highest cross correlation value) were identified as the new wall marker positions. This block-matching procedure was repeated for every eight HRI frames to track the wall position over time and the wall coordinates were updated accordingly. At the same time, the tangential wall angle θ defined by the curvature of the vessel wall was calculated from the slope between two adjacent markers.

D. Flow Velocity Map Derivation
For each of the nine Tx-Rx angle pair image sets, pixelbased slow-time Doppler processing was conducted as previously reported [29], [30] to derive the flow vector maps required for WSR estimation. Briefly, for the same pixel location on each set of the nine unique Tx-Rx beamformed frames, a 134th tap FIR highpass filter (cutoff at 66 Hz) was first applied along the slow-time direction to suppress tissue clutter. Next, the lag-one autocorrelation phase algorithm [23] was performed for an ensemble of 64 samples to derive the mean Doppler frequencies. This treatment yielded MN = 9 mean Doppler frequencies for the same pixel location. The Doppler frequencies were then grouped together to form the measurand vector u and multiplied to the pseudoinverse of angle matrix A to derive the flow vector v according to (1), as detailed in Section II-C. This flow vector estimation was repeated for all pixels to derive the full-view axial and lateral velocity maps (V z and V x ), which collectively form the flow vector map. Subsequent flow vector maps were generated for every eight slow-time samples (87.5% overlap with the previous frame) to produce a finer temporal representation of the flow dynamics.
To mitigate the overestimation of flow velocities at near-wall regions as an effect of clutter filtering [38], regularization (as described in Section II-C) on the derived flow vector maps was performed by first upsampling the velocity maps with spline interpolation by a factor of 5. Then, velocities in tissue regions (determined from the wall positions) were assigned to 0 m/s. Finally, flow that was present within 0.4-mm proximity (corresponds to the axial resolution of a three-cycle pulse) to a vessel wall was discarded and was instead estimated using the penalized least-squares regression algorithm [31], [32]. This regression algorithm also spatially smoothed the velocity maps to suppress inconsistent fluctuations while preserving the overall flow profile, as shown in Fig. 2(g) and (h).

E. WSR Estimation
WSR estimation was carried out by computing the four velocity gradient components (G zx , G xx , G zz , and G xz ) to form the stress tensor, followed by rotating the tensor to align with the wall curvature. For all pixels identified on the bloodtissue interface, the gradient maps of G zx , G xx , G zz , and G xz were derived by convolving the horizontal and vertical Sobel kernels (11 × 11 pixels, covering 0.2 × 0.2 mm area) separately with the axial and lateral velocity maps according to (2). For each of the wall pixels, the corresponding stress tensor matrix τ was constructed from the gradient maps according to (3). The tensor matrix τ was subsequently multiplied to the rotation matrix Q, formed using the corresponding wall angle θ according to (4) to align the stress tensor to the wall curvature. The lower left element of the rotated matrix τ corresponded to the WSR of that wall pixel. The procedure was repeated for all wall pixels for all acquired image frames, resulting in a time series of WSR mapping that evolves with the wall motion and flow dynamics.

F. Visualization of Wall Shear Mapping
To visualize the estimated WSR maps, a triplex rendering scheme as described in Section II-E was used. The background B-mode HRI was constructed by coherent summation of three LRIs from different Tx angles. Next, flow information was rendered on top of the B-mode image using our VPI technique [30] by displaying randomly seeded projectiles within the flow region. The length of the displayed projectiles represented the local instantaneous flow velocities, while their color was pegged to its lateral velocity to highlight recirculation. In between display frames, positions of the projectiles were shifted based on their respective displacements calculated as the product of local velocities with a 2.4-ms interframe time period (sliding step size divided by PRF). The length and the color of the projectiles were updated to reflect the new local velocities. Finally, the WSR estimates were overlaid onto the vessel wall with their color pegged to a bi-hue color scale dedicated to WSR representation. Note that the displayed color hue of WSR is indicative of the flow direction as described in Section II-E and Fig. 4 such that the presence of flow irregularities is apparent. The plotted markers were sufficiently dense to form a continuous line over the vessel wall.

IV. EXPERIMENTAL METHODS A. Hele-Shaw Flow Calibration
To assess the performance of WASHI, a Hele-Shaw flow phantom was designed to exhibit a linearly decreasing WSR pattern. The geometry of this linear-shear Hele-Shaw model was adapted from previously established works [22], [39] and was fabricated to be acoustically compatible by modifying our lost-core casting protocol [40]. In comparison to straighttube models, the linear-shear model offers a range of WSR for the evaluation of WASHI, whereas the former's approach provides only a single WSR value [25], [41]. Reference WSR values of the Hele-Shaw model can be obtained by solving the Navier-Stokes equation.
1) Hele-Shaw Model Design: The flow chamber of the Hele-Shaw model (shown in Fig. 5) was designed using computer-aided design (CAD) software (Solidworks; Dassault Systems, Waltham, MA, USA). The flow chamber has a uniform height h of 3 mm with an expanding width w increasing from 6 to 18 mm. The progressive increment of w can be expressed as a function of its lateral position l for 0 ≤ l ≤ 80 mm where w 0 is the entrance width (6 mm) and L was set to 120 mm such that w would be three times of w 0 at l = 80 mm. The linearly decreasing WSR pattern is an effect of the decreasing flow velocity in response to the gradual widening of the flow chamber. To facilitate the connection of the flow chamber to a flow circuit, an adaptor was appended to the entrance of the chamber, as shown in Fig. 5 (shaded in  yellow). This adaptor has an inlet size of 3 × 3 mm, in which its width gradually widens to match that of the flow chamber. The overall length of the adaptor is 80 mm to deliver a smooth transitional flow from the flow connector to the flow chamber.
2) Phantom Fabrication: To fabricate this phantom, the drafted flow chamber was saved as a stereolithography file (STL) and was converted to print instructions at 0.1-mm layers using a slicer software (KISSlicer ver 1.6). The print instructions were executed on a 3-D printer (Model DX; Cre-atBot 3-D Printer, Zhengzhou, China), which has a 0.25-mm nozzle, to construct the flow chamber replica. The physical replica was later suspended in a casting box (160 × 80 × 40 mm), held by two supporting plates on both ends. Next, the polyvinyl alcohol (PVA) solution was poured into the casting box to form the tissue mimic. The PVA solution contained 10% PVA (341584; Sigma-Aldrich, St. Louis, MO, USA), 3% graphite (5 μm-sized; 282863; Sigma-Aldrich) as ultrasound scatterers, 0.3% potassium sorbate (85520; Sigma-Aldrich) as preservative, and 86.7% distilled water. Three 24-h freeze-thaw cycles were administered for the PVA solution to congeal, yielding a cryogel with acoustic speed and attenuation coefficient of 1535 m/s and 0.229 dB/(cm · MHz), respectively [40]. The flow chamber was subsequently removed from the cast to instate the flow chamber. The Hele-Shaw flow phantom was placed in a four-sided box; the top of the box was open for imaging accessibility, while the outlet wall was removed to facilitate fluid release. The phantom's inlet was connected to a flow pump delivering blood-mimicking fluid (BMF) that has a viscosity of 4.1 mPa·s, matching that of blood. During experiments, the flow pump was set to operate at three constant flow rates of 1.5, 3, and 4.5 mL/s.
3) Imaging Experiments: WASHI was performed to map the WSR, as described in Section III. The ultrasound probe was placed at the centerline of the Hele-Shaw flow chamber, 6 cm from the phantom's inlet as annotated as the dashed line in Fig. 5. Flow within this imaging view is expected to be unidirectional along the lateral axis, with the property of linearly decreasing WSR. To generate the reference WSR values for our Hele-Shaw flow phantom, CFD simulations were performed using the Flow Simulation package of Solidworks (Dassault Systems). The CAD model of the flow chamber was first meshed into approximately 130 000 cubical cells with a cell dimension of 0.13 mm in all directions. The following conditions were defined for the simulation environment: 1) the flow was assumed to be laminar; 2) the fluid viscosity was set to 4.1 mPa·s; 3) the inlet boundary condition was set to be constant volumetric flow with uniform distribution; 4) the outlet boundary condition was set to be atmospheric pressure release; and 5) no-slip condition was enforced at the walls of the flow chamber. The CFD simulation was deemed complete when the average pressures at the inlet and outlet converged. Simulation was repeated for flow rates identical to those of the in vitro Hele-Shaw flow experiments (1.5, 3.0, and 4.5 mL/s) and indicated that WSR in our experimental setup ranged from 100 to 500 s −1 . This WSR range is equivalent to a WSS range of 0.41-2.05 Pa (for the fluid of 4.1-mPa·s viscosity), covering both atherogenic and atheroprotective regions [22].

B. Carotid Bifurcation Phantom Investigation
To evaluate WASHI's ability to map and visualize WSR dynamics under pulsatile flow scenarios, imaging experiments were performed on two anthropomorphic bifurcation phantoms devised previously [42]. The two phantoms consisted of a healthy carotid artery bifurcation and a diseased model with 50% eccentric stenosis at the internal carotid artery (ICA) branch entrance. During experiments, the phantoms were connected to a flow circuit equipped with a programmable flow pump [40] to circulate BMF according to the carotid flow waveform with a pulse rate of 72 bpm and a peak flow rate of 10 mL/s. The ultrasound probe was aligned to the long axis of the phantom to perform WASHI at the planar view of the bifurcations. Data over two complete cardiac cycles were acquired according to the procedure described in Section III-A.

C. In Vivo Pilot Experiments
A pilot study was carried out on a 27-year-old healthy volunteer (with no history of cardiovascular disease) that responded to our open subject recruitment to demonstrate the in vivo feasibility of WASHI. This study has been approved by the Clinical Research Ethics Committee of the University of Waterloo (Protocol No. 31694). WASHI was performed with the ultrasound probe positioned over the longitudinal plane of the common carotid artery (CCA) to capture 2 s of data providing WSR dynamics over two full cardiac cycles. The subject was requested to hold the breath during data acquisition to limit motion artifacts. The CCA branch was imaged here for its greater pulsatility compared to the internal and external branches. This imaging scenario serves as an evaluation for our framework in adapting to the wall motion dynamics, specifically on the effectiveness of the wall tracking algorithm. WASHI processing identical to those applied to in vitro experiments was performed to map and visualize WSR for analysis.
For comparison, a WSR estimation was performed directly on the velocity vector field at the blood-wall boundary determined from the Doppler power map. This step was taken to demonstrate the importance of precise wall segmentation and tracking, so as to investigate their implications on WSR estimation. To generate the Doppler power map, we derived the average power of the wall-filtered slow-time signal (across the 2-s acquisition) for all nine Tx-Rx angle pairs at each pixel location. The power threshold of the Doppler power map was manually tuned, prioritizing the inclusion of nearwall flow pixels, to discriminate flow pixels from tissue. The final map was then converted into a binary mask with its edges identified as the static blood-wall boundary. WASHI processing steps (regularization and WSR derivation) were then implemented on the same in vivo dataset using the static blood-wall boundary as the fixed wall positions. Points of interest along the CCA were analyzed for performance comparison with WASHI.

A. Consistent Estimation for a Spatially Varying WSR Field
WASHI demonstrated its capability to consistently map spatially varying WSR field at high accuracy. When evaluating WASHI using the Hele-Shaw model, WASHI was able to trace the declining WSR in response to the decreasing flow velocity along the lateral dimension seen in the WASHI frames in the top row of Fig. 6 for the tested constant flow rates of 4.5 mL/s [ Fig. 6(a)], 3 mL/s [ Fig. 6(b)], and 1.5 mL/s [ Fig. 6(c)]. Visually, the brightness of the cyan hue, indicating the magnitude of the estimated WSR, dims gradually with the laterally decreasing flow velocity (inferable from the vector projectiles). Further substantiating our observation, a similar declining trend in the estimated WSR along the lateral dimension was noted in the line plots charted on the bottom row of Fig. 6. These line plots were generated using data from 416 WASHI frames (equivalent to 1 s of acquisition) as quantitative assessment; estimates from the upper and lower walls were plotted in blue and red, respectively, with the range for one standard deviation of measurement shaded. Here, we find WSR estimates from WASHI closely approximating the CFD-reference (indicated by the dashed line). It should be noted that, although the estimation bias and standard deviation were found to increase with a higher flow rate, these values decrease once they were normalized to the higher theoretical WSR as the flow rate increased. Using CFD-derived values as a reference, the normalized WSR estimation bias and standard deviation for 1.5, 3.0, and 4.5 mL/s flow rates were, respectively, 7.2% ± 15.2%, 4.9% ± 12.6%, and 1.8% ± 9.3% (average error = 4.6% ± 12.4%).

B. WASHI-Enabled Visualization of Oscillatory Shear in a Bifurcating Vessel
Movie S1 presents the cine loop of WASHI on the carotid bifurcation flow phantom along with an inset of the common branch's pulsed Doppler spectrogram where the cardiac phase is denoted. During systole, WSR surged as blood flow accelerates, reaching its maximum upon peak systole. After passing the systolic peak, the deceleration of flow led to the onset of retrograde flow at the carotid bulb, initiating a shear reversal on the anterior wall of the ICA branch as indicated by the color flip of WASHI from cyan to magenta. As blood flow continues to decelerate, the recirculatory flow zone further develops in size, extending segments of the anterior wall experiencing reversed shear. The reverse shear diminished along with the dissipation of the secondary flow zone when approaching end systole. However, retrograde shear is once again observed with the re-emergence of flow recirculation during the dicrotic wave. This cine loop suggests the presence of oscillatory shear in the bifurcating vessel as demonstrated by WASHI's capability to visualize these characteristics.
To analyze the dynamics of the WSR pattern at the carotid bulb, Fig. 7(a) plots the time evolution of WSR at the anterior wall [position annotated on Fig. 7(b)] alongside the pulsed Doppler waveform. The plot verified the oscillatory nature of WSR observed in Movie S1 as noted by the multiple polarity flips. Fig. 7(b)-(e) shows the snapshots of Movie S1 at various instants of interest with their corresponding time points marked in Fig. 7(a). Fig. 7(b) was captured during the end of peak systole when retrograde flow and reversed shear were first noted. The reversed WSR peaked 24 ms later due to the rapid deceleration of blood flow and recirculation was seen developing in Fig. 7(c). Then, with the systolic phase coming to an end, the WSR transitioned to a forward shear as the secondary flow in Fig. 7(d) has completely dissipated. This episode of oscillatory WSR reemerged to a smaller extent during the dicrotic phase when a recirculation was shown in Fig. 7(e). For the remaining diastolic phase, WSR at the bulb was sustained at low magnitudes (<40 s −1 ) until the next cardiac cycle emerges. Over a cardiac cycle, the recirculatory zone evolved in size; flow separation points (marked with a red arrow for separation point and yellow arrow for reattachment

C. WASHI Showed the Coexistence of High WSR and Low Oscillatory WSR in a Diseased Bifurcation
WASHI cine loop of the diseased bifurcation model is presented in Movie S2 ; pulsed Doppler spectrogram with the range gate placed at the constricted site was also included.
Distinct from the observation made in the healthy bifurcation, reversed WSR was already noted in the stenosed model (at the anterior wall of the carotid bulb) when systole commenced. Here, the retrograde WSR is a result of the recirculatory flow induced by the post-constriction flow jet. The overall WSR increased with the acceleration of blood flow; the maximum WSR (of 800 s −1 ) was recorded during peak systole at the narrowing entrance to the ICA branch. Meanwhile, the WSR of the ICA wall where the flow jet tail was in contact surged to 310 s −1 . Similarly, the retrograde WSR at the carotid bulb also spiked in response to the increased flow velocity. As blood flow decelerates, the overall magnitude of WSR decreased. At this point in time, reverse shear was detected on the posterior wall at the distal end of the ICA branch as another secondary flow zone was seen developing. This secondary flow zone dissipated at end systole, but the recirculatory flow in the bulb region persisted and was further reinforced by the dicrotic wave. The recirculation in the bulb region was not observed to have truly dissipated even toward end-diastole, which was notably distinct from the case in the healthy bifurcation model.
To gain a clearer understanding on the characteristics of WSR at the two secondary flow zones, the evolution of local WSR over a cardiac cycle was plotted in Fig. 8. The estimated WSR at the carotid bulb is plotted in Fig. 8(a), while Fig. 8(b) plots the WSR at the posterior wall at the distal end of the ICA branch [positions marked in Fig. 8(d)]. WASHI snapshots at time points of interest [marked on the Doppler spectrogram sampled at the CCA in Fig. 8(c)] were supplemented in Fig. 8(d)-(i). We noted that the anterior wall at the carotid bulb predominantly experiences reverse shear throughout the entire cardiac cycle, which implies that the recirculatory flow zone never truly dissipated. On the other hand, the first instant of retrograde shear at the ICA distal wall was only spotted after peak systole [see Fig. 8(d)]. The WSR magnitude at the bulb was greatest (−140 s −1 ) during postsystolic flow as blood flow decelerates, which coincides with the maximum secondary flow zone seen at the ICA branch shown in Fig. 8(e). Upon reaching end systole [see Fig. 8(f)], the recirculation in the bulb region persisted to induce low magnitude shear, but the secondary flow in the ICA branch was seen dissipating and the distal wall no longer experiences reverse shear during the dicrotic notch [flow eddy not observed in Fig. 8(g)]. The WSR at the bulb surged once again as flow accelerated during the dicrotic wave and fortified the recirculation [secondary flow zone at the bulb region enlarged in Fig. 8(h)]. Transitioning post-dicrotic wave, blood flow decelerated and returned to forward flow; WSR at the ICA distal branch hovered above zero WSR as noted in the snapshot presented in Fig. 8(i). Throughout diastole, WSR on both lines has a relatively low magnitude, but shear at the bulb region mostly remains in the reverse direction.

D. WASHI Can Adapt to Wall Distension in Vivo
WASHI cine loop acquired at the common carotid is presented in Movie S3 to serve as an efficacy demonstration of our framework in vivo. In general, a resembling trend similar to WASHI cine loops acquired on phantoms could be observed-WSR surges during systole and dicrotic wave. In addition, reverse WSR was briefly noted during the deceleration of blood flow after passing peak systole with the introduction of retrograde flow; the still frame when the secondary flow was detected is included in Fig. 9(a).
Significant arterial wall motion was observed at a greater extent than phantom studies. WASHI has demonstrated its capability in adapting to the tissue motion in Movie   S3 with the framework correctly overlaying the WSR estimates onto the distending wall throughout the entire cardiac cycle. To assess the reliability of the wall tracking component of WASHI, we generated the M-mode plot composed of the background B-mode line and the foreground WSR trace for the line position annotated in Fig. 9(a). In this case, the WSR trace highlights the tracked wall position by WASHI over successive frames. Fig. 9(b) reveals that the WSR trace accurately outlined the arterial wall motion over the entire cardiac cycle, including the sharp distension during systole, indicating that wall motion is reliably accounted for in WASHI.
We examined the estimated WSR at different positions on the CCA, specifically to investigate the discrepancies between WSR on the anterior and posterior walls. Six time traces of WASHI estimates were plotted in Fig. 9(c), (e), and (g) for the positions at the proximal side, midsection, and distal end [exact positions annotated in Fig. 9(a)]. Overall, the time traces show a similar trend, but their magnitudes vary considerably between locations. On the proximal side, the two WSR traces [see Fig. 9(c)] are agreeable, with peak WSRs of approximately 1100 and 900 s −1 on the anterior (dark line) and posterior walls (gray line), respectively. The maximum WSR at the midsection decreased significantly to 400 and 500 s −1 , as shown in Fig. 9(e). Also, reverse shear is noted at t = 0.22 s on the anterior wall, but not at the posterior wall. Finally, we saw a vast difference in the estimated WSR at the distal end [see Fig. 9(g)], with the posterior wall experiencing close to 150% more shear than its corresponding position on the anterior wall during peak systole. In general, WSR of the anterior wall (inner curvature) decreased from 1100 s −1 at the proximal end to 400 s −1 at the distal end, whereas WSR of the posterior wall (outer curvature) remained mostly unchanged on both ends but was reduced at the midsection.
The corresponding WSR measurements at the same lateral positions derived under the static wall positions assumption were plotted in Fig. 9(d), (f), and (g). The WSR time traces exhibited a similar trend as those derived from WASHI and highlighted the differences in WSR between anterior and posterior walls. However, the measured peak WSRs were consistently greater than those measured from WASHI, especially in the midsection region with as much as 102% and 70% increments for the anterior and posterior walls, respectively. At the distal end, a significant difference in peak WSR (70%) was observed for the anterior wall but not at the posterior wall. The peak WSR measurements for both techniques are agreeable at the proximal side of the artery.

VI. DISCUSSION A. Summary of Contributions
Routine clinical monitoring of WSR and, in turn, WSS could improve carotid plaque management since WSS plays a major role in carotid plaque development and its progression [1], [2]. Hitherto, clinical WSR and WSS assessment is limited to heavyweight imaging modalities such as MRI [8], [9] that has a long scan time and provides limited temporal information. In this work, we proposed an ultrasound-based WSR imaging framework called WASHI that offers high spatiotemporal resolution with the potential to be used in a point-of-care setting. WASHI was built upon the broad-view imaging paradigm that can acquire images at high frame rates (over 1000 frames/s), enabling WASHI to derive WSR maps at millisecond resolution. Specifically, WASHI made use of the multi-angle least-squares estimation technique [29] to derive flow velocity maps (see Fig. 1). Coupled with manually identified wall positions that were automatically tracked over time, the near-wall velocities were regularized to suppress overestimation and noise (see Fig. 2). Spatial derivatives were then applied to the velocity maps along the predetermined wall positions using Sobel operations (see Fig. 3) to derive the WSR maps. The derived WSR maps were subsequently visualized using a bi-hue representation to indicate the intensity and direction of the wall shear (see Fig. 4). With the inclusion of vector projectiles, the WSR rendering formed a triplex visualization to intuitively depict the interaction between flow and shear. Compared to other ultrasound WSR estimation methods, WASHI has overcome the limitations of previous works by offering WSR mapping at full view [19] without the need for ultrasound contrast agents [25]. In addition, WASHI was devised to address flow estimation inaccuracies at the near-wall region, which was not accounted for when WSR is directly derived from vector flow images [26], [27].
A series of experiments was conducted to evaluate the performance of WASHI. First, an in vitro setup consisting of a Hele-Shaw flow model that generated a linearly varying WSR pattern (see Fig. 5) was devised to characterize the WASHI's estimation. The results showed that WASHI can accurately derive the WSR (<7.2% bias) for the tested range (see Fig. 6). Next, we assessed WASHI's capability to resolve spatiotemporal dynamics of WSR in flow phantoms of the carotid bifurcation, specifically comparing between a healthy (see Movie S1 ) model against a 50% stenosed diseased (see Movie S2 ) model. While low magnitude oscillatory WSR was observed in the bulb region of both models (see Figs. 7 and 8), focally elevated WSR (over 800 s −1 ) was detected only at the constricted site in the diseased model. More importantly, WASHI captured the multiple occurrences of WSR directional changes over a cardiac cycle. This finding is substantiated by the time trace plots of the local WSR, as shown in Figs. 7 and 8. Finally, the in vivo applicability of WASHI was demonstrated on a healthy volunteer (see Movie S3 ), indicating that WASHI can consistently derive WSR for arteries that are highly distensible (see Fig. 9). This series of work suggests that WASHI is capable of providing consistent WSR information at the full view and at fine temporal resolution (2.4 ms).

B. Enabling Longitudinal Studies of Wall Shear Dynamics
Although the role of WSS in atherosclerotic plaque induction, development, and rupturing has been a long-standing inquisition, the intricate influence of wall shear on these protracted processes remains incompletely answered [43]. Progressive monitoring of the evolving WSS dynamics with longitudinal changes (e.g., plaque morphology) would advance our understanding of WSS's role in the various phases of atherosclerosis. In contrast to existing WSR assessment tools involving MRI and CFD, WASHI can directly acquire shear information noninvasively on a point-of-care setting for longterm monitoring, therefore facilitating longitudinal studies. Moreover, WASHI is suited for monitoring time-sensitive events such as plaque rupture that are restricted by MRI's long scan time to acquire the patient-specific vessel geometry and boundary conditions for CFD simulation [10]. In fact, WASHI enables direct observation of WSR spatiotemporal dynamics during plaque disruption that could not be performed using existing techniques. These observations may bolster our understanding of the role of WSR in plaque destabilization.

C. Detecting Atherogenic Region Using WASHI
It is widely agreed that arterial segments with decreased shear stress favor the development of atherosclerotic lesions. With this notion in mind, metrics, such as the low shear index [44] and time-averaged WSS [45], were introduced to characterize the magnitude of WSS over a cardiac cycle. These metrics can be readily derived from WASHI by linearly scaling the estimated WSR with the dynamic viscosity of the fluid to calculate the WSS; post-conversion WSS scales for all supplementary videos and results figures (see Figs. 6-9) by multiplying the derived WSR with the dynamic viscosity of the fluid (4.1 mPa·s for BMF and blood) were included. Here, a common metric enables a comparison between future WASHI investigations with those in the literature in predicting plaque development.
More recently, the multidirectionality of WSS was investigated as a potential atherogenic factor to predict the plaque development [46]. The multidirectionality of WSS, attributed to transitions between antegrade and retrograde flows, could be observed when flow recirculation is present [47]. Directional transitions not only occur during the formation and dissipation of flow eddies but also frequently observed at the boundaries of a recirculatory flow zone when it changes in size [48]. Multidirectional WSS due to the latter is apparent in our imaging experiments involving the bifurcating phantoms, at the recirculatory flow zones formed at the carotid bulbs. Over a cardiac cycle, the recirculatory zone evolved in size and flow separation points (marked in Figs. 7 and 8 with a red arrow for separation point and yellow arrow for reattachment point) denoting the boundary of the recirculation zone shifts in position-arterial segments within the path of the separation points experienced multiple transitions in WSS directions. The multidirectionality of WSS can be quantified by deriving the oscillatory shear index [49] or relative residence time [50]. Moreover, the co-visualization of WSR and vector projectiles in WASHI enables these flow separation points to be identified, thus providing a more intuitive representation to explain the multidirectional WSS in these regions.

D. Limitations and Future Work
WSR measurements are highly sensitive to accurate bloodwall boundary identification. Specifically, a 0.2-mm deviation could lead to up to 60% error based on the findings of a previous study using contrast-enhanced ultrasound [25]. While the use of contrast agents would enable accurate vessel wall identification, the injection of microbubbles is suboptimal to be included as part of a clinical routine. Although WASHI has demonstrated its capabilities as a contrast-free solution to map WSR, one technical limitation of our framework lies in the need to manually demarcate the wall boundaries for a given frame of the image sequence to facilitate the automatic wall tracking procedure. The intraoperator variability, as well as pitfalls of the wall tracking algorithm (e.g., out-of-plane motion), should be investigated in the future. As a potential solution to eliminate operator dependency and achieve full automation, the Doppler power can be leveraged to determine the blood-wall boundary. We demonstrated the feasibility of this approach as an alternative approach but noted that wall motion cannot be disregarded to ensure that WSR is measured at the blood-wall interface. In addition, a boundary detection performance using conventional power Doppler may suffer from cyclic changes in flow velocities. In the future, this issue could be addressed by employing advanced signal processing techniques such as eigen-based methods [51] or deep learning approaches [52] that would improve the performance in detecting wall boundaries.
The performance of WASHI, in terms of WSR accuracy and efficiency, could be improved by optimizing the parameters in the acquisition and processing pipeline. A parameter space investigation, similar to those performed in our earlier work [30] to optimize vector flow estimation, would be beneficial in improving WASHI by examining the performance tradeoffs of different parameters in Doppler processing, clutter filtering, wall tracking algorithm, flow field regularization, and derivation of the spatial gradients. More importantly, investigations of WSR estimation bias at near-wall regions due to different clutter filtering approaches would be beneficial in improving WASHI. In this work, we have taken the approach to employ a regularization technique that also interpolates the velocities at the near-wall region. We surmize that adaptive clutter filtering techniques [53] coupled with a velocity estimator tailored specifically for slow-flow regions [28] may enable direct velocity estimation at the near-wall regions, in place of interpolation strategies.
To facilitate future investigations, a new WSR instrumentation strategy would be needed. While WASHI measurements were corroborated against CFD-simulated WSR on a Hele-Shaw flow model, this approach is limited to providing reference values only under steady flow conditions. In fact, CFD simulations may deviate from the actual flow due to mismatches in boundary conditions that are more prevalent in the case of pulsatile flow. The evaluation of WASHI's accuracy becomes challenging in these flow scenarios as there is an inherent lack of ground truth. As an alternative, future investigations could be performed with a flow phantom that is compatible with both ultrasound and optical velocimetry techniques to derive reference WSR values [54]. Such an investigative platform will be beneficial in evaluating the robustness of WASHI.
To efficiently perform in vivo measurements of WSR dynamics, real-time implementation of WASHI will be critical. Our framework can be readily implemented on our advanced open research scanner [55] to provide the WSR information in real time. This research scanner is capable of providing channel RF data and is equipped with a GPU to perform the necessary signal processing in real time. As vector flow imaging has been demonstrated in real time on clinical scanners [56], the hurdle in achieving real-time WASHI would now mainly lie at the wall tracking module. This step would benefit from a GPU implementation to accelerate the processing by leveraging parallel computing [57]. Real-time implementation of WASHI will facilitate large-scale in vivo studies in the future, for example, in detecting the changes in WSR due to flow-mediated dilation among different age groups [58].

VII. CONCLUSION
WSS plays an intricate role in the development, progression, and rupture of atherosclerotic plaques. Routine monitoring of WSS would aid plaque management but is challenged by the absence of accessible tools to noninvasively assess WSR. In this work, WASHI has been presented as an ultrasound-based solution that offers full-view visualization of WSR dynamics without the need of contrast agents, aligning with common clinical practice. It is anticipated that longitudinal in vivo studies using WASHI would bolster our understanding on the causal effect between WSR and the phases of plaque development. This body of knowledge can enhance the current standard in the prevention and treatment of atherosclerosis and, in turn, mitigates the risk of potential cerebrovascular events.