In vivo measurement of erythrocyte velocity and retinal blood flow using adaptive optics scanning laser ophthalmoscopy

In vivo measurement of retinal blood flow is obtained by measuring the blood velocity of erythrocytes and lumen diameters of the blood vessels using an adaptive optics scanning laser ophthalmoscope. Erythrocyte velocity is measured by tracking erythrocytes moving across a horizontal scanning line. This approach provides high temporal bandwidth measurements, allowing the fluctuation of blood flow during cardiac cycles to be measured. The technique is most applicable to medium-sized blood vessels.


Introduction
Retinal blood flow dynamics can be influenced by several retinal disorders, including diabetic retinopathy [1], glaucoma [2], and age-related macular degeneration [3]. Accurate evaluation of retinal blood flow has therefore been considered a potentially important requisite for in vivo assessment of retinal status as well as for monitoring pharmaceutical interventions in the eye. While various approaches have been developed for evaluating retinal blood flow in humans, including fluorescein angiography [4], laser Doppler velocimetry [5] and flowmetry [6], Doppler optical coherence tomography [7], stabilized Doppler flowmetry [8], and monitoring of leukocyte velocity [9,10], all of these techniques involve tradeoffs or limitations. For instance, fluorescein angiography requires injection of a contrast agent that is invasive. Techniques based on Doppler shift such as laser Doppler velocimetry and Doppler optical coherence tomography are limited to detection of blood flow in large retinal blood vessels because they measure blood flow by detecting the Doppler shift, which requires a signal large enough to provide reliable results [11]. Monitoring of leukocyte velocity works best for small capillaries near the fovea.
In the past few years, the use of the Adaptive Optics Scanning Laser Ophthalmoscope (AOSLO) has been proven to provide in vivo retinal images with superior lateral resolution [12,13]. These systems produce high contrast, high resolution retinal images by measuring the ocular aberrations with a wavefront sensor and correcting them with a deformable mirror providing near diffraction limited imaging of the human eye. Using the AOSLO, it has been shown that blood velocity can be measured in parafoveal capillaries by tracking leukocytes over time [9,14]. However, because leukocytes constitute less than 1% of the blood volume [15], leukocyte tracking-based blood velocity measurements are limited to capillaries in which the transit of blood cells is in a single file and the leukocyte velocity can represent the velocity of all blood cells [9]. The leukocytes are also flowing through the capillaries intermittently, and therefore this technique is temporally sparse, that is, it cannot measure blood velocity during periods when there are no leukocytes present.
In this study, we present an alternative use of AOSLO imaging that allows direct measurements of blood flow in medium-sized blood vessels by tracking the movement of erythrocytes across an imaging line. Because this technique is based on direct imaging of the light backscattered by erythrocytes, it does not require the use of any contrast dye. Since erythrocytes scatter over a wide range of wavelengths, we use a near infrared light source (840 nm central wavelength), which makes the imaging process comfortable for the subjects. We take advantage of the fact that, when focusing on a blood vessel, we directly visualize erythrocytes as moving bright dots. Figure 1 shows a movie of blood vessels under imaging conditions suitable for tracking both leukocytes and erythrocytes. Leukocytes are visible as large bright pulses, as has been described in detail by Martin and Roorda [9]. Separate from the intermittent large intensity variations arising from leukocytes, we can see much smaller objects continuously moving within the field, and these can be attributed to erythrocytes, which constitute 40% to 50% of the blood volume and appear as bright light-scattering dots. Since erythrocyte velocity is a good indicator of the general blood velocity, measurements of erythrocyte velocity should allow calculation of blood flow by taking into account the lumen diameter of the blood vessel.

The Indiana Adaptive Optics Scanning Laser Ophthalmoscope (AOSLO)
The blood flow measurements are conducted using the imaging system described in detail by Burns et al [13]. In brief, a low-coherence superluminescent diode (Superlum) with a central wavelength of 840 nm and a spectral full-width at half maximum (FWHM) of 50 nm is used as the light source. The pupil of the imaging beam is relayed to an 8 kHz resonant scanner (Electro-Optical Products Corporation) which generates a fast horizontal scan, and then to a mirror mounted on a galvanometer which supplies a programmable vertical scan. The scanned beam is then relayed to the eye, forming a raster pattern of 1.8 by 1.6 degrees angular extent on the retina. The light returning from the retina passes back through the system and is focused onto a retinal conjugate plane where one of eight interchangeable confocal aperture stops is located. An aperture wheel on a computer controlled XY stage (Thorlabs) is used to select the aperture stop and precisely position it. Light passing the aperture stop is imaged onto an RCA avalanche photodiode. The adaptive optics control of the system is maintained using a Shack-Hartmann sensor and a MEMS deformable mirror (Boston Micromachines) with 140 actuators and 4 μm stroke.

