Measuring oxygen saturation in retinal and choroidal circulations in rats using visible light optical coherence tomography angiography

: Visible light optical coherence tomography (vis-OCT) has demonstrated its capability of measuring vascular oxygen saturation (sO 2 ) in vivo . Enhanced by OCT angiography, the signal from microvasculature can be further isolated and directly used for sO 2 extraction. In this work, we extended the theoretical formulation for OCT angiography-based oximetry by incorporating the contribution from motion contrast enhancement. We presented a new method to eliminate the associated confounding variables due to blood flow. First, we performed in vitro experiments to verify our theory, showing a stable spectral derivative within the selected wavelength bands for sO 2 extraction. Then, we tested our method in vivo to measure retinal sO 2 in rats inhaling different gas mixtures: normal air, 5% CO 2 , pure O 2 , and 10% O 2 . Absolute sO 2 values in major arterioles and venules in the retinal circulation can be accurately measured. In addition, we demonstrated the relative changes of sO 2 can be measured non-invasively from choriocapillaris immediately underneath the retinal pigmented epithelium (RPE) in rodents.


Introduction
Optical coherence tomography (OCT) demonstrated tremendous success in imaging three dimensional (3D) tissue anatomy with micrometer-level spatial resolution on a macroscopic length scale [1]. In the past two decades, the technique has been translated into various clinical practices [2]. In particular, OCT has joined the standard of care in ophthalmology, providing rich 3D anatomical images from, for example, various retinal layers, choroid, and optic nerve head. Besides the technical advances on high-speed and high-resolution retinal imaging, functional imaging of various hemodynamic parameters has emerged as an invaluable supplement to anatomical OCT imaging. Perhaps the most important developments of functional OCT in retinal imaging are the measurement of retinal blood flow [3][4][5][6][7][8][9], and label-free retinal OCT microangiography (OMAG) [10][11][12]. Because the neuroretina has a high metabolic demand, the blood flow in the retinal circulation can be an excellent indicator of the viability of the retina and has been associated to various retinal diseases [13][14][15].
Another important hemodynamic parameter of the retinal circulation is the hemoglobin oxygenation saturation or sO2. The sO 2 difference between the arterial and venous blood vessels indicates oxygen extraction from the circulation. When combined with blood flow measurements, oxygen metabolism in the retina can be quantified by gas volume of oxygen consumed per unit time [16,17]. The attempts of using OCT to measure blood oxygenation dated back more than a decade ago, where the spectroscopic contrast around the isosbestic point around 800 nm from oxy-and deoxy-hemoglobin was characterized [18][19][20][21]. Since then, several groups also explored the oxygenation dependent scattering coefficient to create contrast for oxygenated or deoxygenated whole blood. However, none of the efforts has progressed into successful measurements of absolute sO2 in vivo. The reason is that the absorption contrast around 800 nm was overwhelmed by the wavelength-dependent scattering coefficient of red blood cells/whole blood, which have later been theoretically and numerically demonstrated [22,23]. Knowing that hemoglobin has much higher absorbing contrast within visible-light spectral range than within NIR range, Yi et. al. and Robles et. al. explored the feasibility of using visible-light OCT (vis-OCT) to quantify sO 2 from ex vivo and in vitro samples, respectively [24,25]. They further showed in vivo measurement of blood oxygenation [26], and the first demonstration of in vivo retinal oximetry in rodents using vis-OCT combing with the spatially-resolved spectroscopic analysis [27]. Later, to isolate the blood signal from the static tissue in OCT signal, we proposed to combine the speckle contrast and the spectroscopic contrast of blood simultaneously [28] to calculate sO 2 in microvasculature. We demonstrated that sO 2 could be directly calculated from wavelengthdependent OCT angiography without the cumbersome image processing previously needed to locate blood vessels. A similar strategy was recently adapted for both exposed cortical imaging and optical-power-cancelled retinal imaging in rodents [29]. However, the theoretical treatment did not consider the complication of blood flow to the angiography signal, which could confound the blood oxygenation calculation in vivo.
Here, we expanded our theoretical formulation by taking blood flow contrast into account in OCT angiography-based oximetry. Experimental verification was conducted by pumping anticoagulated whole bovine blood through a capillary tube. We tested the influence of different flow rates to the spectral derivative (i.e. the first-order partial derivative of spatially resolved OCT spectrum against wavelength). Further, we proposed a new data processing algorithm to normalize different flow rate and yield a robust sO 2 calculation. We tested our method in rat retina in vivo and detected sO 2 variation upon inhaling different gas mixtures. Moreover we demonstrated, for the first time, a non-invasive measurement of the relative sO 2 changes in choriocapillaris using vis-OCT.

