Optimized Time-Resolved Echo Particle Image Velocimetry– Particle Tracking Velocimetry Measurements Elucidate Blood Flow in Patients With Left Ventricular Thrombus

Numerous studies have shown direct correlations between left ventricle (LV) flow characteristics and cardiac function [1–4]. Hence, resolving the LV flow structure might assist in proper diagnosis and treatment of cardiovascular diseases. Phase-encoded magnetic resonance imaging [5] measures three-dimensional (3D) velocity distributions with good spatial resolution (1 mm), but its temporal resolution (100–200 ms) is insufficient for resolving LV flows. Color Doppler [6] has a higher temporal resolution (20–40 ms) for a single velocity component. Echo particle image velocimetry (PIV) or particle tracking velocimetry (PTV) provides time-resolved (10–20 ms) two-dimensional (2D) velocity distributions at a spatial resolution of 3–5 mm. Applications of echo-PIV for cardiovascular diagnostics have expanded over the last decade [7–10]. It has evolved from the well-established optics-based PIV [11], which consists of acquiring image pairs of illuminated flow field sections containing numerous particles. The images are divided into interrogation windows, and the spatially averaged velocity in each window is calculated by cross-correlating the intensity distribution in successive frames. A single velocity vector is assigned to each window. To obtain reliable cross-correlation based statistics, the typical interrogation window should contain at least 5–10 homogeneously distributed particles. Hence, the uncertainty in the PIV data increases in situations where the particle concentration is low or spatially nonuniform [12]. Contrarily, particle tracking velocimetry [11,13] determines the velocity from the displacement of individual particles, irrespective of their size. Hence, the vector assigned to each particle is not sensitive to the concentration or spatial distribution of particles. Yet, it is sensitive to the precise matching of the same particle in consecutive frames. Consequently, application of PTV, which is broad, typically involves elaborate procedures with multiple criteria for matching corresponding traces. It is quite common [13–16] to use PIV data as one of the matching criteria. In the present application, which involves sparse data, PTV is used for refining the PIV measurements.

Echocardiographic particle image velocimetry facilitates flow measurement where there is no optical access, such as the human cardiovascular system. Contrast agents [17] are intravenously injected and used as flow tracers. Sound “scattered” by them produces PIV-like B-mode images. In general, echo-PIV involves four stages: (i) acquisition and/or reconstruction; (ii) image enhancement to maximize the probability of tracer identification; (iii) velocity calculation; and (iv) validation to remove or replace “outliers.” Ultrasound imaging for PIV measurements was first applied to sediment flows [18]. Medical applications first reported in Kim et al. [19] examine laminar pipe flows and compare correlation-based echo-PIV to optical-PIV. The same group subsequently developed a custom ultrasound system with a high temporal resolution (0.6 ms) [20,21]. They used standard PIV to obtain “guess” displacements, followed by PTV based on multiple matching parameters to obtain velocity distributions. Several other applications of correlation based echo-PIV are reported for in vitro [22] and in vivo [3,8,7,23–26] cardiovascular flows.

Most of the algorithm development efforts reported in these studies have been invested in improving stage (iii), namely the vector calculation procedures. The present study aims to address problems that would occur if contrast echo images, where concentrations and image quality vary substantially, are analyzed as part of routine clinical examinations. Such applications require a robust processing tool that integrates stages (ii)–(iv), where the velocity is calculated for systematically varying processing parameters, and the (most likely) correct vector is obtained through an automated iterative optimization and outlier removal procedures. Furthermore, to extract as much data as possible, the PIV results are adopted as one of the multidimensional set of parameters for measuring the motion of individual traces using PTV. The PIV and PTV data are then integrated to obtain the final result. To demonstrate their potential, these procedures are implemented to analyze contrast echo images obtained from patients with left ventricular thrombus (LVT).

Sengupta et al. [27] have introduced echo-PIV to study in vivo LV flows using porcine subjects, followed later by humans. Subsequently, it has been used to study several LV diseases, such as hypertrophic cardiomyopathy (CM) [4], heart failure (HF) [8], acute myocardial infarction [10] and dilated cardiomyopathy [28]. The vortex strength, shape, and location, as well as momentum and kinetic energy [10], have been calculated from the velocity distributions and used as indicators of cardiac function.

Patients with cardiomyopathy have reduced LV function, resulting in abnormal blood flow and stasis, possibly leading to the formation of an LVT. Most CM patients, however, do not develop LVT. In fact, a prior comprehensive study [29] that has assessed 784 patients with systolic dysfunction, found that only 7% (55 patients) had LVT. The latter are at a very high risk of stroke or systemic embolization [30]. Around 795,000 people in the U.S. alone have a stroke each year and 87% of the strokes are ischemic, with cardio-embolic strokes accounting for 15–40% of them [31]. Hence, as a precautionary measure, most CM patients are prescribed anticoagulants, unnecessarily exposing a large fraction of them to a higher risk of bleeding and severely impairing their quality of life. Therefore, differentiating patients at risk of LVT formation from those who are not is of critical clinical importance.

Left ventricular thrombus is caused by cell injury in the myocardium along with local blood stasis, leading to the activation of the coagulation cascade. It is typically observed near the apex of the LV. The formation and growth of blood clots is a function of the coagulation cascade, flow dynamics, and other biochemical factors. These coupled effects have been studied extensively [32] in the past and are beyond the scope of this study. Moreover, early Doppler-based studies [33,34] have shown that abnormal apical blood flow patterns are correlated with the presence of LVT. The only LVT study [7] where echo-PIV measurements have been performed previously also reports a definite correlation between the LV vortex parameters and thrombus formation. It reports that the LV vortex averaged over the entire cardiac cycle has a lower penetration depth for patients with LVT. However, there are no well-validated criteria relating the flow dynamics to LVT formation that can assist in LVT risk stratification among patients. Hence, the present study focuses on the flow dynamics, and as outlined earlier, the effects of the coagulation cascade and other biochemical factors on thrombus formation are beyond its scope.