Scanning modes
In the regular mode of scanning, the horizontal scanning line is moved vertically and each image is generated on a line-by-line basis. The resulting images have two spatial dimensions (horizontal position and vertical position) and are referred as XY images. In this case, the imaging system forms retinal images at a rate of 15 Hz, which is too slow to reliably measure blood velocity. For this reason we manipulate the voltage driving the vertical scanner with a programmable D/A converter, allowing us to switch to a different scanning mode. In this mode, the scanning line is moved vertically during the first portion of each image frame, forming a partial XY image. The vertical movement is stopped when the scanning line reaches a retinal feature of interest, for example, a blood vessel. The image lines, generated from this point on until the end of the frame formation, consist of image information from the same retinal location retraced repeatedly by the horizontal sweep of the imaging beam at 8 kHz. Even though the scanning line does not actually move vertically during this time period, the frame grabber continues to generate the image line after line. As a result, the latter portion of each frame consists of image lines collected sequentially from the same retinal location. So, rather than having two spatial dimensions as XY images, this partial image has one spatial dimension (horizontal position) and one temporal dimension (acquisition sequence). We refer to the scanning mode with the vertical scan halted as XT scanning and the resulting partial images as XT images. When the whole frame is completed, the scanning line is returned to the vertical starting position and the whole procedure repeats. Schematics of the two different scanning modes can be seen in Fig. 2, together with examples of the resulting image frames.

Blood velocity calculation
For the example in Fig. 2, the video shows bright objects continuously moving with the blood flow from frame to frame. We believe most of these objects represent backscattering from moving erythrocytes, because they are more numerous than leukocytes, and larger than blood platelets. When the vertical scan is stopped over the blood vessel, the bright objects are visualized as diagonal streaks in the XT images. Figure 3 illustrates the formation of the streaks. In this specific case, an erythrocyte is moving from the upper right to the bottom left over time due to blood flow. During the time period required for the erythrocyte to cross the scanning line, it is continuously sampled by the horizontally sweeping imaging beam. Because the horizontal component of the velocity vector has a direction from right to the left, the erythrocyte moves more and more to the left on the sequentially obtained image lines, thus forming a diagonal streak in the XT image. The horizontal velocity component is simply the distance the erythrocyte travels horizontally, divided by the time this takes, i.e., the slope of the streak. After calculating the horizontal velocity component, it is straightforward to compute the erythrocyte velocity by accounting for the angle between the scanning line and the vessel axis.
In practice this computation requires several intermediate steps. First, we remove the sinusoidal distortion of the images caused by the resonant galvanometer with offline software developed in MATLAB (Mathworks). Next, we measure the angle θ (0 < θ < π/2) of the streaks within the central lumen of the vessel. Since we are imaging the blood flow from the top, streaks appearing within the central lumen can result from cells moving at different location in depth within the diameter of the blood vessel, which causes different slopes. However, since cells near the center of the blood vessel typically move with the highest velocity, these cells will produce streaks with the shallowest slopes at a given position within the vessel. We used the shallowest slopes to compute the blood flow velocity at each position. This is a relatively easy discrimination, because when the primary focus is in the middle of the blood vessel, the streaks from this plane are well focused and easy to identify. For each XT image, 3 measurements are obtained and then averaged, producing an estimate of the horizontal component of blood velocity along the central vessel axis: (1) where f is the horizontal scan frequency (8 kHz), and k is the system magnification calculated taking into account the axial length of the subject's eye. Horizontal velocity V p is then adjusted for the angle α between the horizontal scanning line and the vessel axis to calculate the blood velocity V ax along the axis of the blood vessel: (2)