Theory
In this section, we described the extended theoretical formulation for the spectroscopic analysis on vis-OCT angiography-based oximetry. The following motion-enhanced dynamic scattering model is to demonstrate that a linear relationship can be established between sO 2 and OCT angiographic derived spectrum measurements. Without loss of generality, we can start the formulation from a wavelength-dependent A-line signal from OCT angiogram ( ) , AA z λ , where z defines the axial coordinate and λ denotes the wavelength. Note that we will omit the finite axial resolution and sensitivity roll-off for simplicity purpose. At any given point along the A-line, where t μ is the total attenuation coefficient of local tissue attributed to both optical absorption and scattering. For whole blood, it is a linear combination of absorption coefficient a μ , scattering coefficient s μ , and a packing factor W, t a s W μ μ μ = + [22,32]. For other retinal tissue, reduced total attenuation coefficient, As we take the partial derivative of OD against wavelength (spectral derivative), we have If we assume that the wavelength is much smaller than the motion displacement of the blood cells, the phase variation becomes random and the enhancement factor ( ) G ϕ Δ loses wavelength dependence. The spectral derivative of ( ) G ϕ Δ could be considered to be zero and the equation becomes Now let us consider the attenuation contribution from the whole blood and the static tissues. As light propagates in tissue, the cumulative attenuation is  μ was close to zero [32,33]. As a result, the spectral derivative of Eq. (7) is We can further express _ Oxy It is clear that Eq. (11) establishes a linear relationship between the spectral derivative of the OD function and sO 2 , while ( ) serves as a constant scaling factor. During our experiment, l b is a constant determined by the choice of ROI. Its exact value can be estimated by measuring the mean optical path on the OCT angiograph.

Animal preparation
We anesthetized wild-type Long Evans rats using Ketamine/Xylazine cocktail solution (11.45 mg Ketamine and 1.71 mg Xylazine per milliliter of solution, respectively). The solution was administrated via intraperitoneal injection (IP) at a dosage of 87 mg Ketamine and 13 mg Xylazine per kilogram of body weight. A supplementary dose (30% of the original volume) was administrated after 30 minutes of the initial injection to maintain the animal at deep anesthesia. We applied 0.5% Tetracaine hydrochloride ophthalmic solution for local eye anesthesia and 1% Tropicamide ophthalmic solution for pupil dilation. During the entire period of anesthesia, we maintained the core body temperature at 37 °C using an electronic heating pad with feedback control (FHC Inc). The animal was held on a custom-made imaging stage and properly immobilized. A pulse oximeter attached to the rear limb of the animal provided real-time monitoring of its heart rate and systemic arterial oxygen saturation.
The animal was first ventilated with normal air (21% O 2 , 79% N 2 ) at 4.2 Standard Liter per Minute (SLPM). Then we supplied the following gas mixtures in the order of: (1) 5% carbon dioxide (21% O 2 , 74% N 2 , and 5% CO 2 ); (2) pure oxygen (100% O 2 ); and (3) 10% oxygen (10% O 2 and 90% N 2 ). For each gas inhalation, the animal was allowed to stabilize for at least 3 minutes before the measurement was taken. Additional time was allowed if the pulse oximeter reading showed large fluctuation (>2% deviation over 10 seconds). We supplied normal air for 3 minutes between changes. We took one additional measurement after the sequence with normal air ventilation. For each OCT measurement, the corresponding oximeter readings were also recorded simultaneously.
All the experimental procedures were approved by the Northwestern University IACUC and conformed to the Association for Research in Vision and Ophthalmology Statement on Animal Research.