The optimized processing procedures are implemented to obtain time series of instantaneous and phase-averaged velocity and vorticity distributions over 10–15 cardiac cycles using images acquired from four LVT patients. Following Refs. [7,8], the data are used for calculating the time evolution of the penetration depth and strength of the LV vortex. In two cases, the decay and fragmentation of the LV vortex as it migrates away from the mitral annulus results in low apical velocity throughout the cardiac cycle in spite of a high penetration depth. This observation is based on direct measurements of the time evolution of apical velocity, and estimates based on the LV vortex-induced velocity. These findings suggest that direct examination of the time evolution of the entire LV flow structure facilitated by the presented PIV/PTV procedures might be useful for assessing LVT risks. Due to the small number of patients involved in this preliminary study, the present effort should be followed by comprehensive clinical tests.

The original study involves 14 patients with past or current confirmed LVT. They have undergone functional cardiac CT, routine and contrast-enhanced transthoracic echocardiogram (echo). To limit circadian effects and influence of physical activity, all patients have been scanned between morning and afternoon while they were awake and lying down on their left side. The idea of integrating PIV/PTV to the diagnosis has been introduced in the latter phases of the project. Consequently, usable contrast echo data exist only for four patients labeled A–D (Table 1). The plan is to use the present results as preliminary motivation for a more comprehensive clinical study. Patients whose LVT has been cured in the past are labeled as “cured” (LVT-CU), while those that have not shown any improvement despite prolonged anticoagulation are labeled as “chronic” (LVT-CH). Patients who were scanned for this study shortly after LVT detection and are responding to anticoagulation therapy are labeled as “recovering” (LVT-R). This study is approved by the Johns Hopkins Medicine Institutional Review Board.

Contrast echocardiography is performed using the GE Vivid-9E system with 1.5–4.0 MHz probes. The mechanical index, a measure of the bio-effects of the ultrasound beam, varies between 0.12 and 0.40 for this study. It is defined as the ratio of the peak negative pressure (MPa) to the square root of the center frequency (MHz) of the ultrasound wave. Contrast enhancement is produced by bolus intravenous injection (antecubital vein) of 0.5–1.0 ml of contrast agent OPTISON, followed by 4–5 ml of normal saline to flush any remaining contrast in the blood vessels. Depending on the quality of blood circulation in the patient, it takes about 30 s for the contrast to reach a peak concentration in the LV. At this phase (Fig. 1(a)), the echo images are saturated with microbubbles, enabling clinicians to delineate endocardial borders and obtain qualitative information about regional LV wall motion and overall cardiac function. However, due to limited local intensity variations, this concentration is too high for detection of individual microbubbles and velocity calculation. About 30 s later, as the concentration begins to decrease, traces of individual bubbles become discernable and yet distributed over the entire LV (Fig. 1(b)), making it optimal for velocity measurements. As the concentration decreases further in the next 60–120 s (Fig. 1(c)), individual bubbles can still be traced, but their number is too low for mapping the overall flow. Based on examination of numerous datasets, echo-images containing 100–200 microbubble traces in the LV appear to be optimum to process using the integrated echo-PIV algorithms.

In the present study, the velocity measurements are based on time-resolved, apical four-chamber (A4C) views recorded at the optimal concentration for ten or more cardiac cycles. The A4C view is a standard orientation recommended [31,35,36] for the assessment of LV function and wall motion and is part of the imaging protocol in almost all echocardiographic examinations. It is consistently reproducible for a given patient. All the present analysis is based on the A4C view. However, if needed, the same procedure may be easily extended to study the flow in other planes, and multiple views could be combined into a complete picture describing the LV flow or flow in other organs.

A frame rate of 60–80 fps is achieved with scanning depths of up to 13 cm and sector widths of up to 60 deg. The scan line density is maximized at first. The probe tilt angle, sector width, and depth are then adjusted to cover the entire LV, thereby prescribing the frame rate. All images are acquired in uncompressed DICOM mode at full spatial and temporal resolution. The cardiac phase is determined for each image using the available electrocardiogram signal. The R-wave peak from the Q-R-S complex is used as the starting reference for each cycle, and for calculating the mean cycle duration. Although the durations vary among patients, the variations within the 10–15 cycles recorded for each patient are smaller than the time delay between frames. Hence, the mean cycle duration is used for determining the phase of each frame. Depending on the heart rate (HR) and acquisition rate; 82, 45, 65, and 52 frames per cardiac cycle are recorded for patients A, B, C, and D, respectively. The processed data consist of time series of instantaneous velocity maps for 10–15 cycles. All the velocity maps at the same phase of the cardiac cycle, i.e., 10–15 of them, are then ensemble-averaged to produce a phase-averaged time series representing the cyclic evolution of the flow structure. Such averaging, while not essential, reduces the measurement uncertainty and compensates for the sparse echo-PIV/PTV data.

The objective of image enhancement is to separate microbubble traces from other objects such as muscles, bones, and tissue. Echo-PIV images introduce several challenges, as illustrated in Figs. 2(a) and 2(b): (1) Nonuniform spatial and temporal seeding resulting from inherent inability to homogenize the microbubble concentration in the LV. (2) Varying image sizes resulting from variations in bubble sizes and their agglomeration. (3) Images of moving tissue boundaries, such as the septal or lateral walls, are often not significantly different from those of microbubble agglomerates. (4) Spatially and temporally varying background noise. For example, the papillary muscle appears intermittently in the center of the LV in the A4C view. As discussed below, addressing these issues requires a flexible enhancement approach with parameters that vary across the image and between frames.