Vessel lumen diameter measurement and flow rate calculation
The contrast of the blood vessels in a confocal SLO image depends on the interaction of light and the blood vessel. The light that comes back to the detector is a combination of the light scattered from the vessel wall and the backscattered light from the blood flowing inside the vessel lumen [16]. Various methods have been developed to quantify retinal blood vessel diameter [17][18][19]. Typically the blood vessel intensity profile is fitted with a Gaussian function and the vessel width is estimated as the width between two preset thresholds on the Gaussian curve. When computed in this manner, the vessel width does not account for the dark regions in the peripheral lumen, assuming these represent the vessel wall.
With our XT images, we can resolve the blood flow near the vessel wall in unprecedented detail, which allows us to detect cells moving within these dark regions. In our observations, cell movements in these regions are sparse, presumably because most of the erythrocytes move in the central lumen as a result of shear force, forming a plasma layer near the vessel wall [20]. The blood flow near the wall appears darker, because of the lack of light-scattering blood cells as well as the simple geometrical reason that we are imaging an en face representation of a cylindrical vessel from the top, and the total lumen volume sampled on the sides is small. However, with the sparse cell movements resolved as streaks in the XT image, we can better estimate the width of the vessel lumen based on the extreme horizontal positions over which we can see evidence of blood flow, which to our knowledge has never been possible before. As in Fig. 4, the oblique sectional diameter D p is obtained by measuring the distance between extreme horizontal positions of the streaks, and the lumen diameter D v is calculated using the following formula: (3) where α is the angle between the horizontal scanning line and the vessel axis.
Another problem with some measuring techniques is that vessels are obscured by a central bright stripe resulting from specular reflection of light from the vessel walls, which will contaminate data. While this can be controlled using polarization [19], the problem does not arise in our XT images, because stationary structures appear as vertical lines, and we can easily exclude the specular reflection from both the width and velocity estimation (Fig. 4).
Once velocity and vessel diameter are calculated we can calculate the total flow rate. To do this we make the assumption that velocities within the vessel are non-Newtonian. That is, the velocity profile in a blood vessel is blunter than a parabola which only applies to perfect laminar flow [21]. The velocity profile factor, where V ax is the blood velocity along the vessel axis and V s is the velocity averaged over the vessel cross section, is thus lower than the value of 2 expected for perfect parabolic flow. The velocity profile factor changes with the lumen diameter of the vessel. Koustisaris [17] introduced a profile factor function (4) where D v and D c are the blood lumen and erythrocyte diameters, respectively [22].
After applying the profile factor, the average blood volume flow Q is computed as: (5) Combining Eq. (1) through Eq. (5) allows us to compute the flow rate for individual vessels.
The blood flow fluctuates during the subject's cardiac cycle because of the heartbeat. Thus pulsation needs to be taken into account when measuring blood volume flow. In our case we take measurements at different phases of the cardiac cycle and average them across the cardiac cycle.

Human subjects
Three subjects were tested. Subjects were dilated with 0.5% Tropicamide. Pupil position was stabilized by a mouth rest and aligned relative to the imaging system with a three-axis positioning stage. When the vertical scan was running, the exposure of the imaging beam (380 μW, 850 nm) and the wavefront sensing beacon (50 μW, 680 nm) were more than 10 times below the ANSI safe exposure level for continuous viewing at this field size (650 μm by 550 μm). When the vertical scan was stopped, the same optical power was spread on a smaller field size (according to the calculations for the standard, 25 μm by 650 μm), but the exposure was still safe according to ANSI standards. The vertical scan motion was programmed to resume after no more than 3 seconds, which ensured safety and comfort. This research was approved by the Indiana University Institutional Review Board and met the requirements of the Declaration of Helsinki.