Image acquisition
We developed a spectral-domain vis-OCT system to image rodent retina in vivo. A supercontinuum laser source (SuperK EXTREME EXW-6, NKT Photonics) provided broadband illumination. The output beam was delivered to the input arm of a free-space Michaelson interferometer via fiber optics. The beam splitter (CM1-BS013, Thorlabs) of the interferometer divided and redirected the beam to the reference arm and the sample arm. In the reference arm, a set of carefully chosen glass plates compensated for dispersions created by the unbalanced optical components. A continuously adjustable ND filter (NDC-50C-2M, Thorlabs) attenuated the beam intensity to enhance the dynamic range of the interference signal. The sample arm consists of a pair of scanning galvanometer mirror (Nutfield Technology) and a Keplerian telescope (0.2 × magnification) to project rectangular raster scanning pattern on the rat fundus. The spectrometer covered a spectrum range from 520 nm to 630 nm. We estimated the transverse and axial resolutions of our vis-OCT system to be 15 μm and 1.7 μm, respectively. A fast CMOS line scan camera (spL2048-140km, Basler) digitized the interference spectrum for image reconstruction. All acquisition sequence and recording were synchronized by a customized LabView (National Instrument) virtual instrument program. We adopted uni-directional scans in our OCT angiography. Here we denote the fast-scan direction as the x axis and slow scan direction as the y axis. During one B-scan, we translated the laser spot along the positive x axis at 50-kHz A-line rate and a total of 400 A-lines were recorded. Immediately after each B-scan, we returned (flied back) the laser spot to the x axis origin without acquiring any signal. The returning speed was limited to four times the forward scanning speed, which was constrained by the mechanical properties of the galvanometer. We repeated the process five times at each y location before moving to the next one. For one vis-OCT angiography image stack, we had 5 × 512 = 2560 B-scans at 512 different y locations, each B-scan consisted of 400 A-lines. The entire OCT image stack took 25.6 seconds to acquire.

Vis-OCT imaging of retinal microvasculature
We generated vis-OCT angiography using phase-sensitive decorrelation algorithm [11,28,34]. The contrast originated from the flowing blood cells inside vessels as they randomly change the phase and amplitude of back-scattered light. The random phase shift caused decorrelation in OCT signal and can be enhanced via high-pass filtering. However, such enhancement may be confounded by bulk motion, which also contributed to the signal via the same mechanism. These motion originated from normal physiological movements (e.g. heart beat and breathing) and the rotating galvanometer mirrors, both of which were inevitable. We must take measures to minimize their effects so as not to interfere with vis-OCT angiography. We adopted a three-step approach to correct the aforementioned bulk motion decorrelation within each y location: (1) Bulk image shift; (2) Axial global phase fluctuations (AGPF) correction; and (3) Lateral global phase fluctuations (LGPF) correction [35]. Without loss of generality, we used the notation ( , , ) A x z t to represent the recovered complex vis-OCT signal of the five consecutive B-scans at a given y location; x and z represents the lateral and axial location, respectively; t indicates the specific B-scan in time order, whose value ranges from 1 to 5 in our case. The first step was to co-register the five real-valued OCT structure image at each y location, so that motion shift greater than the imaging resolution can be corrected. We calculated the cross-correlation function between the two adjacent B-scan images (t + 1 and t). We identified the maximum value of the cross-correlation function along x and z directions, which indicates when the two images were properly aligned. The second image was then shifted back so the two images were properly registered.
The second step was to correct the global phase fluctuations with two phase modifiers AGPF and LGPF on the co-registered OCT images. The AGPF can be estimated by integrating the phase variation along the A-line axial direction (z axis) as while the LGPF can be estimated by integrating the phase variation along the x axis as We multiplied the two phase modifiers AGFP and LGFP to the corresponding complex OCT signal pixel-wise. The corrected complex valued vis-OCT signal is denoted The B-scan angiography cross-section ( ) , , AG x z t at location y can then be calculated by where  represents the expected value; and ⋅ represents the modulus operation. We iterated the process at each y location and generated a three dimensional image stack ( ) , , AG x y z .
At last, en face angiogram was generated by integrating the first 15 pixels of the highest motion contrast along z direction.