As a first step, a static mask is manually drawn along the approximate LV boundaries. For the manual masking, pixels with a time-averaged intensity above the mean intensity, and a standard deviation falling below a prescribed threshold (∼10%), are not considered. Such bright and slowly varying pixel intensities are typically associated with cardiac walls, tissues, etc. The resulting mask is subsequently used for the entire raw image sequence. The phase-averaged intensity is then subtracted from each image mitigating the moving boundary issue. Two-dimensional segmentation is performed to create a list of detected objects, which includes also noise. To identify the microbubbles among them, we take advantage of the high spatial intensity gradients along the boundaries of microbubbles. The distribution of intensity gradient is calculated (Fig. 2(d)), and only objects bounded by gradients exceeding a prescribed threshold level are retained. The mean of the highest 10% gradient values in an area surrounding the object is used to impose the criteria. The resulting image is presented in Fig. 2(e). The procedure is repeated for several threshold levels, resulting in a series of images, each with its own segmented list containing different number and size of candidate bubble traces. They are used later as part of the optimization process of vector calculation.

Here, R is the full dynamic range of the image (255 for an 8-bit image), α is the fraction of pixels satisfying r ≥ r', and P(ω) is the fraction of pixels with intensity ω within the segmented trace. The value of α is prescribed by the user, and for the echo-PIV data involving only pixels belonging to the segmented traces, α∼0.9. In typical optical PIV, the integration is performed over square areas surrounding the pixel of interest, where the particles occupy a small fraction of the image. Hence, typical values for α are smaller than 0.1.

The velocity calculation is initiated using PIV, followed by optimization, and culminated with PTV refinement. For convenience, an in-house cross-correlation-based code gauCorr [37] is used to perform standard PIV calculations on the enhanced images (commercial codes could also be used). Due to the sparse and nonuniform seeding density, it is difficult to preselect an appropriate interrogation window size with a sufficient number of traces. Furthermore, the spatial and temporal variations in velocity gradients require flexibility in selecting the proper: (i) correlation threshold, (ii) allowed variation between neighboring vectors, and (iii) window size. Hence, a parameter space comprising of the image enhancement (gradient threshold) and the three vector calculation variables is created. The velocity distributions are computed for every parametric combination, resulting in multiple vectors for each grid point. Here, variations in window size are accompanied by corresponding changes to the overlap between windows, resulting in a fixed grid spacing. The horizontal (x, lateral-septal) and vertical (y, mitral-apical) velocity components, u and v, respectively, are a function of location (x, y), time (t), and selected parameters (p-space). This initial velocity field is denoted as (u_raw, v_raw) (x, y, t, p).

Figure 3 is a flowchart illustrating how the most likely correct vector is selected out of the “p” values at each (x, y, t). The initial step focuses on removal of outliers. The process is based on an assumption that if multiple parameters give values of (u, v) falling within a narrow range (to be determined), the correct vector is one of them. For a set of “k” velocity vectors calculated at a single point, {(u_j, v_j)$}j=1k$⁠, the velocity magnitudes and angles are denoted as $m={(uj2+vj2)1/2}j=1k$ and θ={atan(v_j/u_j)$}j=1k$⁠, respectively. The range of values of m and θ at this point, denoted as $Δθ(x,y,t)$ and $Δm(x,y,t)$⁠, respectively, are used as measures of data quality there. Points where $Δθ$ and $Δm$ fall below selected threshold levels of $Δθ*$ and $Δm*$⁠, respectively are considered “good points.” The selected thresholds define grid points with the lowest 35% of $Δθ$ and $Δm$ in the entire dataset as good. For the present data, this percentile is a compromise between stringent quality requirements, and the desire to retain sufficient data for outlier removal. The spatially averaged values of $Δθ$ and $Δm$ denoted by $Δ̃θ(t)$ and $Δ̃mt$⁠, respectively, provide a measure of quality for the entire map. Vector maps with the lowest 15% $ofΔ̃θt$ and $Δ̃mt$ are considered as “good maps.” Corresponding values are used as map thresholds, $Δ̃θ*$ and $Δ̃m*$⁠, respectively. The outlier removal process systematically removes vectors from bad grid points. For a bad point $xi,yj,tk$ to be “eligible” for outlier removal, the following conditions must be satisfied:

• (i)

  Global condition: it must either belong to a good map or be adjacent in time to a good map; i.e., $Δ̃θ(tr)≤Δ̃θ*$ and $Δ̃m(tr)≤Δ̃m*$ for r = k or k − 1 or k + 1.

• (ii)

  Local condition: it must have at least $n*=$10 good neighbors out of the total 26 points possible in space and time; i.e., $Δθxp,yq,tr≤Δθ*$ and $Δmxp,yq,tr≤Δm*$ for at least ten out of $p,q,r=i−1:i+1,j−1:j+1,k−1:k+1,$ not including the point itself.

The local condition ensures that the outcome of the outlier removal is not decided entirely by bad neighbors. Selecting $n*$=10 assures that at least two maps contribute good grid points to the mix, thereby increasing the likelihood that the outcome is spatiotemporally smooth. The global condition ensures that local anomalies, such as a small group of good grid points in an otherwise bad map, do not affect the outcome. Moreover, such a requirement also prescribes the order in which the outlier removal is performed.

For outlier removal (only), each of the neighboring vectors is averaged over p to obtain a single value at each point. Following the universal outlier detection procedure in Ref. [40], the median neighboring vector (u_m, v_m) and the median of the deviation of each neighbor from (u_m, v_m), denoted as (r_u, r_v), are calculated. Of the p vectors at a bad point, only those that deviate from (u_m, v_m) by less than 2r_u or 2r_v are retained. If there are none, all the vectors are removed (set to NaN). Consequently, $Δθ$ and $Δm$ decrease along with their corresponding map-averaged values. This process is performed iteratively as long as it keeps on removing outliers from any of the grid points in the entire dataset. The iterations systematically increase the number of good points and maps for a specific selection of $n*$⁠, $Δ̃θ*$⁠, $Δ̃m*$⁠,$Δθ*$⁠, and $Δm*$⁠, resulting in the initial set of optimal values (u_opt, v_opt)¹(x, y, t, p). Subsequently, the threshold values could be refined/updated and the entire procedure repeated.