Results and discussion
A regular SLO image can be distorted if eye movements occur during image acquisition. For the velocity measurements presented here the role of eye movements is minimized in several ways. First, the measurement time (related to the number of selected image lines) required for a single velocity estimate is about 5 ms on average, although up to 10 ms is required to measure the slowest cell velocities. Second, each frame includes an XY scanning portion, and therefore we can estimate how much the eye has moved during the 1/15 second between successive XY scans. For the results presented here we used data where the movements between frames were limited to a few microns (no saccades). Finally, during the XT portion of each frame we can detect motion of the eye based on the fact that stationary structures should appear as vertical straight lines. If there is a horizontal eye movement, then these lines will deviate from the vertical. If there is a vertical component of motion, and the stationary structures are small, then the scan line will be moved off the structure and the vertical lines will be interrupted. We restricted our data analysis to measurements taken where clean vertical straight lines were present, and interframe motions were minimal. Figure 5 shows a movie of the blood flow within an artery of a healthy subject's retina, consisting of 40 frames acquired in 2.667 seconds. The artery is the vessel at point (a) in Fig.  7, whose lumen diameter is about 65 μm. The movie consists of XT images, along with an animation showing the velocity change. The blood velocity averaged over the vessel cross section is plotted over time, with the data points representing the calculated value from the consecutive XT images. The velocity fluctuates with the expected pulsatile pattern, rising rapidly at the beginning of a cardiac cycle and dropping slowly afterward. The velocity waveform shows a dicrotic notch after the systolic peak, which is very typical in arterial pressure waveforms and is thought to result from the closure of the aortic valve [23]. The fluctuation repeats 3 times in the 2.667 seconds, matching closely the subject's heartbeat period which was measured immediately after data acquisition. Figure 6 shows data from a retinal vessel, where the scanning line crosses the vessel axis at a very small angle. The left part of the XT image is taken at a region near the vessel wall, while the right part is at the center of the vessel. The corresponding calculated velocity profile is shown as a bar plot. The steeper slopes on the left show that cells near the vessel wall are moving more slowly than cells near the central lumen, as expected. The velocity profile can only be measured in blood vessels which cross the scanning line with a small angle, which depends on the blood vessel orientation and is not controlled in the current system. Although the cell velocity profile shown here is close to the left side of a blunted parabolic profile, we do not have sufficient data to test the flow model used to calculate flow in this paper.
To test the consistency of the measurement of blood flow using our method, we measured the total flow in the branches of an artery before and after a bifurcation, at least 5 vessel diameters away from the bifurcation to avoid non-laminar blood flow conditions. The results are displayed in Table 1. Figure 7 shows the location of the blood vessels. The flow rate was measured consecutively, first in the parent vessel and then in the daughter vessels. Each measurement took 5 minutes during which we assume that the actual average blood flow rate remained the same. The calculated time-averaged blood flow rate before the bifurcation at position (a) was 2.30 μl/min and the sum in the two daughter vessels at positions (b) and (c) after the bifurcation was 2.37 μl/min.
The system as currently implemented does have some limitations. First, our velocity measurements are required to be taken in the plane of the blood vessel. If a blood vessel were oriented at an angle in depth, we would underestimate the velocity. This is not a major problem, because most of the vessels in the size range of interest lie approximately in the plane of retina, with connections between vascular layers being primarily vertical and readily visible since they go out of focus [24]. The second limitation relates to the size of vessels we can measure. We have not been able to reliably measure flow in the smallest capillaries, due in part to the low contrast. Similarly for the largest retinal vessels, very high blood velocities are not readily measured due to frequency limitations of the horizontal scanner in our system (8 kHz). This scan frequency represents a fundamental sampling rate for this technique, just as it does for the similar approach used by the Heidelberg flowmeter [25]. As a result, this present system cannot accurately measure the highest blood velocities in blood vessels with diameters larger than 100 μm.
There are two approaches to addressing the limitation to the upper end of the velocity range of this technique. The first is simply to use a faster horizontal scanner. The second approach is to control the angle at which the scanning line crosses the blood vessel. The main limitation is that a flowing erythrocyte must take sufficient time to cross the scanning line so that there are enough samples to estimate the slope of the streaks. The time to cross the scanning line is inversely proportional to the component of the erythrocyte velocity perpendicular to the fast scanning direction. Thus, decreasing the angle α between the scanning line and the vessel axis will decrease this velocity component and increase the time during which a given erythrocyte is visible in the XT image and allow higher velocities to be measured. This also has the advantage that it would allow better measurements of cell velocity and the related velocity profile on more blood vessels with a more complete range of orientation angles, since the third limitation of the current system is that it cannot take accurate measurements in blood vessels which cross the scanning line at a large angle. This would require rotating the scan direction with respect to the blood vessel orientation to select an optimal angle for a given vessel, which is technically feasible using an image rotator or rotation of the galvanometer.

Conclusion
Our AOSLO system with programmable scan can image blood flow in retinal blood vessels noninvasively with high resolution. By tracking the movement of erythrocytes across a vertically stationary scanning line in an AOSLO image, the blood flow velocity can be measured. By including the geometry of the vessel the volumetric blood flow rate can be calculated. Our AOSLO blood flow imaging system has a very high spatial and temporal resolution and can accurately trace the known fluctuation of blood flow during the cardiac cycle. This technique has the potential to fill in the gap between blood flow measurement techniques such a laser Doppler velocimetry, which can measure large vessels, and leukocyte imaging, which can measure small parafoveal capillaries.    Determination of lumen diameter (D v ). D p is the oblique sectional diameter based on the estimation of the extreme horizontal positions of streaks caused by cells moving near the vessel wall, and α is the angle between the scanning line and the blood vessel orientation. Stationary reflections appear as vertical lines. Scale bar is 100 μm.   Locations of arterial blood flow measurements before (a) and after a bifurcation (b & c). Scale bar is 100 μm. Table 1 Calculated results of the vessel branches in Fig. 7 Position of the vessel a b c