Blood ﬂow patterns estimation in the left ventricle with low-rate 2D and 3D dynamic contrast-enhanced ultrasound

Background and Objective : Left ventricle (LV) dysfunction always occurs at early heart-failure stages, pro- ducing variations in the LV ﬂow patterns. Cardiac diagnostics may therefore beneﬁt from ﬂow-pattern analysis. Several visualization tools have been proposed that require ultrafast ultrasound acquisitions. However, ultrafast ultrasound is not standard in clinical scanners. Meanwhile techniques that can handle low frame rates are still lacking. As a result, the clinical translation of these techniques remains limited, especially for 3D acquisitions where the volume rates are intrinsically low. Methods : To overcome these limitations, we propose a novel technique for the estimation of LV blood velocity and relative-pressure ﬁelds from dynamic contrast-enhanced ultrasound (DCE-US) at low frame rates. Different from other methods, our method is based on the time-delays between time-intensity curves measured at neighbor pixels in the DCE-US loops. Using Navier-Stokes equation, we regularize the obtained velocity ﬁelds and derive relative-pressure estimates. Blood ﬂow patterns were characterized with regard to their vorticity, relative-pressure changes (dp/dt) in the LV outﬂow tract, and viscous energy loss, as these reﬂect the ejection eﬃciency. Results : We evaluated the proposed method on 18 patients (9 responders and 9 non-responders) who un- derwent cardiac resynchronization therapy (CRT). After CRT, the responder group evidenced a signiﬁcant (p < 0.05) increase in vorticity and peak dp/dt, and a non-signiﬁcant decrease in viscous energy loss. No signiﬁcant difference was found in the non-responder group. Relative feature variation before and after CRT evidenced a signiﬁcant difference (p < 0.05) between responders and non-responders for vorticity and peak dp/dt. Finally, the method feasibility is also shown with 3D DCE-US. Conclusions : Using the proposed method, adequate visualization and quantiﬁcation of blood ﬂow patterns are successfully enabled based on low-rate DCE-US of the LV, facilitating the clinical adoption of the method using standard ultrasound scanners. The clinical value of the method in the context of CRT is also shown.