In all cases, the data quality varies across the cardiac cycle. For example, phase-averaged distributions of $Δ̃θ$ and $Δ̃m$⁠, denoted as $Δ̃¯θ$ and $Δ̃¯m$⁠, respectively, for one of the cases are shown in Figs. 4(a) and 4(b). It includes values corresponding to the initial raw data (u_raw,v_raw) as well as those resulting from the outlier removal process for thresholds defining 5%, 15% and 25% of the entire data as good. The corresponding threshold levels are (⁠$Δ̃θ*$⁠,$Δ̃m*$⁠) = (43 deg, 1.9 pixels), (49 deg, 2.3 pixels), and (54 deg, 2.7 pixels), respectively. A comparison to the phase-averaged EKG of the same patient (Fig. 4(e)) suggests that the highest $Δ̃¯θ$ and $Δ̃¯m$ for (u_raw,v_raw) occur during the LV ejection/inflow phase, followed by the atrial systole phase. The least variability occurs during diastasis. After outlier removal, the variability is substantially reduced for all threshold levels. The strictest threshold (43 deg, 1.9 pixels) does not classify enough maps as good, resulting in several maps not satisfying the global condition even after iterations. Hence, they are not modified during the outlier removal process, leaving higher values of $Δ̃¯θ$ and $Δ̃¯m$ in parts of the cycle. The other two thresholds result in the same distributions, namely relaxing the level beyond (49 deg, 2.3 pixels) has minimal effect on (u_opt, v_opt)¹(x, y, t, p). Hence, all the subsequent analysis has been performed based on these thresholds. The same type of evaluation has been performed for local goodness thresholds, as summarized in Figs. 4(c) and 4(d), where $Δ¯θ$ and $Δ¯m$ are values phase averaged for all the grid points. In this case, in addition to the raw data, results are shown for thresholds corresponding to 25%, 35%, and 45% of the data. The corresponding threshold levels are (⁠$Δθ*$⁠,$Δm*$⁠) = (14 deg, 0.7 pixels), (21 deg, 1.0 pixels), and (33 deg, 1.5 pixels), respectively. The improvements achieved by outlier removal are evident for the local values as well. Furthermore, selecting $Δθ*$ = 21 deg and $Δm*$= 1.0 pixels, corresponding to 35% of the raw data, appears to be a reasonable compromise. The effect of varying n* from 9 to 11 has been found to be insignificant (results not shown). The overall improvement in data quality is highlighted by the distributions of $Δθ$ and $Δm$ for all the points in the dataset shown in Figs. 4(f) and 4(g). The small fraction of points remaining with high values of $Δθ$ and $Δm$ is associated with sites where either one of the conditions for outlier removal is not satisfied. It should be noted that the specific choice of $Δθ*$⁠, $Δm*$⁠,$Δ̃θ*$⁠, and $Δ̃m*$ that leads to such improvements in the data quality may differ for other datasets. However, their selection can be made by following the process outlined above, i.e., by ensuring that the thresholds are strict, yet allow a sufficient number of maps/points to be classified as good.

The next step involves smoothing of the data. In each point, (u_opt, v_opt)¹(x, y, t, p) is averaged over p. Then, robust gridded smoothing involving all (x, y, t) is performed following the procedures described by Garcia [41] to obtain (u_ref, v_ref)ⁱ⁼¹(x, y, t). This algorithm is based on a penalized least-squares-fit that minimizes the sum of the residuals squared and a penalty representing the roughness of the fitted spline. This tool has been widely used for post-processing of standard PIV data [42–44]. Subsequently, (u_raw,v_raw) is compared to (u_ref, v_ref)ⁱ⁼¹ for each grid point, and only vectors falling within a prescribed deviation, [(u−u_ref)²+ (v − v_ref)²]<2 pixel², are retained to obtain (u_opt, v_opt)ⁱ⁼². The smoothing and comparison steps are repeated until the deviation between successive iterations [(uⁱ_ref − uⁱ⁻¹_ref)²+(vⁱ_ref − vⁱ⁻¹_ref)²] averaged over all points falls below 0.02 pixel². Then, vectors in (u_opt, v_opt)ⁱ that are closest to (u_ref, v_ref)ⁱ are selected as the final PIV velocity field (u_piv,v_piv)(x, y, t). Locations with NaN remain as such, although they could be replaced with the reference values for visualization and guiding the particle tracking phase.

Although (u_piv, v_piv) is optimized, a further refinement can be achieved by tracking the motion of individual particles using PTV. As detailed in prior studies [15,45], when the seeding density is sparse or nonuniform, and when the local flow involves high velocity gradients, the standard PIV correlation-based procedure is inadequate. In such cases, the PIV results along with a series of other criteria (details follow) could be used as a basis for PTV measurements. This approach has been referred to as “super-resolution PIV” [15] in densely seeded flows. Accounting for the specific location of each particle improves the accuracy of the velocity gradients and maximizes the amount of data that can be extracted from the available contrast echo images. Hence, subsequent PTV steps are aimed at refining (u_piv, v_piv) in some areas.

where (G) denotes the threshold for Δ prescribed for every dataset. The candidate with the lowest score is assumed to be the matched trace and included in the unstructured grid (u_ptv, v_ptv), provided it satisfies global constraints: ΔA < 20%, ΔP < 20%, and C > 70%. Traces that are not matched are not considered in the subsequent analysis. Typically, about 20% of the list of segmented particles is successfully matched, with higher fractions found during the diastasis phase. The choice of matching parameters is based on prior experience in optical tracking [13,46–48]. They are particularly suitable for sparse images, where the likelihood of mismatching traces is lower than that occurring in densely seeded optical data. Other parameters that could be incorporated, if needed, include constraints on the velocity magnitude, data continuity, following traces over time, comparison between neighboring matched traces, etc. Extensive manual inspection of the present tracking procedures has confirmed that the above parameters are adequate for the current application. The choice of the numerical thresholds for ΔA, ΔP and C is expected to vary among datasets, due to differences in seeding, frame rate, image quality, etc. As a guiding rule, one may start with a set of strict thresholds (e.g., ΔA < 10%, ΔP < 10% and C > 80%) and relax them incrementally, until a compromise between the fraction of matched particles and accuracy is reached. The final choices should be tested manually over a large sample size of particle pairs. Additional or alternative matching parameters may also be utilized in other situations, as discussed earlier.

