1. Introduction

Retinal blood flow (RBF) is an potential biomarker for diabetic retinopathy [1], age-related macular degeneration [2], and glaucoma [3]. Moreover, systemic diseases, including hypertension [4], diabetes [5], hyperglycemia [6] and neuro-degenerative diseases such as Alzheimer’s dementia [7,8] are also associated with abnormal changes in dynamic RBF.

Methods to characterize ocular blood flow include conventional angiography with fluorescein/indocyanine green (ICG) [9, 10], color doppler imaging [11], laser doppler imaging [12] and laser speckle imaging [13]. Different OCT-based blood flow velocimetry approaches have also been introduced in the past for quantifying erythrocyte flow dynamics, but are currently still lacking a reliable and comprehensive validation in vivo [14, 15]. Another concept, known as “retinal functional imager (RFI)”, possesses the capability of determining erythrocyte velocity by assessing relative concentration of hemoglobin chromophores in both vessels and the capillary background [16]. However, reproducibility of the measurement hence derived has substantial variability [17], and more so for older subjects and glaucoma subjects [18]. Adaptive Optics - Scanning Laser Ophthalmoscopy (AO-SLO) allows for quantification of flow in the microvasculature within a small field of view [19]. In particular, near confocal AO methods [20] and flood AO methods [21] provide precise blood flow values in smaller fields of view. However, none of these methods have yet been established as an accurate and precise gold standard for measuring retinal blood flow velocity in the living human eye.

Erythrocyte mediated angiography (EMA) is an emerging imaging technique which is capable of directly measuring erythrocyte velocity in vivo [22]. In EMA, ICG-labelled erythrocyte ghosts are injected into the bloodstream, enabling blood flow velocimetry by following the motion of individual erythrocytes over time. This makes EMA-based velocimetry an easily understandable and interpretable technique that can be directly visually verified, even by non-experts in the imaging field. One major drawback of EMA is the relatively invasive nature of the technique. Still, however, EMA has the capability to serve as a generally accepted validation tool for in vivo measurement of blood velocity in the retina.

Current methods of velocity determination from EMA data rely on delineation by human experts [23], including erythrocyte identification and tracking, which is a timing-consuming process. This limits the potential adoption and utilization of the technique. As a result, the majority of current EMA applications consist of characterization of vasomotion, the spontaneous oscillation of blood vessels that may cause stasis and erythrocyte pausing [24]. Direct erythrocyte movement measurement has not been hitherto explored in detail. Some studies that characterize cell flow, such as applied particle image velocimetry [25], or spatial temporal plots [26], usually require high frame rates or slower object movement to ensure videos with continuous cell flow for accurate measurement. In comparison, EMA has a large field of view, relatively slower frame rate, and rapid object movement as compared to other technologies [23]. The large cell movements between any two frames underlies our study of an advanced and automated EMA analysis method, with the primary objective being to solve the object tracking problem. Automatic object tracking has an extensive history in the field of signal processing and computer vision [27]. This tracking problem could be considered as determining real object states given specific observation states by solving systemic equations [28]. In this regard, the mathematical tools that are available include the Kalman filter [29], the extended Kalman filter [30], and the particle filter [28]. The discovery of better object representations are another branch of tracking studies [31], especially with the recent advances of deep learning [32, 33]. Particle tracking is particularly important in bioimaging when studying living cell activities [34]. In comparison with images produced from general tracking tasks, biomedical particles lack rich image features and biomedical images are typically burdened by a lower signal-to-noise ratio (SNR), which accordingly brings about greater uncertainty. Researchers commonly separate the steps of particles detection and particle tracking [34], but this method is only beneficial for reliable particle detection and single particle tracking [35]. I. Smal et.al [36] unitized both spatial and temporal information to link detection and tracking under the constant velocity hypothesis. This assumption, however, is invalid when applied to erythrocyte tracking applications. Erythrocyte velocities are varied by changing vessel diameters [37], while the vessels themselves exhibit changes in diameter with respect to their relation to the optic disc [38]. Further, erythrocyte velocities are influenced by the systolic and diastolic cycles [39]. One additional difficulty in EMA, as opposed to other biomedical partical tracking tasks [40], is the potential for erythrocytes to enter/exit the region of interest, as well as frame-by-frame variation of erythrocyte intensity, thereby conferring an increased complexity to tracking. Variation in intensity may be in part due to axial movement of erythrocytes.

In this study, we propose a Monte Carlo-based method that links particle detection and tracking to automatically estimate the velocity of individual erythrocytes moving in small retinal arterioles and venules. To the best of our knowledge, this is the first study to use a probabilistic method to study blood flow. Section 2 briefly reviews EMA procedures and describes details of the velocity estimation method. Section 3 overviews erythrocyte tracking and speed estimation results based on the clinical data. In this section, the automated tracking performance achieved by the proposed method is compared with trained human grader tracking results. Section 4 concludes the study and suggests improvements, applications, and prospective future directions.

2. Methods

2.1. EMA procedures

The detailed procedure of EMA was described by Flower et al [22]. Figure 1 is a schematic diagram illustrating the EMA procedure. In summary, approximately 30 mL of blood were drawn 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.

 

