Automated retinal microvascular velocimetry based on erythrocyte mediated angiography

: Retinal blood ﬂow is an emerging biomarker in ocular and systemic disease. Erythrocyte mediated angiography (EMA) is a novel technique that provides an easily interpretable blood ﬂow velocity quantiﬁcation by directly tracing individual moving erythrocyte ghosts over time in vivo, imaged using a scanning laser ophthalmoscope (Heidelberg Retina Angiograph platform). This tracking procedure, however, requires time-consuming manual analysis to determine blood ﬂow. To overcome this current bottleneck, we developed an objective and automated velocimetry approach, EMA - Automated Velocimetry (EMA-AV), which is based on a modiﬁed sequential Monte Carlo method. The intra-class correlation coeﬃcient (ICC) between trained human graders and EMA-AV is 0 . 98 for mean vessel velocity estimation and 0 . 92 for frame by frame erythrocyte velocity estimation. This study proves EMA-AV is a reliable tool for quantiﬁcation of retinal microvascular velocity and ﬂow and establishes EMA-AV as a reliable and interpretable tool for quantifying retinal microvascular velocity.

. The procedure of Erythrocyte mediated angiography. from each subject. In a sterile environment, erythrocytes were separated via centrifugation and then transferred to a hypotonic solution, leading to the transient opening of pores in the cell membrane, through which ICG entered the cells. Subsequently, the osmolality of the solution was increased to initiate closing of cell membrane pores, leaving the entrapped ICG dye within the erythrocyte ghosts. After microscopic examination of the erythrocytes, and following pupillary dilation of the subject, 1mL blood with autologous ICG-loaded cells were injected intravenously.
As this technique involves the re-injection of autologous erythrocyte ghosts, there is a potential for non-physiologic disturbance of the vasculature. We note that Flower et al completed a detailed assessment of erythrocyte fragility and integrity and found that erythrocyte ghosts made using this procedure are comparable to control erythrocytes [41]. As described by Flower et al [41], and from the experience of our group, this is a well-tolerated procedure without any ocular morbidity.

Image acquisition
For fluorescence angiography imaging, two confocal laser scanning ophthalmoscopes by Heidelberg Engineering GmbH were used, a modified Heidelberg Retina Angiograph 2 (HRA2) and a Spectralis HRA. These devices generate fluorescence angiogram movies by repeatedly scanning over the retina with a laser diode that excites ICG fluorescence. For peak ICG excitation, a diode laser with a wavelength of 786nm was used, and a barrier filter at around 800nm edge wavelength separated excitation and fluorescent light. The modified image acquisition systems limits the vertical field of view (FOV) in order to allow for an increased imaging frame rate at 24.6 frames per second (fps) during acquisition of a 15°× 7.5°(horizontal and vertical) FOV. For achieving an even higher frame rate of 46.1 fps at a 30°× 5°FOV, a sinusoidal scan pattern was employed vertically. The slower scan speed yielded a digital resolution of approximately 11µm, while the faster acquisition yielded twice this digital resolution at approximately 5.5µm, depending on the specific refractive properties of the eye being imaged. The diffraction limited optical resolution of the system is approximately 5.4µm (Airy disk radius). After acquisition, an image sequence was exported from the device software, geometrically corrected in case of sinusoidal scanning, and subsequently analyzed as described below.

Clinical data collection
Subjects for this study were recruited from the ophthalmology and/or optometry service at the University of Maryland. All experiments were conducted in accordance with the Declaration of Helsinki, and were approved by the institutional review board at the University of Maryland. Detailed information is provided in Section 3.2.