Integration and interpolation of (u_piv, v_piv) and (u_ptv, v_ptv) into the final regularly structured vector map is performed using a first-order single value decomposition scheme that has been used successfully before [13]. The analysis is performed over subareas, each containing 3 × 3 neighboring PIV grid points along with the (u_ptv, v_ptv) data in them. In regions without PTV points, the PIV data are retained, and in regions with PTV points, the PIV data are refined. In this process, the weight given to (u_ptv, v_ptv) is doubled when the score given to matched traces falls within the lowest 30%.

Figure 5 shows the distributions of Δ, ΔA, ΔP, C, and S for all the successfully matched PTV candidates for the entire patient B dataset. Each shows the corresponding threshold level. As is evident, Δ and C have clear peaks far from the threshold levels. Conversely, ΔA and ΔP have broad distributions, representing effects of out-of-plane motions and agglomeration, as discussed in Sec. 2.3. Since the true velocity is unknown a priori, the deviation between the velocity obtained from the PIV and PTV could serve as a posteriori assessment of the uncertainty, provided they are independent of each other. A similar approach has been used for assessing the image matching disparity vector distributions in standard PIV [49].

The present values of (u_ptv, v_ptv) are dependent on (u_piv, v_piv) in two steps. First, the location of the search window used to populate the candidate list is derived from (u_piv, v_piv). Considering that the smallest bubble image size is around 2 × 2 pixels, and that the corresponding search window size would be twice the bounding box dimensions, the minimum search area would be 4 × 4 pixels. Hence, the maximum allowed deviation from the PIV result would be 4 pixels, which corresponds to about 40% of the average pixel displacement in the present measurements. Therefore, it represents a relaxed restriction, justifying an assumption that the search area has a weak dependence of (u_piv, v_piv) on (u_ptv, v_ptv). Second, the values of Δ and the associated threshold are inherently dependent on the PIV results. However, the histogram of Δ (Fig. 5(a)) shows a clear peak of Δ = 0.85 pixel and a threshold which is three times the peak value. Hence, the threshold has very limited effect on the selection of particles included in the list of candidates, and the allowed difference between displacements. Similar trends are observed for the other datasets as well, with the Δ peaks ranging between 0.55 and 0.95 pixels, and the threshold is 3 pixels for all of them. Therefore, the most probable value of Δ could be used as a reasonable estimate for the uncertainty in velocity. They correspond to an uncertainty of ±1.4, ±1.8, ±1.6, and ±1.7 cm/s in the instantaneous velocity for patients A, B, C, and D, respectively. Assuming that these uncertainties are caused by random errors, the corresponding uncertainties in the phase-averaged values are $n$ smaller, i.e., ±0.40, ±0.60, ±0.56 and ±0.59 cm/s. It should be noted that the final result is actually a weighted average of the two methods for calculating the velocity, suggesting that the actual uncertainty is even smaller.

A sample instantaneous velocity and vorticity (ω_z = dv/dx − du/dy) distribution overlaid on the original echocardiogram is shown in Fig. 6. Table 1 provides the LVT status, ejection fraction (EF), the Doppler-derived E-wave propagation index (EPI) [50], and the vortex formation time (VFT) [2]. EPI is the ratio between the time integral of velocity at the mitral annulus during the E-wave and the length of the LV. It has been shown [50] as an indicator of LVT formation. Cases with EPI > 1.5 are classified as normal, i.e., the flow has enough momentum to penetrate to the LV apex; EPI∼1 represents borderline apical flow and EPI < 1 is considered high risk for LVT formation. VFT [2] is the ratio of the same velocity integral to the opening diameter of the mitral valve during the E-wave. A value of ∼4 is considered as optimum for vortex formation, and lower values have been correlated with dilated cardiomyopathy patients [2].

Figure 7 shows the distributions of the maximum phase-averaged velocity magnitude at every point. They are overlaid over sample echo images, LVT status, and kinesis levels. Potential implications of these findings in, e.g., assessing the flow-related risk of future clot formation are discussed later. LV wall segments are classified based on the degree of myocardial thickening during systole. Segments that contract normally are labeled normokinetic, while those that do not thicken at all are labeled akinetic. Segments that appear to have subnormal thickening are labeled hypokinetic. To the best of our knowledge, the threshold level of velocity required for effective washing is unknown. Yet, as is evident from Fig. 7, for patients A, B, and C, the apical maximum velocity across all cardiac phases is lower than 5 cm/s in the region of past or current clots and akinesis. For patient D, the maximum velocity is higher than 12 cm/s. In other regions with Hypokinesis or Normokinesis, the velocity ranges between 12 and 25 cm/s.

Figures 8 and 9 compare the distributions of velocity (⁠$u¯2$+$v¯2$⁠)^1/2 and vorticity (⁠$ω¯z$⁠) magnitudes for the same ten selected phases (labeled 0–9), where the overbar denotes phase averaging. A generic EKG signal is provided as a reference. Supplemental videos 1a-b, 2a-b, 3a-b, and 4a-b which are available under the “Supplemental Data” tab for this paper on the ASME Digital collection provide these distributions for all the recorded phases in patients A, B, C, and D, respectively. The parameters of the LV vortex such as maximum vorticity magnitude (⁠$ω¯zmax$⁠), strength, and penetration (location) have been previously used as indicators of cardiac cycle efficiency [8,26] and LVT status [7]. Hence, Fig. 10 tracks the evolution of these parameters for each of the LV vortices along with patient-specific EKG signals for all the cardiac phases. As indicated at the top of Fig. 10, each column represents a different patient.

