## 1.

## Introduction

Vascular dysregulation may play a role in the progression of diabetic retinopathy, age-related macular degeneration, and glaucoma.^{1}^{,}^{2} Whether the role of hemodynamics in disease pathogenesis is primary or secondary, quantitative measures of blood flow and volume represent potential biomarkers for ocular disease detection, progression, and response to treatment. Quantification of ocular hemodynamics has been historically challenging due to the lack of appropriate diagnostics. Optical methods such as laser speckle^{3} and laser Doppler imaging^{4} are limited in their quantitative capabilities, particularly under conditions of multiple scattering. Optical coherence tomography (OCT) has shown promise for quantification of hemodynamics due to its capability to isolate primarily single scattered light. Dynamic scattering, as quantified by OCT angiography,^{5}^{,}^{6} is related to red blood cell content. Doppler OCT may be used to calculate volumetric flow, and this process is facilitated by the method of *en face* integration,^{7} which obviates the cumbersome determination of vessel angle.^{8}^{,}^{9} Rapid two-dimensional scanning is typically required to generate volumetric data sets that assess pulsatility.^{10} However, rapid scanning sacrifices transverse sampling density and may lead to biased flow estimates. One alternative approach is to gate Doppler OCT acquisition to the heartbeat.^{11} Another alternative approach is to use a slow scanning protocol^{12} for assessing total flow. While slow scanning does not directly assess pulsatility, for assessment of nutritive blood flow, time-averaged flow may suffice. Here we argue that by choosing a slow scanning protocol that resamples a single vessel asynchronously with respect to the heartbeat, accurate time-averaged flow measurements are obtained. We show that total arterial and venous flow measurements in the rat retina agree, supporting the accuracy of this method.

## 2.

## Materials and Methods

The rat eye was imaged with a 1300-nm spectral/Fourier domain OCT microscope operating at a speed of 47,000 axial scans per second with a transverse resolution of 7.2 microns and an axial resolution of 3.6 microns (Ref. 5), using a contact lens and goniosol on the cornea. Male Sprague-Dawley rats (350 to 430 g) were anesthetized with 40 to $50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{kg}/\mathrm{h}$ alpha-chloralose, intravenously, maintained at 37°C with a homeothermic blanket and dilated with 1% tropicamide (Mydriacyl). Animals were ventilated with 80% air and 20% oxygen, and arterial blood pressure was monitored through a femoral artery cannula. All procedures were approved by the Institutional Animal Care and Use Committee, where these experiments were performed. We used an angiography method that performs high-pass filtering of the complex OCT signal along the slow axis of repeated cross-sectional images, after motion correction.^{5} Angiographic techniques that use the complex OCT signal,^{4} in contrast to techniques that use the magnitude^{13} or phase^{14} alone, can separate static and dynamic scattering with a linear filtering procedure. Performing filtering along the slow axis of repeated B-scans improves sensitivity to slow flow.^{5}^{,}^{15} The OCT angiography scanning protocol consisted of frames (cross-sectional images) of 512 axial scans repeated twice at each transverse y location,^{5} at a total of 5120 transverse y locations, over a field of view of $1.5\times 1.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{2}$. The interframe time was approximately 11 ms, sufficient for the dynamic scattering component of the complex signal in all vessels to decorrelate. As the scan was performed only twice at each $y$ location, high-pass filtering amounted to a complex subtraction. The total time required for OCT angiography was $\sim 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{min}$.

For displaying retinal angiograms as color *en face* images, volumetric angiographic data sets were averaged and downsampled to 1024 pixels along the slow ($y$) axis and upsampled to 1024 pixels along the fast ($x$) axis. The vitreoretinal interface was detected in all cross-sectional images using a simple automated adaptive thresholding technique. Two depth ranges, determined by offsets relative to the detected vitreoretinal interface, corresponding to the vitreal vasculature and the retinal capillary plexuses in the plexiform layers were determined [Fig. 1(a)]. This method is valid in the rat retina, since there are no rapid variations in inner retinal layer thicknesses, as are found in species with a macula. At each transverse location, the five largest voxel values in the angiographic data set within each range were summed to form images of both the vitreal vasculature and the capillary plexuses. These two images were combined by assigning the vitreal vasculature image to the red channel and the capillary plexuses image to the green channel. The color-coded inner retinal angiographic data are shown overlaid on the OCT intensity image [Fig. 1(c)]. A raster scan protocol was used to quantify blood flow. For our Doppler OCT scanning protocol, the B-scan rate was 10.0 Hz (accounting for the flyback time), and four volumes (256 B-scans with 4096 axial scans each) were acquired at intervals of $\sim 25$ to 30 s. Doppler shifts were calculated by estimating the frequency shift along the fast axis, before converting to velocity axial projections. The asynchronous condition was confirmed by analyzing the arterial blood pressure trace, sampled at the B-scan frequency.