Introduction
Heart-failure is one of the leading causes of death globally. The total percentage of the population with heart-failure is predicted to increase from 2.42% in 2012 to 2.97% in 2030 [1] . In heart-failure, many cases are associated with cardiac dyssynchrony, resulting tion of multidirectional and vortical flow patterns that vary over the cardiac cycle. Analysis of blood flow in the LV is of great value to study the cardiac function. Using computational fluid dynamics (CFD), researchers modelled the normal LV and patient-specific LVs based on magnetic resonance imaging (MRI), simulated blood flow in the LV and then analysed blood flow patterns. Quantitative analysis of these patterns demonstrated difference between healthy subjects and patients suffering from cardiac dysfunctions, revealing the correlation between ventricular function and blood flow patterns in terms of ventricular surgical restoration and hypertrophic cardiomyopathy [4][5][6] .
In order to quantify the blood flow patterns and associate them with ventricular functions, several flow features could be extracted from blood flow patterns. Blood flow patterns comprise a number of flow features such as vorticity, pressure changes and viscous energy loss. In recent years, researchers reported flow features such as velocity, vorticity, rate of pressure change during systole and kinetic energy; these features could be employed to assess the cardiac function by modelling the LV of healthy subjects and patients using CFD simulations [7,8] . Vorticity (the curl of the velocity field) quantifies the circular, swirling motion characterizing a vortex. The formation of a vortex is required for an efficient LV function; the vortex decelerates blood flow during LV filling and redirects the flow for ejection. Typically, a counterclockwise vortex appears in the central region of the LV, which persists from late diastole to the isovolumetric contraction period. The formation of abnormal vortices is associated with dysfunction of the LV. As such, many researchers have investigated vorticity as an indicator for LV condition [9,10] .
Additionally, the maximum rate of the pressure increase during systole (peak dp/dt) was demonstrated as a useful indicator of ventricular performance [11] . During systole, the pressure difference gives rise to flow acceleration towards the aortic valve region for flow ejection.
Moreover, viscous energy loss is also derived from the velocity vector field as an indicator of cardiac function and efficiency. To some extent, the viscous energy loss is related to timing and geometry of the formed vortex [12,13] .
Visualization and quantification of the spatial and temporal distribution of these patterns may provide valuable diagnostic information on the ventricular function [2,14] .
Several techniques have been proposed for the visualization of blood flow patterns. By MRI techniques, researchers proposed direct and indirect ways to visualize blood flow patterns. In the direct way, the latest advances of 4D flow MRI enables timevarying 3D visualization of cardiac blood flow, facilitating the diagnosis of congenital heart disease [15] . In the indirect way, MRI was employed to provide the LV geometry during the cardiac cycle. Reconstruction of the LV geometry was then used for CFD simulation. The CFD simulation enabled the visualization of diastolic blood flow patterns characterized by flow velocity magnitude and Reynolds number [16] . By ultrasound techniques, these techniques are typically based on the speckle-tracking method using ultrafast ultrasound acquisitions. A well-known technique referred to as echo-particle imaging velocimetry (echo-PIV) enables a timeeffective estimation of 2D intraventricular flow patterns [2,[17][18][19][20][21] . The velocity vector field are estimated by taking regional pattern matching between two subsequent frames (speckle tracking) after injection of an ultrasound contrast agent (UCA). In the studies by Kheradvar [2,17,21] , although commercial clinical ultrasound scanners were used for fluid dynamics investigation, a frame rate of about 80 Hz was adopted, compromising the image quality in order to gain speed. For these techniques, the frame rate needs to be sufficiently high to spatially track the fast motion of blood. When this rate is low, as for 3D ultrasound imaging, speckle-tracking techniques break down, as speckle decorrelates significantly between frames [18] .
Alternative flow vector imaging methods are based on Doppler ultrasound without UCA injection [22,23] . In [24,25] , combination of speckle-tracking and color flow imaging (CFI) allows 2D and 3D cardiac vector flow imaging in pediatric cardiology, but high frame rates (or high volume rates) of about 50 Hz are still required.
Although few advanced systems offer 2D ultrafast imaging options, most clinical scanners have limited frame rates, which is a general limitation especially for current 3D imaging solutions. To address this limitation, a novel technique is proposed to estimate blood flow velocity vectors and relative pressures. The method is developed and tested with both 2D and 3D dynamic contrastenhanced ultrasound (DCE-US) at low rate. To this end, we borrow concepts from [26] , where a constant velocity vector field is estimated in the prostate employing low-frame-rate DCE-US; timedelays are estimated by maximizing the cross-correlation between neighboring time-intensity curves (TICs) and the associated velocities are derived by dividing the estimated time-delays by the known inter-pixel distance. Here, we extend this technique to estimate dynamic (time-varying) cardiac flow velocity vectors, and refine the resulting spatiotemporal vector field by a regularization strategy based on Navier-Stokes equation. This approach also enables the estimation of relative pressures. Based on the estimated velocity vector and relative pressure fields, quantitative analysis of these patterns were performed to obtain flow features. As these flow features reflect LV flow efficiency, and CRT aims at improving pumping efficiency, we investigate the changes of these features responding to CRT by comparing these features on 18 patients (nine responders, nine non-responders) before and after CRT. To some extent, this investigation is performed as a preliminary validation of the proposed method.
In this work, using the proposed TIC delay estimation, the obtained results suggest the proposed method to enable adequate visualization and quantification of blood flow patterns based on low-frame-rate 2D and low-volume-rate 3D DCE-US of the LV. The employment of low acquisition rates facilitates the clinical adoption of the method using standard ultrasound scanners. Moreover, many other studies have investigated the relation between cardiac fluid dynamics and cardiac function, especially for those patients who underwent CRT [9,27,28] . By comparison with these study conclusions, our investigation of CRT response with flow features also shows promising results. The clinical value of the presented method concerned with CRT is also shown.
The remainder of this paper is organized as follows. In section 2.1 , the DCE-US data acquisition procedure is described. The algorithmic details of the 2D version of the proposed method are presented in sections 2.2 to 2.7 , and then expanded to the 3D case in section 2.8 . The performance of the proposed method on a clinical dataset is detailed in sections 3.1 to 3.2 . We display how several flow features derived from the estimated velocity and relative pressure fields enable the assessment of patient response to CRT. In section 3.3 , we show the method application to 3D DCE-US data. The results are discussed in section 4 and conclusions derived in section 5 .

Data acquisition
The 2D ultrasound data were acquired at the Catharina Hospital (Eindhoven, the Netherlands). The study was approved by the local ethics committee and all patients provided written informed consent to participate in this study. Transthoracic single beat, 2D DCE-US imaging was performed on 18 heart-failure patients from the apical four-chamber view using left ventricular opacification

ARTICLE IN PRESS
JID: COMM [m5G; November 17, 2020;20:4 ] (LVO) settings with a Philips iE33 ultrasound scanner equipped with an S5-1 sector array transducer (Philips Healthcare, Andover, MA, USA) transmitting with central frequency of 3.5 MHz. The acquisition frame rate was 23 Hz and pixel size was about 0.4 mm in both directions. The scan duration was 60 s following an intravenous injection of a bolus of 10-mL SonoVue® (Bracco, Milan, Italy) contrast agent solution with a concentration of 10 μL/mL, enabling recording the full in-and out-flows in the LV [29] . All the scans were performed during apnea, to minimize translation movements of the heart. In order to avoid bubble disruption, the mechanical index (MI) was set to 0.19. The dynamic range was set to 50 dB to make sure that the full time-intensity curves were obtained without clipping their peak (no saturation effects) and such that the full curves were above the noise floor [30,31] . Imaging was performed before and 3 months after CRT. The criterion to distinguish between CRT responders and non-responders was based on a decrease in the end-systolic volume ≥15% measured by bi-plane echocardiography [32,33] . As a reference, pulsed-wave Doppler through the mitral valve was also performed in all patients. One additional 3D ultrasound dataset was also acquired at the Catharina Hospital. A Philips EPIQ ultrasound scanner equipped with an X5 3D matrix transducer (Philips Healthcare, Andover, MA, USA) transmitting with central frequency of 3.5 MHz was used to perform 3D DCE-US with the apical four-chamber view using left ventricular opacification (LVO) setting on a single patient with anterior myocardial infarction in 2006. The scan duration was 80 s and the volume rate was 3 Hz. The MI was set to 0.29. The acquisition was performed during apnea. The two additional 3D data acquisitions were performed at Umberto I Hospital of Sapienza University of Rome (Rome, Italy) by the same ultrasound scanner with the same scanning protocol, but the volume rate was 8 Hz.

Pre-processing
For the 2D data, first, the LV blood pool was manually delineated frame by frame over the cardiac cycle, achieving a dynamic contour. In order to observe relative pressures across the mitral valve, the segmentation included the left atrium (LA) during diastole, but excluded it during systole, when the valve is closed. Then, the impact of spatially incoherent noise was reduced by convolution of the ultrasound image with a 2D spatial Gaussian filter (standard deviation = 2 pixels) [34] . Moreover, temporal upsampling by a factor 2, achieving a virtual frame rate of 46 Hz, was performed to achieve high-resolution time-delay estimation by performing zero-padding in frequency domain.
For the 3D data at 3 Hz, the impact of noise was mitigated by a 3D spatial Gaussian filter (standard deviation = 1 voxel). An interpolation by a factor 4 along the time sequence was then performed to attain a virtual volume rate of 12 Hz by performing zero-padding in frequency domain. For the 3D data at 8 Hz, an interpolation by a factor 2 was performed after the 3D spatial filtering.

2D velocity vector field estimation
Velocity vectors are estimated by associating time-delay estimates with inter-pixel distance vector measurements within a circular kernel [26] . As shown in Fig. 1 , N points are distributed along the circular kernel, with radius R = 3 mm. TICs at each kernel point are obtained by collecting the evolution of the pixel intensities over time. By maximizing the cross-correlation between each pair of TICs, time-delays are calculated. However, different from [26] in which the blood flow moves in a convective-dispersion way at a relatively constant velocity, here blood flow velocities vary considerably during the cardiac cycle. Maximizing the crosscorrelation between TICs results in average time-delay estimates. As such, we leverage moving time windows (one per cardiac cycle) to collect TIC samples from the same phase of subsequent cardiac cycles, determined by synchronized electrocardiography (ECG) (see Fig. 1 ). In particular, the R peaks in the ECG are detected and used as time reference to determine the start of systole for each cardiac cycle, and the time interval between subsequent R peaks indicates the duration (number of frames) of each cardiac cycle. As a result, the time-windowed DCE-US loops contain similar blood-flow patterns, as shown in Fig. 1 for diastole. Time delays are then estimated by maximizing the cross-correlation between each pair of time-windowed TICs: with E[ ·] being the expectation operator and x the coordinates of the kernel center. In practice, we first collect frames from the same phase of subsequent cardiac cycles using a nine-frame window (centered at current frame), and then apply a circular kernel to collect time-windowed TICs; therefore, C t f refers to all the timewindowed TICs collected at current frame t f . With the time windows moving, dynamic time-delay estimates are obtained.
The distance matrix D represents the inter-pixel distance vectors between all pixel pairs with the size of 2 × N ( N − 1) / 2 . The average velocity vector v ( x , t f ) within the circular kernel can be The circular kernel is automatically moved throughout the whole image frame with 10 pixels interval. These estimation steps are repeated for different positions of the circular kernel over a grid of pixels covering the entire LV. This way, the dynamic (full spatiotemporal) velocity vector field of the LV is constructed.

2D Navier-Stokes regularization
The presence of noise in the ultrasound recording, due to the microbubble kinetics and the image formation process [35,36] , can affect the accuracy of the estimated vector field. Therefore, in order to decrease the impact of noise on the estimated vector field, model-based regularization is implemented. This is based on the Navier-Stokes equation, which is employed to represent the relationship between velocity v and relative pressure P in a viscous fluid. Gao et al proposed a modified Navier-Stokes equation, where the extra force term is replaced by a closeness term, defined as the difference between the regularized velocity vector v and its initial estimates ˆ v for a particular frame [18,37] . The resulting equation is given as where the mass density of blood ρ and the dynamic viscosity of blood μ are assumed equal to 1060 kg/m 3 and 0.0035 Ns/m 2 , respectively, and λ denotes the regularization parameter. The second equation imposing divergence of the velocity equal to zero corresponds to the continuity equation.

ARTICLE IN PRESS
In order to obtain the numerical solution of (4) , we employed a fractional-step projection method that is based on the Helmholtz decomposition theorem [18,[37][38][39] . This method is developed in three steps, which are applied to each frame in an iterative fashion.
In the first step, after omitting the pressure term, the differential equation in (4) is discretized by employment of the explicit Euler scheme at the n th temporal step, The result is given as where v n +1 / 2 (t f ) is obtained as an intermediate velocity vector at the (n + 1 / 2) th temporal step for the current frame t f . In order to construct the explicit discretization, the time interval t should be set to a small value and also satisfy the Courant-Friedrichs-Lewy condition [38] . Note that t is an optimization parameter which should not be confused with the actual time elapsing during the cardiac cycle.
In the second step, we decompose the obtained intermediate to the incompressible flow velocity and a gradient of the scalar function proportional to t∇P n +1 according to the orthogonal decomposition theorem [37] . This decomposition is presented as It is important to notice that the variable P n +1 is a numerical quantity that is functional for the iterative scheme and does not correspond to the actual relative pressure. However, the term P n +1 in (6) is still unknown.
At last, the term ρv n +1 (t f ) in (6) can vanish under the assumption that the velocity vector at the n th temporal step is assumed to be incompressible. Hence, the Poisson equation for the derivation of P n +1 is formulated by taking the divergence of (6) : where a successive overrelaxation method is employed to solve the equation [40] .
At current frame, aforementioned steps are repeated until the following convergence criterion is satisfied: where || · || 2 indicates the L 2 -norm.

2D relative-pressure estimation
With the obtained regularized velocity field v , the relative pressure P can be derived from the solution of where the first and second terms on the right-hand side are related to the flow inertia. The transient inertia representing the local acceleration refers to ρ ∂ v ∂t , which describes the velocity variation over time at a given location; the term ρ( v · ∇ ) v results from the convective inertia representing the convective acceleration, which describes the velocity variation over space at a given (current) frame; the term μ∇ 2 v represents the viscous resistance that typically appears near the chamber wall due to friction (shear stress) of the moving fluid with the wall; the term f e denotes the body forces. The body forces have no effect on the fluid for a closed system in which the buoyancy force and the gravitational force cancel out [41] . The local acceleration can be derived from the velocity difference between subsequent frames, while the convective acceleration and viscous resistance can be calculated using the obtained velocity vectors at the current frame.
Finally, similar to the third step in section 2.4 , the pressure Poisson equation is obtained by taking the divergence of (9) and then solved using the same method explained in section 2.4 , which enables the estimation of relative pressure based on the obtained velocity field.

ARTICLE IN PRESS
Here, the overall degree of vorticity is estimated as the mean vorticity during diastole over the full ventricular area.

Relative-pressure change
The maximum of the relative-pressure time derivative, peak dp/dt, is estimated in the LV outflow tract, in front of the aortic valve during systole phase.

Viscous energy loss
The friction between rapid inflow and ventricular wall generates large viscous energy loss (EL), which can be calculated as [19] : where A represents the LV area. In the practical implementation, we calculated the EL value by firstly obtaining the squared spatial differential of the velocity components at each sample point within the LV, which were then integrated.
For (11) and (12) , the spatial and temporal differences are calculated using central difference approximation. All the analysis was implemented in Matlab (Mathworks, Natick, MA, USA), running on a mobile workstation (Intel Core-i7-4710MQ @ 2.50 GHz).

Statistical analysis
Statistical analysis aimed at evaluating the ability of the estimated flow and pressure features to assess the response to CRT. Eighteen patients undergoing CRT received 2D cardiac DCE-US. Of these patients, nine were later identified as responders and 9 as non-responders. The adopted criterion to determine response was a decrease in the end systolic volume by over 15%. The evaluated features were the estimated vorticity, peak dp/dt and viscous energy loss. These features were compared before and after CRT treatment within each patient group. Moreover, the ability to discriminate between the responder and non-responder group was evaluated by calculating the relative variation of each feature before and after CRT. A Wilcoxon rank-sum test was employed to evaluate the statistical significance of the differences in the feature distributions before and after CRT and between groups. Statistical significance was defined as a single-tailed p -value < 0.05.

3D velocity vector fields
The proposed 2D velocity vector estimation and regularization strategy, described in section 2.3 and 2.4 , is also extended to 3D data. To this end, the 12-point 2D circular kernel is replaced by a 100-point (N = 100) 3D spherical kernel i.e., where C t are the TICs collected from the windowed volume sequence. In order to uniformly distribute points on the spherical kernel, the set of angle coordinates φ and θ is defined based on the spherical Fibonacci mapping method [42] . Mean vector fields are computed in systole and diastole by proper windowing of the TICs based on the synchronized ECG recording. The ECG R peaks are detected and used to determine the transition time between diastole and systole at each cardiac cycle. The transitions between systole and diastole are determined by the time when the smallest LV volume is found, based on the ultrasound data. Once systole and diastole are determined for each cardiac cycle, two sets of TICs can be constructed that represent systole and diastole, respectively. These are used to estimate the mean vector fields for the two phases. The estimated 3D velocity vectors were visualized with color-coded streamlines, where the colors, blue to red, indicate the start and end points of streamlines, respectively. Fig. 2 presents the estimated velocity vectors before and after Navier-Stokes regularization for a CRT-responder patient at baseline and three months after CRT at four phases during one cardiac cycle. By comparing the velocity vectors before and after regularization, the quality of the estimated velocity vectors is clearly improved by imposing physical fluid dynamic conditions throughout regularization, showing more pronounced mitral inflow and vortex flow. Qualitatively, the movement of the blood-flow pattern during diastole and systole could be described as follows. In diastole, the vector field shows blood flow from the LA to the LV, where it forms two vortices at the left and right sides of the mitral inflow (jet). From late diastole to the isovolumetric contraction phase, a typical flow pattern characterized by a counterclockwise vortex appears in the middle of the LV moving towards the apex. This produces a deceleration of the blood flow, which is redirected in preparation to the ejection during systole. In systole, the velocity vectors show flow towards the aortic-valve, where blood in ejected through the aortic valve into the aorta.

2D Visualization
A color map is realized for visualization of the relative pressure ( Fig. 3 ); green represents low pressure and red represents high pressure. Although the distribution of relative pressure is different at baseline and three months after CRT, the way the relative pressure varies and the relationship between relative pressure and blood flow still could be depicted: during early diastole, a pressure gradient can be observed across the LA and LV. This pressure difference between the atrium and ventricle, caused by LV relaxation and deformation, accelerates blood from the LA into the LV. This high flow filling the LV results in a decrease of the pressure gradient between the LA and LV, after which, the initial pressure differences are reversed, resulting in deceleration and redirection of blood flow. During late diastole, intraventricular pressure variations are low, while a vortex forms and moves to the middle of the LV. As expected, lower pressure can be observed for areas that are closer to the center of this vortex, During systole, an increase in the pressure gradient is observed as a result of the LV contraction, with high pressure in the apical region compared to low pressure in the aortic-valve region. This gradient produces blood acceleration towards the aortic valve. With the ejection of blood out of the LV, the high-pressure area moves from the apex towards the aortic valve along the LV wall, while low pressure appears in the apical region. This process results in the formation of a reversed pressure gradient at the end of systole. Estimated velocity vectors before and after Navier-Stokes regularization at four phases during one cardiac cycle for a CRT-responder: top row is the baseline result before regularization; the second row is the baseline result after regularization; the third row and bottom row are the three months after CRT results before and after regularization, respectively. The end-systolic volume before CRT is 97.8 mL, and after CRT is 71.8 mL.
data. The early diastolic wave (E-wave) and atrial contraction wave (A-wave) can be observed with both methods, and the shapes of their velocity profiles are similar, as demonstrated in Fig. 4 a and  4 b. From this, the E/A ratio can be determined, amounting to about 0.7 in this case. Clinically, an E/A ratio < 1 indicates LV impairment. By comparing Fig. 4 c and Fig. 4 a, three months after CRT, the pulsed-wave Doppler measurements show that the velocity magnitude increases, especially in the early diastolic wave, although the E/A ratio remains still < 1. By comparing Fig. 4 b and 4 d, these changes in the estimated velocity vectors across the mitral valve are in agreement with the pulsed-wave Doppler measurements. The correlation coefficient between the pulsed-wave Doppler ve-locity and the regularized velocity estimated with the 3-mm circular kernel is 0.76 (R-squared-value = 0.58) for the baseline data, and 0.73 (R-squared-value = 0.53) for the after-CRT data. However, we observe an underestimation of the velocity magnitudes by the proposed method ( Fig. 4 b and 4 d), which depends on the size of the chosen circular kernel ( Fig. 4 b), displaying a trade-off between the ability to measure high flows (larger kernel), and the spatial resolution of the estimates (small kernel). The estimated velocity magnitude increases when the radius of the kernel increases from 3 mm to 6 mm. Note that a very large kernel (and thus low spatial resolution) leads to significant spatial averaging of the high-flow jets, also resulting in underestimation of the velocity magnitude, 6 Fig. 6. Flow features for CRT non-responders at baseline and 3 months after treatment: (a) vorticity; (b) relative-pressure change (peak dp/dt) and (c) viscous energy loss. as shown by the estimated velocity magnitude decrease when the radius of the kernel increases from 6 mm to 9 mm in Fig. 4 b. Table 1 reports the patient characteristics in terms of Doppler and LV end-systolic volume measurements before and after CRT. Fig. 5 presents the estimated flow and pressure features (vorticity, relative-pressure change and viscous energy loss) before and after CRT treatment in the CRT responder group. The vorticity shows a significant increase after CRT with a p -value of 0.04. A significant increase in the peak dp/dt was also observed after CRT with a p -value of 0.02. The viscous energy loss also shows a decreasing trend after CRT, although the difference was not statistically significant with a p -value of 0.40. In general, this result indicates that after CRT a more efficient ejection occurs in the responders.

Assessment of CRT response
In Fig. 6 , the flow features before and after CRT treatment in the CRT non-responder group are presented. In general, there were no significant differences for vorticity, relative-pressure change and viscous energy loss with p -values of 0.10, 0.29 and 0.69, respectively. Fig. 7 demonstrates the relative variation for these three flow features between the CRT responder and non-responder groups. The difference between these two groups were significant with p -value of 0.002 and 0.0067 for vorticity and relative-pressure change, respectively.

3D velocity vector fields visualization
The ability of the proposed method to estimate 3D flow patterns is illustrated in Fig. 8 . For the 3D data sampled at 3 Hz, during diastole, blood flow fills the chamber from the mitral valve ( Fig. 8 a, top) towards the apex, which is indicated by dense streamlines with blue color around the mitral valve region. Besides this main flow, a vortex in the central region of the LV can also be observed. During systole, a large vortex near the apex redirects blood flow towards the aortic valve ( Fig. 8 b, top-left). Fig. 8 c-8 e and Fig. 8 f-8 h demonstrate the estimated 3D velocity vectors for two datasets sampled at a higher volume rate of 8 Hz. During early diastole, the mitral-valve inflow filling the LV is observed. In Fig. 8 c, we can also observe a small vortex forming with the mitral inflow. In the late diastole phase, a large vortex forms in the LV. The aortic-valve outflow can clearly be observed in systole, especially in Fig. 8 e.
By observing the shape of the vortex, this consists of a series of concentric circles with varying radius. In order to analyze the vortex in space, we extracted the velocity vector patterns plane by plane in the late diastole phase for the two datasets sampled at 8 Hz, as illustrated in Fig. 9 . The concentric circles present different radius and orientations in different planes.

Discussion
The new method proposed in this paper enables estimation of velocity vectors and relative pressures from low-frame-rate DCE-US data. The estimation of velocity vector fields in the LV has been investigated by many researchers based on high frame-rate US [18,24,25,43] , and the vortex pattern has been the focus of several studies [2,10,44] . Also according to our study, blood flow patterns including the formation of vortices provide valuable informa-  . 7. Relative variation between baseline and 3 months after treatment for each flow feature in the responder and non-responder group: (a) vorticity; (b) relative-pressure change (peak dp/dt) and (c) viscous energy loss. Significant differences ( p -value < 0.05) are indicated by an asterisk.

ARTICLE IN PRESS
tion that can possibly aid the diagnosis of ventricular dysfunction. The possibility of performing the analysis at low frame rates facilitates the clinical translation of these methods and enables 3D implementation with lower rate. Although standard clinical DCE-US equipment was employed in several fluid dynamic studies using the echo-PIV technique [2,17,21] , the frame rate still needs to be set as high as possible. Moreover, compared with high-framerate US acquisitions, standard DCE-US acquisitions compromise the frame rate to achieve improved contrast enhancement, which is especially relevant for accurate myocardial perfusion assessment by stress echocardiography [45] . Standard DCE-US also supports clinicians with accurate LV endocardial contouring, which is also robust to poor acoustic windows in difficult-to-image patients [46,47] .
In the context of 2D visualization, the estimated time-varying velocity vector fields are qualitatively in agreement with flow patterns that are expected based on LV physiology. First, LV rapid filling during diastole evidences higher flow through the mitral valve associated with the formation of two vortices on either side. These two vortices move towards the apical region along the ventricular wall. At the end of early diastole, the formation of a main counterclockwise vortex appears close to the apical region. During late diastole, this vortex persists, while a second blood inflow occurs as a result of atrial contraction. During systole, the blood flow is directed towards the aortic valve, where it is finally ejected into the aorta. The vortex formation improves the efficiency of the ventricular function, producing flow with minimal viscous energy loss. Several published clinical studies also show different heart-diseasespecific vortex patterns [48][49][50] . The vortex role varies according the phase in the cardiac cycle [49,51,52] . Our findings show a good agreement with a published CFD simulation, in which the vortex formation and its interaction with the ventricular structures over the full cardiac cycle were fully investigated [53] .
In our study, we estimated and investigated flow patterns in heart-failure patients. Consistent with heart-failure pathophysiology, we observed a vortex that persists from late diastole to the isovolumetric contraction period and is characterized by low vorticity values. After CRT, improved vortex shapes that yield significantly higher vorticities were noticeable in CRT responders. In line with pathophysiological expectations, we also measured a reduced viscous energy loss and significantly increased pressure gradients in those patients who responded to therapy. Moreover, relative variations before and after CRT in estimated vorticity and relativepressure change showed significant difference between CRT responders and non-responders. Many other studies have investigated the relation between cardiac fluid dynamics and cardiac function, especially for those patients who underwent CRT. Nakao et al reported that CRT response in patients with dilated cardiomyopathy can be predicted in the presence of a greater summed maximum vortex flow assessed by vortex flow mapping using cine magnetic resonance imaging [27] . Cimino et al evaluated the char-acteristics of LV flow dynamics 6-months after CRT between responders and non-responders using echo-PIV patterns. More compact and uniform vortex formation in CRT responders and more irregular and nonuniform vortex formation in non-responders were observed, along with a significant change of energy dissipation both in responders and non-responders [28] . Goliash et al applied the echo-PIV technique to investigate the effects of acute interruption and reactivation of CRT on cardiac fluid dynamics [9] . As mentioned in this study, heart-failure patients showed a weaker early diastolic counterclockwise vortex formation with significantly reduced vortex intensity compared to normal hearts. Hong et al also showed that quantitative vortex parameters such as vortex relative strength and vortex pulsation correlation have significant higher values in normal hearts than abnormal hearts using the PIV technique [10] . With CRT aiming at improving heart efficiency, the heart condition of CRT responders after CRT treatment may be assumed to be comparable to normal. As such, these studies support our findings based on the difference in flow patterns between CRT responders and non-responders, and between healthy and heartfailure subjects, indicating that the proposed method enables assessment of the LV function and condition. However, the clinical implications of vorticity and viscous energy loss in the assessment of the cardiac function are still not well understood and are being investigated, as the pathophysiological relevance of cardiac fluid dynamics is still unclear [51] .
We then showed that the proposed method can be extended to 3D, dealing with even lower frame/volume rates. The feasibility of our method in 3D was demonstrated by evaluating the flow fields during systole and diastole in three patient datasets. The initial results are encouraging and consistent with physiological expectations, indicated by the mitral-valve inflow in diastole and a large vortex redirecting blood flow towards the aortic valve in systole. By observation of 3D velocity vector patterns, the mitral inflow is not parallel to a certain 2D imaging plane, which reveals that the estimated velocity vectors in 2D are actually the projection of the true 3D velocity vectors onto the imaging plane. For the 2D case, out-of-plane flow affects the velocity estimates and introduces ambiguity in the estimation of the velocity vector magnitude. By comparison of the 3D velocity vector patterns with the 2D patterns in two separate patient datasets, the 2D vortex is intrinsically the projection of a longer 3D vortex consisting of a series of concentric circles with varying radius in different planes. In addition, this 3D vortex can also be slightly twisted in the sagittal and transversal planes, which may generate an angle between the vortex and the 2D ultrasound scanning plane, resulting in inaccurate vorticity estimates. The published CFD simulations also demonstrates difference between 2D and 3D vortex [54] . Considering that vorticity can represent a key flow feature to quantify CRT response, the spatial vorticity of a 3D vortex can provide more accurate and global information to describe the ventricular efficiency, facilitating the  assessment of CRT response. This evidences the need for 3D imaging.

ARTICLE IN PRESS
There are however several limitations to this study. Segmentation of the LV and LA was performed manually, which is timeconsuming. Automatic segmentation will be investigated in future work [55] to speed up the process, facilitating the clinical implementation of the method. Here, the contour of automatic segmentation should be highly accurate because the viscous energy loss is proportional to the sum of the squared spatial differential of the velocity components within the LV, and related to the area of the LV. As such, inaccuracies in LV segmentation affects the accuracy of the viscous energy loss estimates. Moreover, the contour is used as boundary in the regularization. Since the regularization is sensitive to boundary conditions, the contour of automatic segmentation may also influence the regularization step of the proposed method. Although the segmentation can be improved in the future, delineating the LV chamber benefits from the use of DCE-US and the resulting LV opacification, especially in patients that are difficult to image as the acoustic windows of their LV are poor [46,47,56] .
Comparing the peak velocity estimated by the presented method with pulsed-wave Doppler velocity, taken as the ground truth, a peak velocity of 0.3 m/s was estimated, against around 0.6 m/s estimated by Doppler for the presented baseline data, and  a peak velocity of 0.45 m/s was estimated, against 0.8 m/s calculated by Doppler for the after-CRT data. This underestimation may be ascribed to the low-frame-rate, insufficient to catch the fast dynamics of blood flow. This may especially affect the velocity estimation in the diastolic injection phase, possibly producing underestimation of the peak velocity. Although the employment of larger kernel sizes leads to higher velocity estimates, which approach the Doppler estimates, this also results in spatial blurring of the estimated features. The estimated velocity magnitude would again decrease when the kernel is oversized due to blood flow averaging over space. The circular kernel can ensure the rotational invariance, but it may degrade the flow estimation when flow moves towards

ARTICLE IN PRESS
JID: COMM [m5G; November 17, 2020;20:4 ] a specific direction. Moreover, the regularization may also degrade the peak velocity because it is likely to smooth the flow field based on the Navier-Stokes model. The temporal resolution is mainly determined by the adopted temporal window size, whose lower limit is determined by the acquisition frame/volume rate. Large window sizes enable more accurate delay estimation by providing more data points; this is counterbalanced by lower temporal resolution, producing timedelay estimates that represent a larger temporal fraction of the cardiac cycle. An extreme case occurs in 3D, where due to the low volume rate, the estimated vector fields are average representations of the full systolic and diastolic phases.
A common issue during ultrasound acquisitions, especially in echocardiography, is probe motion due to freehand scanning and myocardial tissue motion, which may introduce artifacts in the extracted TICs because of spatial variations along the time series, resulting in inaccurate velocity estimation. Motion compensation can be considered to mitigate this problem [57] . Even so, out-of-plane motion in 2D cannot be estimated, and 3D acquisition is required for proper motion compensation.
Moreover, in solving the Navier-stokes regularization, blood was assumed to be a Newtonian fluid with constant viscosity during the whole cardiac cycle. However, it has been shown that the non-Newtonian property of blood may be relevant during the deceleration period of the cardiac cycle, when the velocity and velocity gradients approach zero. As a result, different flow patterns, especially for diseased LV, may form when accounting for non-Newtonian blood properties [58] . CFD simulations also indicated that the mitral valve and aortic valve could affect the flow pattern, especially the vortex formation [59] . Therefore, including the mitral and aortic valve in the proposed method might lead to improved estimation performance.
In general, the proposed method has only been tested on a small group of patients focusing on a specific cardiac pathology. In the future, extended validation should be carried out, also considering other pathologies that may be reflected in LV flow and pressure patterns. Once the proposed method is proven, other clinical applications in addition to the assessment of CRT response could be investigated, such as identifying optimal site of CRT implant [60] , optimizing the CRT device setting [28,60] , and diagnosing several specific cardiac diseases such as dilated cardiomyopathy and hypertrophic cardiomyopathy [61] .

Conclusion
Using the proposed TIC delay estimation, windowed in subsequent phases of the cardiac cycle, velocity vectors are successfully measured based on low-frame-rate 2D and low-volume-rate 3D DCE-US of the LV. The employment of low acquisition rates facilitates the clinical adoption of the method using standard ultrasound scanners. Being based on DCE-US, the proposed method also facilitates the investigation of patients showing poor LV acoustic windows, who are otherwise difficult to image in daily clinical practice. Using Navier-Stokes regularization, an improved estimate is obtained, from which relative pressures can also be derived. The obtained results suggest the proposed method to enable adequate visualization and quantification of blood flow patterns. The clinical value of the method in the context of CRT is also shown.