Fig. 1 The procedure of Erythrocyte mediated angiography

Download Full Size | PDF

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.

2.2. 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.

2.3. 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.

2.4. Erythrocyte velocity estimation

2.4.1. 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 ${\left\{I_0^i\right\}}_{i=1}^T$. 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 ${\left\{I_0^i\right\}}_{i=1}^T$. Output images are denoted as ${\left\{I^i\right\}}_{i=1}^T$. Here T indicates the total number of frames in the video and $I_{\mathit{avg}}=\frac{{\sum }_{t=1}^T{I}_0^t}{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.

 

Fig. 2 Image preprocessing

Download Full Size | PDF

 

Fig. 3 The typical temporal average images for (a) 24.6fps video (b) 46.1fps video. The gray-scale value of images were stretched to the range of 0–255 to aid visual appreciation.

Download Full Size | PDF

2.4.2. 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.

 

Fig. 4 EMA-AV framework for erythrocyte velocity estimation.

Download Full Size | PDF

2.4.3. 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):

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):

2.4.4. 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: ${\left\{w_{k-1}^i,x_{k-1}^i\right\}}_{i=1}^N$, where $x_{k-1}^i$ is the state of particle i in time t_k−1, and $w_{k-1}^i$ is its corresponding weight (which must be normalized to describe the possibility of being a real object). The detailed definition of $w_k^i$ in our application is given in Section 2.4.6. The posterior PDF at time t_k−1 can be expressed as Eq. (4):

At time t_k, the posterior PDF can be approximated by a new set of weighted MC particles ${\left\{w_k^i,x_k^i\right\}}_{i=1}^N$, as shown in Eq. (5):

2.4.5. 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_{k-1}^i}$ to each MC particle, $x_{k-1}^i$, and transform the continuous integration into a discrete integration. $v_{x_{k-1}^i}$ is sampled from the uniform distribution in Eq. (7). A one-to-one mapping between the set ${\left\{w_{k-1}^i,x_{k-1}^i,v_{x_{k-1}^i}\right\}}_{i=1}^N$ and the set ${\left\{w_k^i,x_k^i\right\}}_{i=1}^N$ 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 ${\left\{w_{k-1}^i,x_{k-1}^i\right\}}_{i=1}^N$, 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_{k-1}^i$ was set as $\frac{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]. ${\overline{w}}_{k-1}^i$ 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 ${\left\{w_{k-1}^i,x_{k-1}^i\right\}}_{i=1}^{N_c}$, 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_{k-1}^i}$, 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 ${\left\{x_k^i\right\}}_{i=1}^{N_c}$, and their weights, ${\left\{w_k^i\right\}}_{i=1}^{N_c}$, 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_{k-1}^c$, 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/

3. Results and discussion

3.1. 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.

 

Fig. 5 Sample tracking results from the 46.1fps videos. (a) A venule is displayed as the target vessel. (b) An arteriole is displayed as the target vessel. The topmost images in (a) and (b) are variance maps for the corresponding video sequences, with the vessel skeleton labelled in red.

Download Full Size | PDF

 

Fig. 6 Sample tracking results from the 24.6fps videos. (a) A venule is displayed as the target vessel (b) An arteriole is displayed as the target vessel. The topmost images in (a) and (b) are variance maps for the corresponding video sequences, with the vessel skeleton labelled in red

Download Full Size | PDF

3.2. 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.

3.2.1. 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²) 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.

Table 1. Summary of statistical results analyzing mean erythrocyte velocities as determined by trained human graders and EMA-AV.

View Table | View all tables in this article

 

Fig. 7 Linear regression plots for mean erythrocytes velocities (for 11 vessels) between trained human graders and EMA-AV. (a) 24.6fps videos (b) 46.1fps videos

Download Full Size | PDF

 

Fig. 8 Bland-Altman plots for mean erythrocytes velocities (for 11 vessels) between trained human graders and EMA-AV. The 95% confidence intervals of the mean difference and limits of agreement are shown in the plot. (a) 24.6fps videos (b) 46.1fps videos

Download Full Size | PDF

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].

3.2.2. 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.

Table 2. Frame based ICC values of 11 vessels from 5 eyes in 24.6fps EMA videos. Under the vessel index column, E represents eye index, A represents arteriole index, and V represents venule index.

View Table | View all tables in this article

Table 3. Frame based ICC values of 11 vessels from 5 eyes in 46.1fps EMA videos. Under the vessel index column, E represents eye index, A represents arteriole index, and V represents venule index

View Table | View all tables in this article

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 2.1. In the 24.6fps videos, erythrocytes were more likely to aggregate and to split. For example, 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

Download Full Size | PDF

3.2.3. 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.

 

Fig. 10 Bland-Altman plots for average top 10% erythrocytes velocities (11 vessels) between trained graders and EMA-AV. The 95% confidence intervals of mean difference and limits of agreement are shown in the plot. (a) 24.6fps videos (b) 46.1fps videos

Download Full Size | PDF

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² value for mean erythrocyte velocity was 0.856, and the ICC value was 0.9553. For trained graders, the intra-setting R² 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].

4. 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² 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.