The *en face* integration method presented earlier,^{7} when applied to a raster scanning protocol, does not explicitly account for pulsatility. A raster scanning protocol samples different locations in a vessel at different times with respect to the heartbeat. In order to make the discussion more conceptually tractable, we introduce the simplifying assumption that the overall vessel velocity profile shape does not vary during the heartbeat. Supporting evidence for this assumption can be found in a recent publication by Santisakultarm et al.^{16} If the velocity profile shape does not vary appreciably during the heartbeat, the measured velocity is a separable function of space and time. In other words, the velocity at any location in a vessel can be described as the product of a spatial function describing the velocity profile in that vessel and a temporal function p(t) describing pulsatility in that vessel, where $p(0)=1$:

## (1)

$${\mathrm{v}}_{\mathrm{z}}(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t})={\mathrm{v}}_{z}(\mathrm{x},\mathrm{y},\mathrm{z},0)\mathrm{p}(\mathrm{t}).$$It should be noted that the model above does not assume a particular velocity profile shape. The desired average flow is then given by the following expression:^{7}

## (2)

$${\mathrm{F}}_{\mathrm{avg}}=\int {\int}_{\mathrm{x},\mathrm{y}}\frac{1}{\mathrm{T}}[{\int}_{\mathrm{t}=0}^{\mathrm{T}}{\mathrm{v}}_{\mathrm{z}}(\mathrm{x},\mathrm{y},{\mathrm{z}}_{0},\mathrm{t})\mathrm{dt}]\mathrm{dxdy}=\int {\int}_{\mathrm{x},\mathrm{y}}{\mathrm{v}}_{\mathrm{z}}(\mathrm{x},\mathrm{y},{\mathrm{z}}_{0},0)\mathrm{dxdy}\times \overline{\mathrm{p}},$$## (3)

$$\overline{\mathrm{p}}=\frac{1}{\mathrm{T}}[{\int}_{\mathrm{t}=0}^{\mathrm{T}}\mathrm{p}(\mathrm{t})\mathrm{dt}]$$## (4)

$${\widehat{\mathrm{F}}}_{\mathrm{avg}}=\sum _{\mathrm{n}}[{\int}_{\mathrm{x}}{\mathrm{v}}_{\mathrm{z}}(x,{\mathrm{v}}_{\mathrm{y},\text{scan}}{\mathrm{t}}_{\mathrm{n}},{\mathrm{z}}_{0},{\mathrm{t}}_{\mathrm{n}})\mathrm{dx}]\times {\mathrm{v}}_{\mathrm{y},\text{scan}}{\mathrm{T}}_{\mathrm{s}}=\sum _{\mathrm{n}}\mathrm{g}({\mathrm{t}}_{\mathrm{n}})\times \mathrm{p}({\mathrm{t}}_{\mathrm{n}}).$$For convenience, we have defined

## (5)

$$\mathrm{g}(\mathrm{t})=[{\int}_{\mathrm{x}}{\mathrm{v}}_{\mathrm{z}}(\mathrm{x},{\mathrm{v}}_{\mathrm{y},\text{scan}}\mathrm{t},{\mathrm{z}}_{0},0)\mathrm{dx}]\times {\mathrm{v}}_{\mathrm{y},\text{scan}}{\mathrm{T}}_{\mathrm{s}}.$$The function g(t) represents line integrals of the velocity axial projection profile along the $x$ direction. The time-dependence arises from the fact that scanning is performed in the $y$ direction. Sampling a function in the time-domain leads to replication of the spectrum in the Fourier domain at intervals of the sampling frequency ${f}_{s}=1/{T}_{s}$. Moreover, multiplication leads to convolution in the Fourier domain and summation of samples leads to evaluation of the Fourier transform at zero frequency. In order to gain more insight into the above expression, we replace the summation of samples by the Fourier domain equivalent.

## (6)