Image preprocessing
Prior to application of the proposed method for speed estimation, all video frames must be preprocessed, as illustrated in Fig. 2. To accomplish this, firstly, images are de-noised by applying a 5 × 5 Gaussian filter. The Gaussian filter is a good approximation for the point spread function of the fluorescence spots, which can also improve object visibility [42]. The image registrations are conducted based on the intensity-based automatic image registration algorithms in MATLAB 2018b, using the 'imregister' function and 'multimodal' configurations. After de-noising and registration, the images are denoted as The vessels of interest are delineated and labelled by trained graders in a manual fashion prior to automated velocimetry analysis. In order to remove background noise and simultaneously enhance erythrocyte signals, the temporal average of the whole sequence, I avg , is subtracted from each individual raw image in the sequence Here T indicates the total number of frames in the video and I avg = T t =1 I t 0 T . A representative temporal average image for 24.6fps videos and 46.1fps videos is displayed in Fig. 3, where the vessel shape is conspicuous and can be used for labelling the vascular network.

Overview of the EMA-AV framework
EMA-AV is a Monte Carlo (MC) based method, and a diagrammatic representation of the method is shown in Fig. 4. Beginning with frame t k−1 , the EMA-AV method is constituted by two general steps: MC particle initialization and MC particle movement. With certain movement stopping criteria, the new frame, t k , is used for initialization. The entire tracking process is considered complete when all video frames have passed through the step of initialization. To explain the EMA-AV method, Section 2.4.3 -2.4.5 explicate the mathematical framework behind EMA-AV, and Section 2.4.6 describes the EMA-AV procedures in detail.

Stochastic dynamic system
The tracking problem can be conceptually understood as the task of recovering the real dynamics of object (erythrocyte) states (positions) from a sequence of noisy observations, which can be modelled by a set of stochastic dynamic equations, as shown in Eq. (1):  x where x k represents the state of the object at time t k . f k−1 is a transition function describing the evolution from x k−1 to x k , as a first-order Markov process. α k is an independent and identically distributed (iid) random process representing the process noise. h k defines the relationship between the state, x k , and the measurement, z k . θ k represents the iid measurement noise. In essence, the tracking problem aims to estimate the posterior probability density function (PDF), p(x k |z 1:k ).
The solution to this problem is, in a general sense, constituted by two steps: prediction, and update [28]. In the prediction step, x k can be estimated from p(x k−1 |z 1:k−1 ) and transition density π k |k−1 (x k |x k−1 ), based on the Chapman-Kolmogorov equation shown in Eq. (2): When the new measurement is obtained, the posterior PDF can be updated according to Bayes rule, based on Eq. (3), where g k (z k |x k ) is the likelihood function describes how likely the noisy observation signal can be the signal from targeted object:

Sequential Monte Carlo (SMC)
The SMC method establishes a general framework for a numerical solution to the dynamic system problem [40]. It utilizes a large number of random weighted MC particles as an approach to optimal estimation of the posterior PDF. We extend the idea to estimate the object movement velocity (Section 2.4.5). The description of the general SMC method is given as following: Suppose at time t k−1 , the posterior density is approximated by N weighted MC particles in: , where x i k−1 is the state of particle i in time t k−1 , and w i k−1 is its corresponding weight (which must be normalized to describe the possibility of being a real object). The detailed definition of w i k in our application is given in Section 2.4.6. The posterior PDF at time t k−1 can be expressed as Eq. (4): where δ(x) is the Dirac function. At time t k , the posterior PDF can be approximated by a new set of weighted MC particles , as shown in Eq. (5): where In general cases, q k ≡ π k |k−1 . To avoid particle weight degradation [28], a potential problem where MC particles will be less likely to represent the real signal when tracking, a re-sampling process is ordinarily required to preserve more large weight particles and discard more small weight particles. The re-sampling step selects N particle(s) from , with the probability of selecting particle i being w i k [43].

