Transthoracic ultrasound localization microscopy of myocardial vasculature in patients

Myocardial microvasculature and haemodynamics are indicative of potential microvascular diseases for patients with symptoms of coronary heart disease in the absence of obstructive coronary arteries. However, imaging microvascular structure and flow within the myocardium is challenging owing to the small size of the vessels and the constant movement of the patient’s heart. Here we show the feasibility of transthoracic ultrasound localization microscopy for imaging myocardial microvasculature and haemodynamics in explanted pig hearts and in patients in vivo. Through a customized data-acquisition and processing pipeline with a cardiac phased-array probe, we leveraged motion correction and tracking to reconstruct the dynamics of microcirculation. For four patients, two of whom had impaired myocardial function, we obtained super-resolution images of myocardial vascular structure and flow using data acquired within a breath hold. Myocardial ultrasound localization microscopy may facilitate the understanding of myocardial microcirculation and the management of patients with cardiac microvascular diseases.

The relation between Shannon theory and the tracking can be revealed by simple cases, where multiple MBs move in a straight vessel, as shown in Fig. M1a.To correctly pair MB 2 by the nearest-neighbour method, the distance between the positions of MB 2 in two frames should be smaller than the distance between MB 2 in the second frame and MB 3 in the first frame.
where dm is the moving distance, the ratio of flow speed to frame rate, and db is the MB distance which is inverse to the MB concentration.According to the above equation, a higher frame rate is needed to track MBs at higher concentrations and faster speeds.Motion model can predict the MB position along the flow direction, as shown in Fig. M1b.To correctly pair MB 2 in this case, the distance between the MB 2's predicted position and actual positions should be smaller than the distance between the MB 3's predicted position and MB 2's actual position, where dp is the MB moving distance predicted by a motion model.With the motion model, tracking method can correctly pair faster MBs at the same frame rate and MB concentrations, compared to the nearestneighbour method.However, it doesn't mean the larger the dp is, the higher accuracy the tracking can achieve.To avoid pairing MB 2 in Frame 2 to MB 1 in Frame 1 in the case shown in Fig. M1c, dp should satisfy Combining Eq. ( 2) and (3), we can obtain With the frame rate f, the equation can be rewritten as where vm is the flow velocity and vp is the predicted velocity and vm/db can be regarded as the spatial repetition frequency, defined by how many times an MB passes through the positions of the other MBs in one second.Then, the temporal frequency, frame rate or pulse repetition frequency f, should be at least double the spatial repetition frequency to correctly pair the MB by nearest-neighbour.Using the motion model, the spatial repetition frequency of the MB movement vm/db can be reduced by vp/db.
Fast-moving MBs at high concentrations can be more accurately tracked using motion model than only using the nearest-neighbour, especially when the frame rates are restricted by physical imaging depth and the number of angles used for compounding.In our previous work [1], we proposed a feature-motion-model framework to achieve accurate and efficient MB tracking, which combined MB image feature and Kalman motion model in a cost function and paired MBs by finding the total minimum cost in a graph-based assignment.The Kalman motion model's state vector x, translation matrix F and observation matrix H is defined as

Fuzzy initialization
where x and z are the lateral and depth location of MBs in the image; vx and vz are the MB moving speeds along the corresponding directions.vx and vz, i.e. vp, were set to 0 for the newly appeared MBs as we could not get the MB velocity before pairing.As demonstrated in the last section, wrong pairing is still likely to happen for new MB with this initialisation where motion model does not take effect by assuming new MB to be static.Therefore, we proposed fuzzy initialisation and combined it with our previous framework in the way shown in Fig. M2a.The fuzzy initialisation was done for the newly appeared MBs using three consecutive frames, as shown in Fig. M2b.The fuzzy initialization was done with the below steps.
Cost matrix construction.To initialise the M new MBs in the first frame, a 3D cost matrix for paring the M new MBs with the MBs in the following two frames was built.Various information can be used for calculating the cost.In this study, the cost consisted of three components: innovation of motion model, intensity differences, and intensity variances along the candidate tracks.
We estimated vx and vz for each new MB in the first two frames, and therefore the innovation of motion model could be calculated by the distance between the observed MB locations and the predicted MB locations in third frame where the superscripts, m ≤ M, n ≤ N and q ≤ Q, denote the MB indexes in the three frames.When giving a search window, model innovations for pairs out of the searching window were given as infinity.
The intensity difference dI among MBs in three frames was calculated as where I is the MB image intensity.
The maximum intensity projection (MIP) of contrast enhanced ultrasound (CEUS) sequence presents the vasculature in a low resolution.While pairing two MBs, the trajectory can be along or across two vessels.It can be assumed that pixel intensities on the trajectory along one vessel have smaller variance than those on trajectories across vessels.Therefore, the intensity variances of MIP pixels on each trajectory can be used for calculating the cost of pairing and defined as where the numerator is the standard deviation of pixel intensities on a trajectory and the denominator is the average.To save computation, the trajectory between one pair was defined as a straight line; seven points were equidistantly sampled on the trajectory between MB pairs in the first two frames and then assigned to candidate pairs among three frames which shared the same pairing between the first two frames; the pixel intensity of the points on each trajectory was obtained from the nearest MIP pixel.Three components were in the same size of 3D matrix and combined to a 3D matrix by the first-order principal components of the principal components analysis (PCA).
Finding the minimum.The minimal cost was found for each new MB independently, instead of finding the total minimum of all the pairs to save the computation.The innovation distance d m and velocity vector v m between MBs in the first two frames corresponding to the minimal cost were saved for each new MB.