$${\widehat{\mathrm{F}}}_{\mathrm{avg}}=\sum _{\mathrm{n}}\mathrm{g}({\mathrm{t}}_{\mathrm{n}})\times \mathrm{p}({\mathrm{t}}_{\mathrm{n}})=\frac{1}{{\mathrm{T}}_{\mathrm{s}}}\times [\sum _{\mathrm{m}}(\mathrm{G}*\mathrm{P})(\mathrm{f}-\frac{\mathrm{m}}{{\mathrm{T}}_{\mathrm{s}}}\left)\right]{|}_{\mathrm{f}=0}.$$In the above expression, G(f) and P(f) denote the Fourier transforms of g(t) and p(t), respectively, and * denotes convolution. P(f), the Fourier transform of the pulsatile function p(t), consists of a spectrum of repeated delta function spaced by ${f}_{0}=1/T$. Equation (6) is also closely related to the Poisson summation formula. We note that

## (7)

$$\mathrm{G}(0)={\int}_{\mathrm{t}}\mathrm{g}(\mathrm{t})\mathrm{dt}=\int {\int}_{\mathrm{x},\mathrm{t}}{\mathrm{v}}_{\mathrm{z}}({\mathrm{x},\mathrm{v}}_{\mathrm{y},\text{scan}}\mathrm{t},{\mathrm{z}}_{0},0){\mathrm{dxdt}\times \mathrm{v}}_{\mathrm{y},\text{scan}}{\mathrm{T}}_{\mathrm{s}}.$$By comparing Eqs. (2) and (6), it is obvious that in order for the flow estimate to be accurate, the following must hold:

## (8)

$${\widehat{\mathrm{F}}}_{\mathrm{avg}}\approx \frac{1}{{\mathrm{T}}_{\mathrm{s}}}\times \mathrm{G}(0)\times \overline{\mathrm{p}}=\int {\int}_{\mathrm{x},\mathrm{t}}{\mathrm{v}}_{\mathrm{z}}(\mathrm{x},{\mathrm{v}}_{\mathrm{y},\text{scan}}\mathrm{t},{\mathrm{z}}_{0},0)\mathrm{dxdt}\times {\mathrm{v}}_{\mathrm{y},\text{scan}}\times \overline{\mathrm{p}}={\mathrm{F}}_{\mathrm{avg}}.$$To understand the conditions under which the above approximate equality holds, we neglect the bandwidth of G(f) and consider only the effects of sampling the pulsatile waveform. This assumption is justified if sampling is sufficiently dense in the $y$ direction. (Allowing a non-negligible bandwidth for G(f) could account for sparse sampling density.) P(f) consists of a fundamental component (heartbeat frequency of ${f}_{0}$, determined from the arterial blood pressure waveform) and harmonics at $2{f}_{0}$, $3{f}_{0}$, etc. We assume that the highest harmonic with significant energy is equal to $N{f}_{0}$. This spectrum is replicated at intervals of ${f}_{s}=1/{T}_{s}$ due to the sampling operation. Thus the Nyquist criterion requires that the sampling frequency ${f}_{s}>2N{f}_{0}$. Here we argue that more generally, provided that no harmonic of the replicated pulsatile spectrum falls at zero frequency, the flow estimate is approximately equal to the actual time average flow, as expressed in Eq. (8). Thus the more general condition is that the sampling frequency is asynchronous with respect to the pulsatile frequency ($k{f}_{s}\ne n{f}_{0}$ for all harmonics $n=1,2,\dots N$ and positive integers $k$). If this condition is met, even though the pulsatile waveform may be undersampled and thus cannot be reconstructed, the time-average can be determined by averaging sampled values. The time-domain equivalent of this criterion is that different phases of the periodic heartbeat are all adequately sampled.

Flow in each vessel was calculated individually using the *en face* integration procedure described above. Only one location (*en face* plane) for each artery or vein branch was chosen for flow calculation. Locations were chosen with a sufficiently small angle to avoid phase wrapping and fringe washout artifacts, but with a sufficiently large vessel angle to yield a measurable Doppler shift.

## 3.

## Results and Discussion