SMC for erythrocyte velocity estimation: EMA based automated velocimetry (EMA-AV)
In SMC, the transition function f is based on prior knowledge of object velocities [36]. In the tracking procedure, the function can gradually filter out uncertain signals with low signal weights. However, with regards to our erythrocyte speed estimation application, the transition function (object velocity) required estimation frame by frame. To that end, we borrowed the idea of weighted MC particles to consider signal uncertainties, sacrificing the filtering ability of the SMC method. The weighted average velocity of MC particles is a better speed estimation. To ensure that the MC particles represent multiple erythrocytes, we introduced spatial clustering of MC particles in Section 2.4.6. Each MC particle cluster represents an erythrocyte, whose velocity, v k−1 , at time t k−1 , needs to be estimated based on p(v|z k−1 , z k ). Once v k−1 is determined, the transition density, π k |k−1 , may be deduced without considering the process noise, α. Utilizing the weighted MC particles in an erythrocyte, p(v|z k−1 , z k ) can be computed as follows: In order to compute the integration in Eq. (7), we allocate v x i k−1 to each MC particle, x i k−1 , and transform the continuous integration into a discrete integration. v x i k−1 is sampled from the uniform distribution in Eq. (7). A one-to-one mapping between the set {w i can be established naturally. Equations (6) and (7) can also be simplified as: 2.4.6. Elaboration on the steps of EMA-AV As shown in Fig. 4, the EMA-AV method has two main steps: MC particle initialization and MC particle movement. MC particle initialization ascertains the weighted particle set {w i , as mentioned in Section 2.4.4. The N weighted MC particles represent N pixel signals from an unknown number of erythrocytes. Their weights represent how likely the signal is to be originating from erythrocytes. Firstly, the particles were selected in a random and uniform fashion within the vessel of interest, and w i k−1 was set as 1 N . In accordance with the gray-scale value in I k−1 (observation), the image from the pre-processing step, the MC particle weight values could be updated. The updating of the particle weight values was accomplished according to Eq. (10) based on the Gaussian intensity distribution assumption of likelihood function g k (z k |x k ), as mentioned in Section 2.4.3 [40].w i k−1 were subsequently normalized according to Eq. (5): Here, I max and σ were both hyper-parameters. If the grayscale value of a particle was larger than I max , the non-normalized weight of the MC particle was maximized. In the experiment, the I max was set as 100, and sigma was set as 30. To exclude invalid frames associated with a blink, if the largest particle grayscale value was lower than an established hard-set threshold, this frame was discarded.
The re-sampling process was conducted in following the precaution stated earlier in Section 2.4.4, namely, to avoid MC particle weight degradation. After re-sampling, low weight MC particles were usually discarded. However, there was still a potential for several MC particle clusters to be located within the vessel. These MC particle clusters, which could represent various independent erythrocytes, were amenable to have their movement velocities computed independently.
In order to define the MC particle clusters mathematically, a mean shift clustering algorithm [44] was applied to separate the particles as dictated by their two-dimensional (2D) coordinates. As current, cluster analysis is an important branch in the domain of pattern recognition [45]. However, as opposed to classic centroid-based clustering methods, such as the "k-mean algorithm", mean shift algorithms do not require an explicit preliminary designation of the number of clusters. This was appropriate considering the uncertainty of the exact number of erythrocytes in each frame. All the MC particles in each mean shift cluster represent an erythrocyte. In this study, the bandwidth of the mean shift algorithm was set to 5.
Spatial clustering could generate C MC particle clusters. Given cluster c, the set of weighted particles could be denoted as {w i where N c was the number of particles in cluster c. As previously introduced in Section 2.4.5, each particle could be assigned a random velocity value v x i k−1 , and the particles could move accordingly within the 2D region of the vessel of interest. If all the particles exited the target vessel, this corresponding particle cluster was removed from the cluster list. After movement, the particle positions were denoted as {x i k } N c i=1 , and their weights, {w i k } N c i=1 , were re-calculated according to the gray-scale value in the next frame, I k . Re-sampling and spatial clustering operations were conducted in a similar fashion to the former step. If spatial clustering generated a new particle cluster (cluster split), the movement of cluster c ended. This was designed to address the concern of erythrocyte aliasing, which could have otherwise affected the experiment results (during manual labelling, trained human graders also excluded cases where there was any uncertainty as to the identification of cells from frame to frame). If no cluster split was observed, the cluster's movement velocity, v c k−1 , would be recorded following Eq. (9), and the movement continued.
When the movements of all C particle clusters had ended, the initialization and movement steps in Fig.4 needed to be restarted from frame t k . As introduced in 2.4.2, the entire tracking process was considered complete when all video frames had passed through the step of initialization. The EMA-AV framework is implemented in Python 3.6. It requires approximately ten minutes to process an EMA video (with 300 frames) using an Intel i7-7700 CPU (by comparison, more than eight hours are required for manual labelling by trained human graders). The python class of EMA-AV has been uploaded to http://taolab.umd.edu/research/

Tracking results visualization
In order to present the erythrocytes tracking results in an intuitive manner, the typical tracking chains in the 46.1fps and 24.6fps EMA videos are displayed in Figs. 5 and 6. In each frame, one venule and one arteriole are shown. The practical purpose of these images are to allow readers to appreciate the reliability of automatic erythrocyte speed estimation. Examples of angiograms with tracked erythrocytes using EMA-AV are included in the supplemental materials (Visualization 1 and Visualization 2). The topmost images in Figs. 5 and 6 are variance maps of corresponding EMA videos, with the target vessel skeleton labelled. The variance maps display increased contrast of smaller vessels compared to temporal average frames, as defined in Section 2.4.1, but have more noise for image preprocessing.

Statistical analysis
In total, 11 vessels (5 venules and 6 arterioles) from 4 subjects (5 individual eyes) were used for the erythrocyte tracking algorithm validation. All vessels were imaged by the modified HRA2 and Spectralis, in both 24.6fps and 46.1fps modes, within the same day. Trained human graders were appointed the task of labelling and tracking erythrocyte movements on a frame by frame basis. Manual labelling was found to be a time-consuming process which thereby limits large scale applications and experiments. However, based on previous studies, grading results from different graders showed the intra-class correlation coefficient (ICC) [46] value to be 0.99 [23]. Therefore, in this study, results obtained by trained human grader tracking served as a basis for comparison with the data yielded by EMA-AV.

Mean erythrocyte velocity analysis
Mean erythrocyte velocity is an important index to study overall vessel blood flow [22], which is a tool that may promote the understanding of blood flow alteration in eye disease. In EMA, the mean velocity was computed based on all tracked erythrocytes. In comparing the results obtained by EMA-AV with those from the trained graders, several modalities of data analysis were employed (ICC, coefficient of determination (R 2 ) and the mean coefficient of variation (CV)), and these are summarized in Table 1. Statistical analysis was performed in RStudio (version 3.5.2). The linear regression plot is displayed in Fig. 7. In Fig. 8, Bland-Altman plots are presented to visualize the difference between the results obtained by EMA-AV versus the trained graders. Based on these data, in addition to the results of a paired t-test, it is evident that no statistically significant difference exists in mean velocity results as determined by trained graders and EMA-AV. This conclusion supports the proposition that EMA-AV may serve as a reliable and accurate alternative to trained human grader reporting of mean erythrocyte velocities as based on EMA.
EMA-AV displayed a high degree of correlation with the trained graders results in both 24.6fps and 46.1fps videos. 24.6fps results showed higher correlation than 46.1fps results, based on experimental data. A potential explanation could be that the mean signal to noise ratio in the 24.6fps videos was 11.32, which was 33% higher than the SNR in 46.1fps videos (SNR=8.51). In Fig. 8(b), it is shown that the outlier data points, with regards to the limits of agreement, were from the video with an SNR of 4.33. In such conditions, even single particle tracking in simulated videos cannot be accomplished [35].

Frame-based erythrocyte velocity analysis
By means of EMA, it is also possible to observe blood particle component dynamics under normal physiologic conditions. Moreover, EMA facilitates the study of the relationship between retinal ocular blood dynamics and ocular disease [22]. In order to evaluate individual cell tracking performance, frame based ICC values between trained graders and EMA-AV were computed. In deciding whether to include or exclude data reported by the graders and EMA-AV for ICC computation, it was first assessed whether the two tracking media had produced mutually agreeable results. If both the trained graders and EMA-AV mutually reported an erythrocyte tracked from frame t to frame t + 1, its velocity was then used for ICC computation. If either a trained grader or EMA-AV reported an erythrocyte tracked (but there was no mutual agreement), then the velocity would be excluded from ICC computation. A potential reason for lack of agreement was due to a tendency for human graders to overlook erythrocytes in low contrast situations. EMA-AV would reject particle movement when erythrocyte cluster splitting had occurred, as displayed in Fig. 4. To clarify, all particles were included when computing the mean erythrocyte velocity in Section 3.2.1.
Frame-based ICC values of 11 vessels, in both the 24.6fps and 46.1fps videos, are displayed in Tables 2 and 3. Based on the information displayed in the tables, it could be concluded that ICC values vary between vessels, potentially due to higher velocities related to vessel size or also possibly resolution. Further, it was possible for the variable ocular physiologic status to influence EMA image quality and thus alter frame-based ICC values. Compared to 46.1fps videos, 24.6fps videos had higher SNR values than 46.1fps videos, as previously stated in Section 3.2.1, while the spatial resolution of 24.6fps videos was inferior that of 46.1fps videos, as mentioned in Section if there was a potential erythrocyte cluster in frame t and two in frame t + 1, and the two clusters in frame t + 1 had a significant grayscale value difference. In this situation, the trained grader displayed a tendency to link the cluster in frame t with the closer cluster in frame t + 1, while, conversely, EMA-AV tended to link it with the brighter cluster in frame t + 1. These tendencies may be understood as systematic differences in EMA methods, which, as discussed in Section 3.2.1, did not necessarily influence mean velocity estimation performance statistically. EMA-AV can generate a plot of the cardiac systolic and diastolic cycles based on frame by frame erythrocyte velocities, and an example is shown in Fig. 9. The limitations of this method to generate such a plot are that the velocity of each frame is determined by individual erythrocytes and not gated with an ECG signal at the current stage. Further investigations in frame by frame erythrocyte velocities with ECG gating may permit understanding of systolic and diastolic cycles and their importance in the relationship between retinal blood flow and cardiovascular (or other systemic) disease(s). Fig. 9. An example cardiac systolic and diastolic cycle plot (venous). In order to de-noise the signal, the results were smoothed every 0.5s.

Comparison between 24.6fps EMA and 46.1fps EMA
As mentioned in Section 2.1, the primary distinction between 24.6fps and 46.1fps Heidelberg SLO imaging was the change in scanning pattern. For both trained grader and EMA-AV results, when assessing all 11 vessels, there was no statistically significant difference with regards to the mean erythrocyte velocity between the two frame rate settings (human: p = 0.179, EMA-AV: p = 0.277). However, for each vessel, when comparing the average top 10% fastest velocities, the values in 46.1fps EMA videos were significantly higher than the values in 24.6fps EMA videos for both human graders (p = 0.986) and EMA-AV (p = 0.963). The top 10% fastest velocities could be of greater importance than mean erythrocyte velocity when assessing faster velocities such as those occurring during cardiac systole. Moreover, based on the aforementioned data, it could be observed that 46.1fps EMA was capable of capturing more accurate faster velocity data. The EMA-AV algorithm was capable of effectively ascertaining the average top 10% velocities in both 24.6fps video and 46.1fps videos, and, when compared to the trained graders, the ICC values were 0.98 in both frame rates. The corresponding Bland-Altman plots are shown in Fig.  10. As shown in the figure, the average top 10% velocities derived from the trained graders were generally higher than those from EMA-AV. An explanation for this phenomenon is that in Eq. (7), v min and v max were required to be pre-defined to restrict the search range of the SMC method.
To further refine the intra-setting variability of the two different imaging frame rates, one arteriole was excluded whose diameter was > 70µm, as this larger arteriole would be too fast for velocity detection at 24.6fps, and its CV value between 24.6fps and 46.1fps videos was larger than 0.15. The vessel diameter range for the remaining 5 arterioles was 30µm ∼ 55µm, and 22.5µm ∼ 50µm for the 5 venules. The vessel diameters were computed based on the temporal average image from each video sequences. Our study found the mean velocities in small vessels from EMA-AV were more robust to different frame rates as compared to trained graders. For EMA-AV, the intra-setting R 2 value for mean erythrocyte velocity was 0.856, and the ICC value was 0.9553. For trained graders, the intra-setting R 2 value was 0.8057, and the ICC value was 0.9435. 3.3. EMA-AV: pros, cons, and comparison to other methods EMA-AV offers the benefit of direct measurement of erythrocyte ghosts in vivo. This comes with several potential benefits and drawbacks. Of these, the most promising benefit is the potential to accurately and precisely characterize absolute erythrocyte velocity in the retinal microvasculature. The invasive nature of this technique is a clear drawback of this procedure, as it involves autologous injection of ICG-labelled erythrocyte ghosts. As compared to other forms of invasive angiography, however, EMA-AV has a favorable risk profile. ICG is relatively safe and less likely to result in allergy or nausea and vomiting as compared to fluorescein angiography [47]. Furthermore, the amount of ICG in ICG labelled erythrocytes is approximately 1/700 of that in conventional ICG, which may make it safer. Adaptive Optics Scanning Laser Ophthalmoscopy (AO-SLO) and other Adaptive Optics techniques offer a precise and noninvasive method of determining blood flow in a small field of view generally using equipment that is not commercially available. EMA-AV is clearly more invasive, but allows for determination of blood velocities in a wide field of view using a modified commercially available device. OCT-Angiography currently allows for mapping out of blood vessels, and determination of relative blood flows using techniques such as Variable Interscan time analysis (VISTA) [48]. More established such as Color Doppler Imaging (CDI), laser speckle imaging, and laser doppler flowmetry allow for noninvasive determination of relative blood flows as compared to absolute blood flows and have generally shown higher variability. For example, the Canon Laser Blood Flowmeter has been shown to have measured coefficients of variation for blood flow of up to 39.7% [49].

Conclusion
Erythrocyte mediated angiography (EMA) is a promising method for quantifying blood flow via direct observation of erythrocyte movement. In an effort to overcome the time-intensive process of trained human grader erythrocyte tracking, this study proposes an automatic and objective method of analyzing EMA videos, called EMA-AV (EMA-Automated Velocimetry), which is capable of ascertaining erythrocyte velocities. EMA-AV applies the concept of the sequential Monte Carlo method to determine erythrocyte velocities. In two frame rate settings (24.6fps and 46.1fps), EMA-AV derived mean erythrocyte velocities and frame by frame velocities were highly correlated with trained grader tracking results, as supported by the aforementioned statistical analyses, which included an ICC, R 2 and a CV.
Erythrocyte velocities obtained by EMA-AV have potential applications in studying the relationship between retinal blood flow dynamics and the cardiac cycle, including systolic and diastolic changes. 46.1fps EMA videos were more accurate for describing faster systolic erythrocyte velocities than 24.6fps EMA videos. Moreover, EMA-AV demonstrated the capability of capturing faster velocities as accurately as trained graders in both frame rate settings. Despite the comparisons between trained human graders and EMA-AV being highly dependent on the EMA method and image sampling characteristics, a reliable automatic analysis method could still benefit the further validation of EMA.
Future development of EMA may involve both hardware and software improvements. With regards to hardware, an increase in scanning rate with a minimal consequential reduction in SNR is a potential means of broadening the erythrocyte velocity estimation range, which could further enhance the robustness of EMA. In the current state of EMA image analysis, trained human graders manually label the target vessel in the pre-processing step, and, as a result, EMA-AV may be regarded as a semi-automatic modality. To further increase automation, image segmentation techniques may be used. In the interim, there is scope for enhancement of image particle description by augmentation with artificial intelligence techniques.
In conclusion, EMA-AV is a reliable prototype method for EMA video analysis and erythrocyte velocity estimation. It is expected to broaden the applications of EMA and facilitate quantitative analyses of retinal microvascular dynamics.