Repetition within different sizes of searching windows.
To increase the robustness of the initial paring, the above two steps were repeated for different search window sizes.In this study, the searching window sizes were varied between 0.3 to 1.2 times the maximum MB moving distance between adjacent frames (assuming maximal blood flow velocity of 150 mm/s).As a result, for the r th searching window, d m,r and v m,r were saved for the mth new MB.
Initialising Kalman state vector.An intuitive way to combine multiple v m,r into one for the new MB is averaging.However, the distribution of v m,r may not be Gaussian, and some extreme cases might significantly affect the average.As the distance d m,r can be regarded as errors of the linear motion model, the normalised histogram, shown in Fig. M2b, represents the probability of the certain value of errors when assuming linear MB movements for this data.Therefore, the velocity of the new MB can be calculated by the below weighted average to reduce the effect of extreme cases where w m,r is the number obtained by indexing the histogram by d m,r corresponding to the v m , and R is the number of searching window sizes used.A new bubble might not be paired with any other bubble in a small searching window.In this case, the searching window was excluded from the average in Eq. (10).Only the new bubbles were paired in more than two searching windows were initialised by Eq. ( 10), and otherwise were set as static.
In the above initialisation, the speed components in the state vector were obtained by the weighted average across the different searching window sizes.As a result, a velocity vector from an initialised MB might not point to any MB in the next frame.Therefore, we named the above procedures as fuzzy initialisation.
After being initialised, new MBs in current frame and MBs existing in previous as well as current frames were paired with all the MBs in the next frame by finding the total minimum of cost through the graph-based assignment.The cost for the graph-based assignment was defined as the ratio of MB image intensity difference to the probability obtained through motion model by p(ok|xk|k−1) where the subscript n|n' denotes the estimate of the value in frame n given observations up to and including frame n'; ok is measured MB state vector, consisting of lateral and depth coordinates; Sk is the innovation covariance.

Kalman motion model parameters
Covariance of measurement noise R, covariance of model prediction noise Q and covariance of estimation noise P are parameters that need to be set for the linear Kalman motion model and sometimes can significantly reduce tracking accuracy with inappropriate values.R can be set by the localization uncertainty of the imaging system, which can be measured by in vitro experiment and estimated for in vivo acquisition.Although Q can be manually chosen according to user experiences, we proposed a method to estimate Q from data to avoid human interaction and reduce effort in tunning parameters when dealing with various vasculature and haemodynamics.
As presented in the last section, a linear motion model is used to predict MB movement in the third frame using the first two frames.Steps from 1 to 3 in the last section can be used for to estimate Q for each three adjacent frames.Instead of saving distance between the prediction and observation for each MB, the discrepancies along lateral e x,m = (x q +x m −2x n ) and depth e z,m = (z q +z m −2z n ) were saved for each MB in the first of the three frames.By assuming most of pairing after these three steps are correct, the discrepancies can be used to approximate the measurement pre-fit residual of Kalman filtering, the distance between predicted and measured MB locations.Therefore, we used the variance of the discrepancies to calculate the pre-fit residual covariance.
where  9: $ and  9; $ are the variance across e x,m and e z,m respectively.As variance estimation becomes more accurate with a larger number of data, discrepancies were obtained from ten groups of three consecutive frames picked from the image sequence, and a variance  9 $ was calculated with all the discrepancies along the lateral and axial directions and used for  9: $ and  9; $ .
The S in one frame was defined by where Pk−1|k−1 and Pk|k−1 are the updated and predicted estimation covariance matrices respectively.Using optimal Kalman gain, Pk−1|k−1 is independent from the measurements and can achieve an asymptotic value as  < = ( < −  <  8 ( <  8 + ) '  < ) 8 +  (14) Therefore,  < is determined by the defined R and Q and used for Pk-1|k−1 to make MB Kalman filter independent from how many frames the MB have appeared in the imaging plane.The time continuous noise model [2] was used and thus where q is the variance of error when assuming constant velocity in the linear Kalman motion model.As analytical solution of Eq. ( 14) was not found, Q was estimated using Eq. ( 15) while q was calculated by the bisection method in Algorithm 1.
Algorithm 1 Find q Require: R, M O M, [ " ,  > ] searching range of q while  > −  " >  do q ←( " +  > )/2 Calculate  < with Q and R Calculate S with Eq. ( 13) if || > | O | then Ø || increases with q in the searching range  > ← q else  " ← q end if end while Return q

Fig. M1 |
Fig. M1 | Three cases to pair multiple MBs moving in one straight vessel.a, a case where MBs are to be paired without prediction of motion.b, a case where MBs are to be paired by motion model with predicted movement less than the moving distance.c, a case where MBs are to be paired by motion model with predicted movement more than the moving distance.Number in each disk denotes the identity of an MB.db is the distance between MBs in one frame.dm is the moving distance of the MB between two frames.dp is the moving distance predicted by motion model for MBs.

Fig. M2 |
Fig. M2 | a, a framework combining fuzzy initialisation and graph-based assignment for MB tracking.Fuzzy initialisation is done for new MBs detected in frame n by using MBs in frame n, n+1 and n+2.After initialisation, all MBs in frame n are paired with MBs in frame n+1 by our previously proposed feature-motionmodel tracking framework.b, fuzzy initialisation diagram.Moving velocities of a newly appeared MB is estimated with MBs in next two consecutive frames by finding minimum of pairing cost under different searching window sizes.Weighted average of the velocities obtained under different searching window sizes is used as the velocity components in the Kalman state vector of the new MB.v m,r is measured velocity between the first two frames corresponding to the minimal cost found for the mth new MB under the rth searching window size.Innovation distance d m,r in the third frame corresponds to the minimal cost found for the mth new MB under the rth searching window size.