The swirling strength [51] criterion (λ_ci) is one of the popular parameters used as a means of differentiating between vortices and broadly distributed vorticity in shear flows. For a 2D velocity distribution, where only one vorticity component is known, it refers to the imaginary part of the complex eigenvalue of the velocity gradient tensor. Typically, a distinct vortex exists when the magnitude of λ_ci exceeds a flow-dependent threshold level. In this study, only regions with λ_ci > 1 s⁻¹, corresponding to 10% of the maximum level are considered. The boundary of each vortex is defined by points where $ω¯z=0.5ω¯zmax$⁠. Sample lines defining the boundaries of two vortices are provided in Fig. 9 for patient B (phase 8). Only vortices with an area (A_v) larger than 2% of the LV end diastolic area (A_LV) are considered. To differentiate between multiple vortices at a given phase, each is colored and marked differently (the same marker is used in Figs. 9 and 10). The values of $ω¯zmax$ are normalized by HR. The vortex strength (Γ) is calculated by integrating vorticity over A_v for each structure, normalized as $Γ̂$= Γ/(HR*A_LV). In cases with multiple structures, the combined values of $Γ̂$ are shown in solid black lines. The vortex penetration depth (Y_v/Y_LV) is the distance of the vorticity peak from the mitral annulus normalized by the LV length [26].

As shown by a recent numerical study [52], the near-wall residence time and thrombin concentration decrease with increasing vortex-induced velocity. This velocity is a function of both the vortex strength and its distance from the wall at each phase. For a two-dimensional point vortex, the velocity magnitude is $Γ/2πr$⁠, where $r$ is the distance from the vortex center. Hence, a simple way of estimating the apical flow induced by all the (n) vortices present in each phase, based on the previously introduced parameters, is to calculate the evolution of $v^ind=∑i=1nΓ^i/[2π(1−YVi/YLV)]$⁠. Here, Σ indicates summation over n, the subscript i refers to a specific vortex, and $1−YVi/YLV$ is the corresponding dimensionless parameter characterizing the distance between the vortex and the apex. The profiles of $v̂ind$ are also presented in Fig. 10.

The 2D estimates of $v^ind$ ignore viscous effects and the inherent three-dimensionality of cardiac flows. Alternative formulations that account for viscous effects can be found in Refs. [53] and [54]. However, the present data can also be used for characterizing the apical flow directly. For example, one could average the velocity component tangent to the wall over a 1 × 1 cm (nine vectors) area located near the apex, $vapx$⁠. The evolution of this parameter nondimensionalized as $|v̂apx|=|vapx|/(HR*ALV/YLV)$ is also provided in Fig. 10. It appears that the trends of $|v̂apx|$ agree with those of $v̂ind$ as far as difference in magnitude between the present patients and some of the variations along the cycle are concerned, e.g., for patient D. The two velocity parameters also have similar magnitudes. However, other trends in cyclical variations differ, e.g., those of patients A and B, including the unrealistic values of $v̂ind$ for patient B during phases 0 and 9 caused by $YVi/YLV$ ∼ 1 for one of the vortex fragments. Hence, in spite of its simplicity, the 2D inviscid vortex model appears to provide (an imperfect but) reasonable estimate for the apical velocity magnitude and cyclical variations. The extent of washing might also be affected by the duration of exposure to a certain velocity magnitude. Hence, Fig. 11 presents a histogram of $v̂apx$ versus the corresponding exposure time during a cardiac cycle for the four patients.

As is evident (Figs. 8–10), the entire flow structure varies significantly with phase and among patients. Even the peak inflow and outflow timings vary. The initiation of inflow, characterized by a sharp rise in $ω¯zmax$⁠, appears during the atrial systole (A-wave) for patient B and during early diastole (E-wave) for patients C and D. The location, shape, timing, size, and even the number of structures associated with the LV vortex vary. For patient B and C, the migrating vortex occupies a fraction of the LV, while a persistent large vortex dominates the LV for patient D. These variations affect the distributions of velocity magnitude, including the direction of horizontal velocity during the inflow and outflow phases. Once formed, the vortex propagates toward the apex (Y_v/Y_LV increases) while growing in strength (⁠$Γ̂$⁠), reaching maximum values at the isovolumetric contraction (IC) phase (9-0). These trends are consistent with those reported for healthy and HF patients with a preserved EF [8]. Subsequently, $Γ̂$ declines to its minimum value after systole before the next inflow phase. The change in circulation (Δ$Γ̂$⁠) between IC and the initial inflow for patients B-D is also indicated in Fig. 10. Consistent with the reported trends [8], the present magnitudes of $ΔΓ̂$ = 2.0, 2.1, and 0.8 increase with the corresponding EF = 20%, 30%, and 15% for patients B, C, and D, respectively.

It should be noted, however, that the trends of EF and Δ$Γ̂$ alone do not show a consistent difference between patients B and D, who are recovering from a clot, and patient C, who has a chronic apical clot (Fig. 5). Indeed, EFs of 36.4%±11.9 and 33.2%±9.5 have been reported [7] for non-LVT and LVT cases, respectively, rendering EF-based LVT risk diagnosis unfeasible. However, in accordance with the low inflow speed (Fig. 8), the vortex penetration for patient C (Y_v/Y_LV < 0.3) is significantly lower than those of patient B and D (Y_v/Y_LV > 0.5) throughout the cardiac cycle. Between myocardial infarction patients with and without LVT, the study [7] shows that when averaged over the entire cycle, Y_v/Y_LV < 0.45 is a significant independent condition for the presence of LVT. The present patient C clearly does not satisfy this condition, but patients B and D do. Patient A, who has recovered from LVT before the present study, also satisfies this condition. However, even though the LVT statuses for all the present patients are consistent with the Y_v/Y_LV-based condition, the apical flows in patients A and B are still poor in the previous clot location throughout the cycle (Figs. 7 and 8). Hence, having a high Y_v/Y_LV might not be sufficient for effective washing and full recovery of the wall. This issue is discussed in Sec. 4 for each of the patients. It considers their history, LVT status, EF, EPI, and VFT along with trends observed in echo-PIV results over the entire cycle.