Figure 1(a) shows a projection of the angiographic data set over 100 microns in the transverse direction. As shown in Fig. 1(a), the inner retinal vasculature is divided into three distinct layers (labeled 1, 2, and 3). Vitreal arteries and veins supply and drain the inner retina (1). Dense capillary beds are found in both the inner plexiform layer (2) and outer plexiform layer (3), with the latter being the most dense. This finding is consistent with the high oxygen consumption in the inner and outer plexiform layers of the inner retina, and the high ATP requirements of synaptic activity and maintenance of membrane ion gradients.^{17} Our angiographic imaging technique enables overlay of the vasculature [Fig. 1(a)] on the standard intensity image [Fig. 1(b)], coloring the capillary plexus green and the vitreal vessels red [Fig. 1(c)]. The inner plexiform layer vasculature [labeled 2 in Fig. 1(c)] corresponds to sublamin a, a layer responsible for the off-center ganglion cell contacts. Figure 1(d) shows peripheral vasculature in the retina, using the coloring scheme described above. In large vitreal vessels, multiply scattered photons are detected from the range of path lengths assigned to the capillary plexus and colored green. Thus, the multiple scattered photons (green channel) combine with the single scattered photons (red channel) to make the larger vitreal vessels appear yellow-orange, whereas they should be red based only on their anatomical location. Figure 1(e) shows wide-field angiographic imaging of the optic nerve head. The optic nerve head is at the bottom of the image.

Sampling of a pulsatile function must adhere to the Nyquist criterion in order to accurately reconstruct the pulsatile waveform. However, time-averaged flow measurements do not require accurate reconstruction of the pulsatile waveform. Our argument above suggests that accurate average flow measurements can be obtained, in the absence of pulsatility information, by ensuring asynchronous sampling with respect to the heartbeat. Figure 2(a) shows an arterial blood pressure trace (blue, oversampled at 1 kHz) along with sampled values at 10.0 Hz (red circles). The Fourier transform magnitude of a 25-s time course, an approximation of P(f), shows that the fundamental frequency ${f}_{0}$ is 5.9 Hz, and the highest frequency in the blood pressure signal is $\sim 59\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{Hz}$ [Fig. 2(b)]. Figure 2(c) shows a plot of the mean of the sampled blood pressures at different sampling rates. The positive and negative deviations of the mean of samples from the true mean in Fig. 2(c) occur when the asynchronous sampling condition is not met. At ${f}_{S}=10.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{Hz}$, the mean of the samples was 98.3 mm Hg, while the true mean was 98.5 mm Hg, and the asynchronous condition is met at this sampling frequency. Slight heart rate variability would broaden the delta functions in P(f) and, thus, further restrict the range of sampling frequencies where portions of the replicated spectrum do not fall at zero frequency. At the other extreme, a highly variable heart rate would satisfy the asynchronous condition independent of sampling rate as different phases of the heartbeat would automatically be sampled.

Figure 2(d) shows Doppler OCT velocity axial projections overlaid on arteries and veins in the grayscale OCT angiogram, while Fig. 2(e) shows a histogram of flows in arteries and veins. The venous flow is into the optic nerve head (positive $z$ direction) and hence is positive, while the arterial flow is out of the optic nerve head (negative $z$ direction) and hence is negative. Some residual line artifacts can be seen in arteries in Fig. 2(d) due to pulsatility. Due to sufficiently asynchronous sampling, we assume that the flow measurements presented in Fig. 2 can be taken to represent the temporal average of the flow. The magnitude of total retinal arterial blood flow was $6.48\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{L}/\mathrm{min}$ (standard deviation over four volumes of $0.72\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{L}/\mathrm{min}$), while the magnitude of total retinal venous blood flow was $6.16\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{L}/\mathrm{min}$ (standard deviation over four volumes of $1.15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{L}/\mathrm{min}$). The mean values are comparable to the time-average of previous pulsatile blood flow measurements in the rat, performed under isoflurane anesthesia.^{10} In steady state, the total retinal arterial and venous blood flows are expected to be equal. The fact that arterial and venous flows are approximately equal supports the accuracy of our method.

## 4.

## Conclusion

High-resolution angiography and total average blood flow measurements in the rat retina were demonstrated using OCT. Combined with asynchronous sampling, the method of *en face* integration was shown to enable robust, angle-independent, and simple time-averaged flow quantification. Sectoral blood flow has advantages over total retinal blood flow,^{10}^{,}^{18} including, for instance, that it is possible to compare sectoral flow to sectoral nerve fiber loss or visual fields in glaucoma.^{19} With measurements of oxygen saturation in individual retinal vessels, it will be possible to determine sectoral oxygen extraction, which, together with sectoral flow, will help to determine the local metabolic rate of oxygen consumption in the retina.

## Acknowledgments

We acknowledge support from the National Institutes of Health (R00NS067050, R01EB001954), the American Heart Association (IRG5440002), and the Glaucoma Research Foundation Catalyst for a Cure 2. We also thank Maria Angela Franceschini as well as Bruce Rosen for general support.