sO 2 estimation
We performed spectroscopic angiography to obtain the spectral profile of motion enhanced contrast of blood vessels [10,11,28], which was used to calculate the absolute blood sO 2 . Instead of using the full spectrum, we employed short-time Fourier transform (STFT) to perform spatial-spectral dual-domain analysis. We empirically chose four Gaussian-shaped STFT windows centered at 565 nm, 568 nm, 571 nm, and 574 nm, each with 10 nm fullwidth-half-maximum (FWHM) bandwidth. The full spectrum was windowed using these narrow-band STFT windows, respectively. For each center wavelength chosen, we generated the corresponding narrow-band OCT angiography using the method described in the above section. The calculated wavelength-dependent OCT angiography is denoted , OD z λ . We performed a linear regression to calculate the slope of the straight line that went through the center of the four OD points. We used this slope to estimate the partial derivative of the OCT angiograph spectrum against wavelength.
To convert the spectral derivative into absolute sO 2 values, we performed a calibration based on pulse oximetry spO 2 readings where we assumed that the major retinal arteries had the same sO 2 as systemic arterial sO 2 . Based on the calibration on each animal, we can measure the absolute sO 2 values in the major retinal arteries and venules around the optic disk.
We summarized the data processing steps in Table 1. Table 1. Summary of data processing procedures to extract sO2.
Operations procedure to extract sO 2 for any given blood vessel 1 To generate wavelength-dependent OCT angiography using STFT and motion enhanced dynamic scattering contrast according to Eq.  To demonstrate that the spectral derivative method is not affected by flow velocity, we designed and performed the following in vitro experiment. We fabricated a capillary phantom using a polystyrene tube (inner diameter: 125 µm) filled with anticoagulated bovine whole blood (Quad Five). We fully oxygenated the blood by exposing it to normal air at room temperature for 1 hour prior to the experiment. The phantom was embedded beneath 2.5% agarose gel to mimic a blood vessel in the retina. We controlled the blood flow rate in the tube using a syringe pump (PHD 2000 Infusion, Harvard Apparatus) at various pre-set values within physiological range. In order to image the in vitro sample, we slightly modified our OCT system. An objective lens (LSM03-VIS, Thorlabs) replaced the telescope system to create a focusing spot about 5.5 µm in diameter. The same data acquisition and processing protocols were applied. Figure 1(a) compares the cross-sectional view of OCT structural images and motionenhanced images. It can be clearly seen that surrounding agarose gel has strong scattering background in conventional OCT. In comparison, the motion-enhanced images largely suppressed the static background, while the moving blood within the tube was enhanced due to dynamic scattering. Even when the flow rate was set to zero, intrinsic random Brownian motion was sufficient to enhance the contrast (the left most images in Fig. 1(a)). On the sectioned en face view of the OCT angiograph showing the capillary phantom, we took four rectangular region of interest (ROI's) of equal size along its longitudinal axis ( Fig. 1(b)). We extracted the wavelength-dependent motion contrast spectrum within the selected ROI's to evaluate the influence of flow velocity. As expected, the shapes of these spectra were comparable, showing characteristic absorption spectrum of oxygenated whole blood. The spectrum calculated from zero-flow phantom is shown in Fig. 1(c).

Phantom experiment verification
We extracted two parameters to characterize the obtained motion contrast spectrum, the averaged intensity (I m ) and spectral derivative (∂I/∂λ) from 565 nm to 574 nm. Their values were plotted against increasing flow rate in Figs. 1(d) and 1(e), respectively. We performed liner regressions to estimate the slopes of both curves. I m showed a clear decreasing trend with a negative slope, which is statistically significant (p < 0.001). On the contrary, the slope of the ∂I/∂λ curve did not deviate from zero (p = 0.74). The result showed that ∂I/∂λ was flowindependent across the tested range, suggesting that sO 2 estimation based on ∂I/∂λ would not be affected by different flow rate.   Figure 2(a) is the conventional structural OCT image. The anatomical features of the retina can be clearly distinguished. The cross sectional OCT angiogram is shown in Fig. 2(b). We further overlaid the structure image with the angiography in Fig. 2(c), showing locations of each vessel within the retina. Besides the vasculatures in retinal circulation, we also found strong blood motion contrast directly beneath the RPE, which suggested the existence of blood flows in choriocapillaris. Beyond that, visible light attenuated rapidly. We then compared the en face view of rat fundus obtained by two OCT modality. Figure  3(a) illustrates the conventional structural OCT image. Although the superficial vessels can be clearly observed, we cannot distinguish the microvasculature. In comparison, Fig. 3(b) showed motion-enhanced dynamic scattering contrast of the same data set, where the hyperreflective signals from the static tissues are largely suppressed. Due to limited lateral resolution, choriocapillaris appeared as a single bright layer and was excluded through segmentation. Further, we calculated the sO 2 of the major retinal vessels in the retinal circulation during normal oxygen inhalation, as the color-encoded angiographic image in Fig.  3(b). We can clearly distinguish the arteries and veins based on their individual sO 2 . Finally, we rendered the 3D vascular structure in Fig. 3(c). Besides the dominating brunches of central retinal arteries and veins, laminated structures of micro-vasculatures within nerve fiber layer (NFL) and outer plexiform layer (OPL) are also distinguishable. The fly-through en face microangiogram is shown in the movie.

Absolute sO 2 measurement in retinal circulation
We next performed in vivo verification of our spectral derivative method. We calibrated the absolute sO 2 values from the major retinal arteries against a pulse oximeter readings under different ventilation conditions (normal air, 5% CO 2 , pure O 2 , and 10% O 2 ), as shown in Fig.  4(a). The calibration process corresponded to the Step 5 in Table 1. The result shows a strong positive correlation (R 2 = 0.87) between the two independent measurements, indicating that the calibration offered reasonably good estimation of the true sO 2 . Based on the established calibration curve, we measured the major retinal arterial and venous sO 2 in respond to the change of inhalation gases. In the course of the experiment, we changed the inhalation gases in the order of normal air, 5% CO 2 , pure oxygen, and 10% oxygen. As shown in Fig. 4(b), the arterial and venous sO 2 were significantly different (p < 0.05), except when pure oxygen was inhaled and the shunting between arterioles and venules elevated the venous sO 2 . Besides the arteriovenous sO 2 difference, we observed a moderate increase (p = 0.06) in both arterial and venous sO 2 when we supplied 5% CO 2 gas. This increased sO 2 is expected because 5% CO 2 has the effect of vasodilation, which improves tissue oxygenation [36,37]. Both the arterial and venous sO 2 increased significantly (p < 0.05) under pure O 2 and decreased significantly (p<0.01) under 10% O 2 .

Relative sO 2 change in choriocapillaris
Given the fact that we had strong flow-contrast enhanced signals from choriocapillaris directly underneath the RPE, we attempted to measure the relative changes of sO 2 in choriocapillaris. Due to limited lateral resolution, we could not distinguish individual choriocapillaris; instead, we averaged the signal over the entire area of choriocapillaris. We performed image segmentation to remove any area shadowed by superficial blood vessels, and yielded one averaged spectrum ( ) , c AA z λ , where z c was the axial position in the choroid. In order to eliminate the contamination from the retinal circulation in the inner retina, we normalized the choriocapillaris OCT angiograph spectrum by the intensity spectrum from the layer directly above RPE and calculated the spectral derivative. We used the same calibration curve obtained from retinal circulation to obtain the relative change of sO 2 with respect to normal air breathing, as in Fig. 5. When we supplied 5% CO 2 gas mixture, the average sO 2 in choriocapillaris was observed to increase moderately by 2% (p = 0.29). When we supplied pure O 2 and 10% O 2 gas mixture, the average sO 2 was observed to increase by 3.5% (p < 0.001) and decreased by 5% (p < 0.001), respectively.

Discussion
In this work, we presented an extended theoretical formulation for OCT angiography-based oximetry. This method used the flow contrast in OCT angiography, based on which we normalized the wavelength-dependent angiography and provided a robust measurement of sO 2 . The approach differs from our previous vis-OCT oximetry publication [27], where only static OCT amplitude information was used. The introduction of dynamic contrast eliminates the need for cumbersome vessel segmentation, and has arguably higher sensitivity and noise resistance than our previous method. More importantly, the new approach also enables us to measure the relative change of sO 2 in chororiocapillaries, which are generally inaccessible in structural vis-OCT images. We observed the change of sO 2 in both retinal circulation and choriocapillaris under different inhaling gas mixtures, which include normal air, 5% CO 2 , 100% O 2 , and 10% O 2 . Our in vivo measurements were verified by independent and concurrent pulse oximetry.
In our dynamic scattering model, we assumed that the motion enhancement term (G) is wavelength-independent. The assumption requires that the phase fluctuation is random between 0 to 2π without any preference. The condition is automatically satisfied when the displacement of scattering particles during the scanning dwell time are much larger than the wavelength. In our case, the red blood cell travels about 0.13 mm during the 0.01 s dwell time, which suffices the requirement.
We calculated the first derivative of the spectral range from 560 nm to 580 nm to establish its correlation with sO 2 . It differs from the least square fitting method over much larger wavelength range (530 nm-600 nm) as previously reported by our and other groups [28,29]. The major benefit is that the permitted incidence power in the current study could be reduced by about 2/3, which is ~0.3mW in the current experiments. By removing unnecessary exposure, the safety margin for vis-OCT for human can be significantly improved and the discomfort can be dramatically reduced. By sacrificing the spectral information from the absorption contrast, our method requires independent oxygenation measurements for calibration purpose. Fortunately, pulse oximetry utilizes the pulsatile profile and provides the reading of the overall arterial oxygen saturation, which should be equal to the values in the major retinal arterioles immediately exiting optic nerve head. This allows us to calibrate our spectral derivatives as we showed in Fig. 4(a) in every experiment. Such a calibration is easily achievable using a pulse oximetry that is widely available in the clinical settings.
Interestingly, we were also able to measure sO 2 variations in choriocapillaris right beneath the RPE under different inhalation conditions. The choroidal circulation is very different from retinal circulation, both anatomically and functionally [38]. The choroidal blood is supplied through posterior ciliary arterioles and drained through vortex vein. To our concern, the major difference between the two circulations is that the choroidal circulation has very small arteriovenous sO 2 difference, whereas the retinal circulation has a large arteriovenous sO 2 difference [39,40]. This results in an almost uniform sO 2 level in the choriocapillaris under normal physiology conditions, which allows us to measure sO 2 over a relatively large area of choriocapillaris.
However, there are a few complications that prevent us providing absolute sO 2 measurements in choriocapillaris. First, the true sO 2 level in choriocapillaris is still unknown. The most accurate reference is the direct measurement of oxygen tension using an invasive oxygen sensitive microelectrode [38]. Although it is commonly assumed that the choroidal circulation has the systemic oxygenation as in the arterioles, the microelectrode measurements turned out to be lower than expected [38]. Without a reliable reference, we were not able to calibrate our spectral derivative as we did for the retinal vessels. Alternatively, one might consider using the least square fitting method with the full spectral range of the hemoglobin absorption contrast to characterize choriocapillaris sO 2 . However, the optical properties of dispersed single blood cells in capillary network are not entirely clear, which would bring another uncertainty in the least square fitting method. Second, light passes through entire retina before choriocapillaris. The accumulated attenuation is also wavelength dependent, which would affect the spectral derivative. Nonetheless, we found that this is the first time a non-invasive measurement can be implemented to measure relative choriocapillaris sO 2 variations.

Conclusion
In summary, we presented a new theoretical formulation for vis-OCT angiography-based oximetry. The flow contrast was used and characterized by in vitro experiments. After proper spectral normalization, we demonstrated, for the first time, robust in vivo measurements of absolute sO 2 variation in the retinal circulation and relative sO 2 variation in the choriocapillaris. The new method requires much reduced illumination power, which greatly benefits clinical applications of vis-OCT.