Patient A had nonischemic cardiomyopathy and a clot was detected two years prior. The patient had since been treated with anticoagulation continuously, resulting in resolution of the LVT. The patient passed away seven months later from HF complications. The low velocity over the LV is consistent with the poor VFT and EF. The acquired echo-PIV images only cover the septal half of the LV due to its very large size thereby capturing the LV vortex only during diastole (Figs. 9 and 10, phases 2–9, Supplemental video 1a-b which is available under the “Supplemental Data” tab for this paper on the ASME Digital collection). The restricted domain still enables the evaluation of the flow characteristics during phases 2-9 that may be sufficient to comment on apical washing. However, complete views would be preferable. Although the clot has resolved, the vortex strength (⁠$Γ̂$⁠) and apical washing (⁠$v̂ind$ and $v̂apx$⁠) are very low throughout the cycle. As evident from the histogram of $v̂apx$ in Fig. 11, the apical velocity is extremely low (<3 cm/s) throughout the cardiac cycle in spite of the high vortex penetration (Y_v/Y_LV). Hence, even though the LVT has resolved, the analysis indicates high risk for LVT relapse without anticoagulation, as indicated also by the very low EPI.

Patient B had alcoholic nonischemic cardiomyopathy, and an LVT was detected near the apical side of the septum, following which anticoagulation was started. Eighteen days later, at the time of this study, the LVT was no longer present. The LV contains multiple vortices for a significant fraction of the cycle (Figs. 9 and 10, Supplemental video 2a-b which is available under the “Supplemental Data” tab for this paper on the ASME Digital collection). As indicated by the blue line (marked with ⋄) in Fig. 10, when a new vortex rolls up near the mitral annulus (phase 7), a remnant of an older one (gray line, marked with ○) is located at higher penetration (Y_v/Y_LV). Both vortices migrate toward the apex and eventually merge (green (marked with ▸) and gray (marked with ○), phases 2–3) while shedding a short-lived small structure (yellow, phases 3–4). A new vortex rolls up shortly afterward (phase 7). The pertinent observation from the evolution of these multiple vortices is that their strengths are high when they are located near the mitral annulus, but by the time they reach the apex, they decay by more than 50%. Consequently, the washing (⁠$v̂ind$ and $|v̂apx|$⁠) in the area corresponding to the past clot is lower than other parts of the LV. Yet, this induced velocity is never zero, as evident from Fig. 11, suggesting that the apex is being washed persistently with ∼4 cm/s of flow. Since the minimum velocity required for effective apical washing is unknown, it is unclear whether the washing is adequate for preventing thrombus relapse should the anticoagulation be discontinued, as indicated also by the EPI.

Patient C had ischemic cardiomyopathy and suffered multiple heart attacks in the past. An LVT was detected three years prior and in spite of being on anticoagulants since, all subsequent tests, including the present study have shown a persistent clot, with an EF = 30% higher than the other three patients. The flow analysis provides a consistent picture as shown in Figs. 7–11 and Supplemental video 3a-b which is available under the “Supplemental Data” tab for this paper on the ASME Digital collection. The vortex location and most of the LV flow is confined to the mitral annulus throughout the cycle. The vortex remnants that do penetrate into the LV are very weak and ineffective in apical washing (⁠$v̂ind$ and $v̂apx$⁠). As evident from Fig. 11, the apical washing is extremely low (<2 cm/s) throughout the cycle, indicating very high risk of LVT and providing a likely explanation as to why the LVT did not resolve despite adequate anticoagulation and reasonably fair EF. The high change in vortex strength (Δ$Γ̂$⁠) compares well with the EF.

Patient D had a history of hypertension and hyperlipidemia who presented with breathlessness. Found to have ischemic cardiomyopathy, reduced EF of 15% and severe stenosis of the left anterior descending artery, the patient was treated with percutaneous coronary intervention. The clinical judgment was that the persistent low EF and sluggish blood flow associated with cardiomyopathy led to clot formation. At the time of hospitalization, the clot size was 2.9 × 1.4 cm, but 26 days later, it diminished to below 0.5 cm during the present study. The flow analysis indicates the presence of a single large and strong LV vortex that penetrates effectively toward the apex where the clot is present. This induces the highest apical flow (⁠$v̂ind$ and $v̂apx$⁠) and very effective washing (Figs. 7–11, Supplemental video 4a-b which is available under the “Supplemental Data” tab for this paper on the ASME Digital collection). In fact, the apical velocity never goes below 4 cm/s at any point in the cardiac cycle (Fig. 11). The values of Y_v/Y_LV, A_v/A_LV, EPI, and VFT are comparable to those of healthy LVs [2,26,50]. Hence, the flow analysis predicts that the LVT will likely resolve and the patient will be at low risk of LVT relapse if anticoagulation is discontinued. The low Δ$Γ̂$ is also consistent with the EF. Remarkably, four months after this study, the patient called in and reported a full clot recovery and normal EF.

Hence, from these preliminary results, it appears that the distributions of $v̂ind$ and $v̂apx$ are promising metrics to evaluate the complete effect of the LV vortex on LVT formation. The former combines the effect of LV strength, presence of multiple vortices, as well as their location and duration. The histograms in Fig. 11 provide a quantitative assessment for the postulated differences in the extent of apical washing among the four patients. Such histograms could be used for obtaining a single parameter characterizing the efficacy of washing that combines, e.g., the velocity magnitude and its duration or the fraction of the cycle with velocity exceeding a yet-to-be determined threshold velocity. Determining such a value requires a substantially larger cohort of patients at varying levels of LVT status than the present four cases.

This section summarizes several limitations of this study. First, inherently, all 2D planar (optical and echo) PIV measurements determine the distribution of velocity components in the sample plane, i.e., they represent a 2D projection of a 3D flow. As long as the tracers remain in the measurement plane during consecutive exposures, the in-plane projection of the displacement can be measured. The present sample plane (A4C view) contains the dominant velocity component of the mitral inflow [36], reducing, but not eliminating, the likelihood that the out-of-plane motion in certain regions is strong enough that the tracers escape between exposures. Using PTV, in particular, increases the likelihood that the measurements are based on tracers that remain in the sample area in consecutive exposures. In regions with strong out-of-plane motion, one has to realign the plane or reduce the time delay between exposures to obtain the in-plane velocity. If needed, contrast echo data could be recorded in multiple planes/orientations, and then phase averaged and integrated to resolve 3D features of the flow structure.

Second, the image acquisition rate of clinical ultrasound machines when covering the entire LV is typically limited to 60–80 fps, restricting the maximum measureable velocity. This issue has already been discussed in prior studies [55,56]. For example, for an interrogation window size of 32 pixels, magnification of 0.3 mm/pixel, and displacement of half window, the velocity is 0.38 m/s. While one could extend the window, or shift the correlated windows, the spatial resolution would be compromised. The thickness of the ultrasound beam, presumably ∼ 1.7 cm [57], also limits the spatial resolution of the data. The maximum velocity may affect the data in the center of the LV, hence, the values of $Γ̂$ and $v̂ind$⁠. However, it is less likely to affect the measurement in low speed regions, emphasizing the importance of measuring the apical washing directly using $v̂apx$⁠. It should be noted that new ultrasound machines with acquisition rates exceeding 200 fps are becoming available, presently only as research tools.

Third, inherently, contrast traces must be present for performing velocity measurements. Hence, one may encounter challenges to obtain useful data in regions with poor circulation, where the contrast concentration is low. For all the present patients, the images used for analysis are acquired after the peak concentration phase, when the number of tracers drops to a level that does not flood the entire LV and can be traced individually. In cases with poor circulation one has two potential options: (i) use data obtained during period of much higher concentration to characterize the flow in regions of poor circulation, and then combine the resulting phase-averaged data with that obtained in better seeded areas, and (ii) rely on the LV vortex induced velocity to estimate the extent of local washing.

Fourth, the number of patients included in this study is small, excluding any comparisons with patients with healthy LVs or CM patients without LVT. Moreover, the clinical conditions of the four patients considered are also significantly different. Hence, the time history of flow structure and magnitudes between patients differ substantially, even after normalizing the data using the HR and LV dimensions. This limitation prevents us from reaching general conclusions about the significance of different flow features beyond those affecting apical washing.

Fifth, the uncertainty estimates for the instantaneous velocity measurements are not as small as those typically achieved in optical PIV [45]. For example, in the instantaneous map shown in Fig. 6, the uncertainty is about 10% of largest velocity, and smaller than only 94.5% of the velocity vectors. However, it decreases upon phase averaging. For the present data, which are based on averaging 10–15 cycles, the uncertainty in velocity magnitude is three times smaller, and the percentile of vectors larger than it increases to well above 99.5%. Hence, diagnosis based on phase-averaged data improves the reliability of the conclusions.

Finally, the present data do not have the resolution to measure the distribution of LV wall shear stresses. Hence, it cannot be integrated into coagulation models. However, it could be used for LVT risk assessment. Based on the significantly larger database discussed in Son et al. [7], and the limited cases considered in this study, it appears that the LV vortex parameters are correlated with the presence of an LVT. Clearly, more comprehensive studies with larger patient cohorts coupled with validated flow simulations and biochemical models are required to sufficiently elucidate these effects.

This study introduces an integrated series of processing tools aimed at calculating time-resolved 2D velocity distributions in the left ventricle over the entire cardiac cycle from contrast echo images obtained during routine clinical scans. The integration encompasses image enhancement, PIV, and PTV to optimize the quality and extent of data. These procedures have been implemented to study the LV flow in four patients with a history of past or current LVT. Distributions of phase-averaged velocity and vorticity reveal the evolution of the LV vortex properties, including peak vorticity, strength, area, and penetration depth. The inflow is associated with a sharp increase in $ω¯zmax$⁠, but the vortex strength and penetration continue to increase, peaking at the IC phase. Subsequently, the strength decays to a minimum at the end of systole. The present data satisfy the condition on average vortex penetration (Y_v/Y_LV) for the presence of LVT as reported previously [7]. They also confirm that Δ$Γ̂$ increases with EF [8]. However, high vortex penetration (Y_v/Y_LV) alone does not assure high apical flow due to decay and fragmentation of the LV vortex as it migrates away from the mitral annulus. Hence, the total vortex-induced apical velocity (⁠$v̂ind$⁠) is introduced and compared with direct measurements of the tangential velocity near the apex (⁠$v̂apx$⁠) as parameters to quantify the apical washing. It is advantageous to present both in the present introductory study since each provides different ways to quantify the apical washing. In fact, for two patients with high Y_v/Y_LV and resolved LVT after anticoagulation, the extent of washing in the vicinity of the past thrombus is low over the entire cycle. Hence, indicating high risk of thrombus relapse if anticoagulation is discontinued. This observation is consistent with their low EPI. Only one case involves a high vortex induced velocity in the thrombus area, in spite of low Δ$Γ̂$ and EF, but in agreement with a high EPI.

Although limited in scope and in need for verification in a larger clinical cohort, the present findings demonstrate the advantage of examining the time evolution of the entire LV flow structure, especially in regions prone to LVT. This kind of information can help guide physicians in making anticoagulation therapy decisions for cardiomyopathy patients who are at risk of LVT. Such optimized processing tools could be developed into clinical tools that can extract flow structure from routinely acquired contrast echo images. Their application may be extended to study a broad range of pathologies beyond the LV.

The authors would like to thank Professor Rajat Mittal from JHU for motivating this study; Susan Phillip, Kati Kashmer, Erin Jensen, and Nicole Soellner from the JHU Echo lab for help with acquisition and Professor Gianni Pedrizzetti (University of Trieste) for valuable inputs on echo-PIV processing. Data are publicly available through the Gulf of Mexico Research Initiative Information & Data Cooperative (GRIIDC) available at the website link² (doi:10.7266/N7HH6H4H).

• Gulf of Mexico Research Initiative.