1. Introduction

Retinal oxygen metabolism is a key factor in the pathogenesis of leading blinding diseases, including diabetic retinopathy (DR) and age related macular degeneration (AMD) [1–4]. There are growing evidences suggesting that functional variations precede structural alterations [2, 4–8]. Thus, retinal metabolic rate of oxygen (rMRO₂), a hallmark of retinal metabolism, may be a key biomarker for early disease screening, evaluating disease severity, and monitoring response to therapeutic intervention. Two main parameters, retinal blood flow and retinal oxygenation, are required to quantify rMRO₂. Retinal blood flow can be calculated given the diameter of retinal blood vessels and their flow velocity, both of which can be reliably quantified using the state-of-the-art ophthalmic imaging tools [9–11]. Therefore, the major challenge that hinders the incorporation of rMRO₂ into ophthalmic care is the lack of a well-established method for accurate retinal oxygen measurement in clinics.

Though several existing technologies have been explored to quantify retinal oxygen, they are currently not suitable for measuring rMRO₂ during routine clinical visits. For example, the most accurate method in ophthalmic research uses intra-retinal microelectrodes to directly measure the oxygen tension in the retina [8, 12–14]. Unfortunately, the invasive nature of this procedure limits its practice, and only approximate measurements can be obtained at the retinal surface in human subjects during surgical procedures, such as vitrectomy [12]. Another disadvantage of microelectrode measurement is that it is nearly impossible to map the oxygen tension over a larger area in the retina. In order to identify spatial distribution of oxygenation in individual retinal blood vessels, imaging-based methods are required. For example, phosphorescence lifetime imaging can map the oxygen tension in mouse retinal vessels based on oxygen-dependent emission quenching of the injected phosphorescent probe [15, 16]. Unfortunately, the need to introduce toxic fluorescent probes into the systemic circulation makes it inappropriate for human uses. Magnetic resonance imaging (MRI) can detect dynamic retinal oxygenation response non-invasively in humans following inhalation challenge [17, 18]. However, MRI has low resolution and provides only qualitative measurement. Most importantly, this approach relies on dynamic contrast variation after oxygen challenge and cannot provide steady-state measurements [18].

Recently, researchers focused on image-based extraction of hemoglobin oxygen saturation (sO₂) rather than direct measurements of oxygen tension. sO₂ describes the percentage of hemoglobin binding sites occupied by oxygen molecules. Together with blood hemoglobin concentration, sO₂ can reveal the amount of oxygen carried in the blood, fulfilling the requirement for rMRO₂ calculation. Optical measurement of sO₂ is possible because the two forms of hemoglobin, oxy- and deoxy-hemoglobin (HbO₂ and HbR), have distinct optical properties. As the optical absorption coefficient of whole blood is the linear combination of these two forms of hemoglobin, multi-wavelength spectroscopy can be used to retrieve relative HbO₂-HbR concentration ratio and thus sO₂. Fundus photography [4–7], scanning laser ophthalmoscopy [19], and hyperspectral imaging [20, 21] have all been investigated to measure retinal circulation sO₂ in humans. These approaches rely on the backscattered photons from blood vessels and share similar imaging principle. The major drawback of these approaches is that they can only provide two dimensional en face images, which are susceptible to various interferences and noises in quantifying sO₂. Monte Carlo simulation of photon transport in the retina shows that quantified sO₂ are susceptible to (1) multiply scattered photons, (2) variations in blood vessel diameter, and (3) melanin concentration in the retinal pigment epithelium (RPE). All these factors collectively contribute up to 20% error in measuring sO₂ [22]. In contrast, photoacoustic ophthalmoscopy (PAOM), a three-dimensional (3D) imaging modality that is sensitive to optical absorption contrast, is robust at removing the influence from the background and multiple scattering. However, PAOM requires an ultrasound detector in direct contact with the eyelid [23–25] and no human retinal imaging has been demonstrated yet.

Visible light optical coherence tomography (vis-OCT) is an emerging non-invasive imaging modality suitable for retinal oximetry [26, 27]. Similar to traditional OCT, it uses a broadband light source to provide high-quality anatomical imaging in 3D based on low coherence interferometry. Instead of using near infrared (NIR) light as in all commercial OCT systems, vis-OCT takes advantage of visible light spectrum for improved resolution, and arguably enhanced scattering contrast from certain tissues [28]. Most importantly, sO₂-dependent hemoglobin optical absorption has the strongest spectral contrast within the visible-light wavelength range, which further facilitates sO₂ quantification [29]. If vis-OCT retinal oximetry can be successfully developed in humans, it will have the following advantages: first, 3D imaging and coherent detection minimize the influence from surrounding tissues to warrant high sO₂ quantification accuracy; second, vis-OCT has access to a more continuous optical absorption spectrum of blood, which facilitates inverse calculation since spectral fitting provides much higher sample density than discrete multi-wavelength measurements in other modalities; third, because conventional OCT is routinely used in ophthalmic clinics worldwide, adding functional imaging capability to OCT will enable a quick clinical adoption [30]; and fourth, OCT is a versatile technology which has already demonstrated accurate measurement of diameter and blood flow velocity within retinal vessels [9, 24]. Therefore, vis-OCT, providing both sO₂ and blood flow measurements, has the potential to quantify rMRO₂ by a single imaging modality, which has not been previously possible.

We are focusing on translating vis-OCT from benchtop to clinical applications. Previous studies have shown the capability of vis-OCT to quantify retinal sO₂ in rodents under both normal conditions and during oxygen challenge [27, 31, 32]. In addition, dynamic rMRO₂ measurement was also demonstrated [27]. Recently, we developed a human vis-OCT prototype, which retrieved structural retinal images with fine details and sub-micrometer axial resolution [33]. In this study, we seek to quantify sO₂ within central retinal circulations in healthy human volunteers. Our vis-OCT oximetry revealed that retinal arteries and veins in close proximity to the optic nerve head had at least 20% difference in sO₂. Moreover, we observed sO₂ variations in different levels of arterial branches. To overcome the low signal-to-noise ratio (SNR) in humans as compared with rodents due to reduced illumination power and stronger attenuation from larger eyeballs, we developed a statistical fitting-based approach for sampling OCT intensity. Then we showed that we could use this approach to obtain true OCT intensities in noisy images and discussed vis-OCT oximetry in details.

2. Methods

2.1 Hunan subject recruitment

Our studies were approved by the Northwestern University Institutional Review Board (IRB), and adhered to the tenets of the Declaration of Helsinki. All procedures took place in the Ophthalmology Department at the Northwestern Memorial Hospital. Healthy volunteers were recruited during their routine clinical visits. All recruited volunteers were provided informed consent before participation. Prior to imaging, all subjects underwent a complete eye examination performed by an ophthalmologist. Inclusion criteria required that the subject should have normal or near-normal vision. If the subject wore corrective lenses, they were eligible if the refractive error was in the range of + 2 to −6 diopters. Subjects with known retinal diseases, significant cataract, or vitreous cloudiness were excluded. Subjects underwent pupil dilation during their routine eye examination before vis-OCT imaging.

2.2 Data acquisition

The core of the human vis-OCT prototype is a Michelson interferometer using a supercontinuous laser (SuperK EXTREME EXW-6, NKT Photonics). The detailed specifications and capabilities of the system were described in details in our previous publication [33]. Briefly, the spectrum of the illumination light ranged from 506 nm to 621 nm, giving us the previously-characterized axial resolution of 0.97 µm in air. The square field-of-view (FOV) covered roughly 5 × 5 mm² on the retina. We switched between two raster scanning protocols, (1) 128 × 8 (A-lines × B-scans) during alignments and (2) 512 × 256 for acquisition, to optimize the tradeoff between imaging time and scanning density. The A-line rate was 25 kHz. In most imaging sessions, the laser intensity on the cornea was 220 µW. However, in rare cases, it was moderately increased to 280 µW to compensate for increased optical attenuation in older subjects. Even at the higher intensity, it was still well within the safety limit under the American National Standards Institute (ANSI) laser safety guidelines [34].

We imaged either one or both eyes of each volunteer. In the cases where only one eye was imaged, the right eye was preferred. We gently positioned the subject’s head on a motorized chin rest, which was translated using a remote control. During alignment, the subject was instructed to follow a fixation target using the fellow eye, while the operator monitored the real-time vis-OCT scans to locate a proper region-of-interest (ROI). Once we identified the desired ROI, we asked the subject to keep his/her eyes wide open and refrain from blinking during acquisition. Each acquisition took 5.2 seconds, and the entire session lasted less than five minutes per eye.

2.3 Vis-OCT oximetry

We applied our previously established spectroscopic algorithms to retrieve sO₂ [31]. Its accuracy and stability have been characterized using mathematical modeling, verified using ex vivo phantom experiments, and further confirmed during in vivo rodent studies [27, 35]. Briefly, we applied a series of short-time Fourier transforms (STFT) to split the entire OCT interferogram into multiple spectral sub-bands. Each sub-band was then individually reconstructed for its corresponding vis-OCT image stack. We manually assigned a ROI for each identified vessel along its longitudinal axis. We sampled its spectroscopic vis-OCT intensities from all sub-bands within a 25 µm thick slab located in the middle of the vessel. To inversely calculate sO₂, we employed a modified Beer-Lambert’s law described in our previous work [31] and established a linear relationship between the measured spectroscopic vis-OCT intensity and sO₂. Thus, we can retrieve blood sO₂ from the target vessel by least-squares fitting of the experimental data to the known attenuation spectra of oxygenated and deoxygenated whole blood.

2.4 Statistical sampling strategy on noisy vis-OCT images

Retrieving unbiased true vis-OCT intensity from the selected ROI is the key to accurate sO₂ calculation. Unfortunately, satisfying this requirement is often challenging in clinical settings, where SNR is usually sub-optimal due to safety and comfort limitations. The most straightforward way to cancel the effects of noise is to calculate the arithmetic mean value of a larger ROI. Unfortunately, the noises in OCT intensity images are not zero-mean due to the modulus operation applied as shown below, precluding the use of simple arithmetic average for obtaining unbiased estimation in low SNR situations.

Here we present a statistical fitting method to retrieve unbiased estimation of the true OCT intensity. This approach is based on the statistical distribution model of OCT intensity in noisy images. Our derivation assumes spectral-domain OCT (SD-OCT), but can be adapted for swept-source OCT (SS-OCT). Each spectrometer detector element recorded interferogram intensity D[k_i], where i indicate the sequence of the elements (i = 1, 2, …, N; N is the total number of elements); k_i is the corresponding optical wavenumber for the i^th element. Discrete Fourier transform translate the spectral domain interferogram into complex OCT A-line profiles d[x_i] in the spatial domain,

First, we must address the question that how stochastic noise translates from the spectral domain to the spatial domain. We base our derivation on a shot-noise limited OCT system, which is an accepted model to study SNR in OCT [36, 37]. Shot-noise originates from the discrete nature of photons, which arrive at the CCD detector at random interval during the exposure time. This stochastic process can be modeled by Poisson distribution. When the photon number is large, it approximates normal distribution. Thus, we write the general expression of spectrometer signal for each detector element as,

Apparently, d_s[x_i] is the noise-free A-line and d_n[x_i] is the additive noise term. First, we examine the OCT background noise when ${\text{d}}_{\text{s}}\left[x_i\right]=0$. Decomposing Eq. (3) using Euler’s formula leads to

For simplicity, we define $Re=\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\text{D}}_{\text{n}}\left[k_n\right]\cos \left(k_n{x}_i\right)$ and $Im=\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\text{D}}_{\text{n}}\left[k_n\right]\sin \left(k_n{x}_i\right)$. The OCT intensity can thus be expressed as the modulus of the complex profile

It becomes clear that the statistical distribution of OCT background intensity |d₀[x_i]| depends on the statistical distribution of Re and Im. Here, we only provide explicit analysis on Re, as the two terms are symmetrical and readers can apply the same approach to reach the conclusion for Im. Based on the notion that D_n[k_i] are independent and follow normal distribution, we have the expected value of Re

Here we introduce a scaler, $C=\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\cos }^2\left(k_n{x}_i\right)$, to simplify the expression to $\text{Var}\left\{Re\right\}=\frac{C}{N}{\text{σ}}^2$. When N is large, C is approximately 0.5. In addition to mean and variance calculations, we know that Re follows normal distribution, as it is the weighted sum of independent normally distributed random variables. Thus, we denote $Re~\text{Normal}\left(0, \frac{C}{N}{\text{σ}}^2\right)$. Similarly, we also have $Im~\text{Normal}\left(0, \frac{C}{N}{\text{σ}}^2\right)$. One can further prove that $Re$ and $Im$ are uncorrelated, as their covariance are zero due to the orthogonality of the trigonometric functions. In addition, as the vector pair $Re, Im$ satisfies the definition of multivariate normal distribution, it leads to the conclusion that that $Re$ and $Im$ are indeed two independent random variables, even though they share the same origin.

Upon inspection of Eq. (5), it becomes clear that OCT background intensity |d₀[x_i]| follows generalized chi distribution with two degrees of freedom, or Rayleigh distribution. The probability density function (pdf) of $\left|{\text{d}}_0\left[x_i\right]\right|=t$ is,

Now, we extend the same analysis chain toward the non-background area, i.e., the features in the OCT intensity image. Similarly, we decompose the noise-free A-line into real and imaginary parts, ${\text{d}}_{\text{s}}\left[x_i\right]=a+jb$, where a and b are both deterministic functions of x_i. $\left|{\text{d}}_s\left[x_i\right]\right|=\sqrt{{\text{a}}^2+b^2}$ is the desired true OCT intensity. We can re-write Eq. (4) and Eq. (5) as

It is again apparent that $\left(a+Re\right)$ and $\left(b-Im\right)$ are two normally distributed random variables and we have $\left(a+Re\right)~\text{Normal}\left(a, \frac{C}{N}{\text{σ}}^2\right)$ and $\left(b-Im\right)~\text{Normal}\left(b, \frac{C}{N}{\text{σ}}^2\right)$, respectively. Thus, |d[x_i]| follows non-central generalized chi distribution with two degrees of freedom, also known as Rice distribution. The pdf for $\left|{\text{d}}_0\left[x_i\right]\right|=t$ is

I₀ is the modified Bessel function of the first kind. We further introduce a scaler ${\epsilon }^2=\frac{C}{N}{\text{σ}}^2$ to simplify the expression. It has the same expression as the variance of the above Rice distribution. Another scaler $\nu =\sqrt{a^2+b^2}$ has the same expression as |d_s[x_i]|. Notably, ν is independent from the noise term σ. Thus, we have successfully separated OCT signal from noise, while the value of ν represents the true OCT intensity. For small signal $\nu \to 0$, $g\left(t\right)$ degenerates to $f\left(t\right)$ as expected. Thus, Eq. (11) can be used as a general expression to model the statistical distribution of OCT pixel intensity. By performing a statistical fitting on Eq. (11) using experiment results, we can determine ν and ε, thus retrieving the true OCT intensity.

Figure 1 illustrates the above theory using numerical simulation. We started from a sinusoidal interference spectrum with a length N = 2048 [Fig. 1(a)]. Then we introduced random white Gaussian noise additively to the interference spectrum. The standard deviation (S.D.) of the noise was chosen to produce an SNR of 30 dB (non-DC peak intensity vs. background S.D.) in the reconstructed OCT signal. The shaded area depicted ± S.D of the additive noise. After applying discrete Fourier transform on the simulated spectrum, we retrieved the complex value corresponding to the non-DC peak in the transformed result. It was subsequently plotted in the complex plane [Fig. 1(b)]. We repeatedly generated the spectra and the associated non-DC peak data points 10,000 times to reach a statistically stable solution. The entire process was analogous to sampling a homogenous region in a vis-OCT image at a specific depth, such as within a blood vessel. After all data were collected, we tested the spatial distribution of these complex points and found they follows bivariate normal distribution as discussed above. In addition, both the real and imaginary components followed univariate normal distribution. We then calculated the modules r of these points, which follows Rice distribution. Its pdf was depicted by the gray solid line in Fig. 1(c). For comparison, we conducted the entire simulation again using an interference spectrum with smaller modulation amplitude. The reconstructed signal had an SNR of 10 dB and the distribution of the modules r showed the unsymmetrical nature of the Rice distribution (the gray dashed line). In both cases, we performed statistical fitting using the obtained modules r distribution curves. The parameter ν from the fitting agreed with the constructed modulation amplitude of the simulated interference spectrum, verifying that ν was the unbiased estimation of true OCT intensity.

 

Fig. 1 Statistical characteristics of noise in OCT. (a) Simulated interferogram with added Gaussian noise from a single reflective surface (red solid line). The black solid line is the noise free interferogram. The shaded area depicts S.D. of the added noise, which corresponds to SNR = 30 dB for the reconstructed image; (b) The distribution of the corresponding non-DC peak values on the 2D complex plane after discrete Fourier transform reconstruction. The insert shows the spatial distribution of the complex values in the 2D plane, which were tested to follow bi-variate normal distribution. (c) The distribution of the modulus r. The solid line is generated using the data from (a) and (b). The dashed line showed a non-symmetrical case where SNR was much worse (10 dB).

Download Full Size | PDF

The flowchart in Fig. 2(a) summarizes the steps for retrieving true OCT intensity from human vis-OCT images. Following these steps, we subsequently tested our theory on a moderate-quality human vis-OCT image [Fig. 2(b)]. Here we define regional SNR (rSNR) as the ratio between the sampled mean OCT intensity and the background standard deviation within a select ROI. The latter was estimated using B-scans that don’t contain any identifiable anatomical structures at the given depth from the same volumetric scan. The vessel area (box 1) has rSNR of ~5 dB. We plotted the probability density of OCT intensities of three regions in Figs. 1(c)-1(e) in gray lines, respectively. It was then fitted using Rice distribution, and the corresponding theoretical curve were plotted on top of the experiment data in red lines. In all three cases, the coefficient of determination (R²) was at least 0.96, indicating a close fit. In Fig. 1(c), we observed minor discrepancies between the fitted OCT intensity from parameter ν (red arrow) and the arithmetically averaged mean (gray arrow). In both Fig. 1(d) and 1(e), the fitted OCT intensity (parameter ν) were close to zero (< 0.01) as expected, because both regions were backgrounds without any features. However, arithmetic means were biased (gray arrows) due to asymmetry of the modulus operation.

 

Fig. 2 (a) Flowchart showing the steps of retrieving true OCT intensity using statistical fitting approach. For OCT oximetry, the same ROI is used for all narrow-band images. (b) A typical human vis-OCT B-scan image showing three sampling regions: one from within a blood vessel (1) and two from the background at different depths (2) and (3). (c) The pdf of OCT intensity (gray line) from box (1) and its statistical fitting (red line). The gray and red arrows indicate the arithmetic mean and fitted OCT intensity, respectively. (d) and (e) The pdfs of pure OCT noise intensities from the background at their respective depths in gray lines (boxes 2 and 3). Red lines are the corresponding statistical fittings. The gray arrows indicate arithmetic means. The fitted OCT intensities were both approximately zero.

Download Full Size | PDF

It is worth noting that even though we developed the model assuming that the added noise was independent among detector elements, it could also be applied to situations when added noise share certain correlations. Specifically, there existed cross-talk between adjacent detector elements as introduced by the limited optical resolution of the spectrometer optics. In addition, the supercontinuum laser source used in our vis-OCT system contained high levels of relative intensity noise (RIN), which had strong 1/f dependency. It led to the observation that A-line regions closer to the zero delay line had higher background noise. If simple arithmetic mean were used, it would appear that the region closer to the zero-delay line (box 2) had higher averaged OCT intensity than the one further away from the zero-delay line (box 3), which were demonstrated in Fig. 1(d) and Fig. 1(e) (gray arrows). However, our statistical-fitting model was unaffected by this influence. The differences in noise level manifest itself in terms of the variance of the distribution, or the width of the pdf curves plotted, which did not affect the estimated OCT intensity.

3. Results

3.1 Statistical sampling can reduce oximetry estimation error

To demonstrate how the aforementioned statistical sampling approach benefit vis-OCT oximetry, we carried out the following numerical simulation. We repeatedly generated sinusoidal inteferogram spectrum with different levels of added Gaussian noises, yielding various SNR from 30 dB down to 0 dB. Then, we sampled the corresponding OCT intensity using both arithmetic mean (gray solid line) and statistical fitting (red solid line). The error in each measurement was interpreted as relative S.D. against the known preset true values. As shown in Fig. 3(a), when SNR was good (>20 dB), both methods had negligible deviations from the preset values. However, as SNR decreases, the estimation errors of arithmetic mean approach increased exponentially. In contrary, the statistical fitting approach provided more stable estimation in the 10 dB to 0 dB range.

 

Fig. 3 Comparison of OCT intensity sampling errors on the accuracy of sO₂ estimation. (a) The relative S.D. of extracting OCT intensities using arithmetic mean (black line) and statistical-fitting approach (red line) under different SNRs; (b) Numerical simulation result on sO₂ estimation error at different pre-set sO₂ level and relative S.D. of sampled OCT intensity. (sO₂: 100% = 1).

Download Full Size | PDF

Undoubtedly, the sampling error in the OCT intensity will propagate through the least-squares fitting process to affect the accuracy of sO₂ measurements. Unfortunately, due to the non-linearity of the non-negative fitting algorithm used, it is challenging to analytically formulate the direct correlation. However, we believe that a numerical simulation, though greatly simplified, can shed light on such connection and fulfill our analytical needs. Specifically, we simulated the transmission spectrum of blood at different sO₂ level (0.4 to 0.95) using Beer-Lambert’s law. For each sO₂ level, we perturbed the transmission spectrum by adding additive Gaussian noises. The relative S.D. of the noise was inferred from the findings in Fig. 3(a) and ranged from 0 to 0.3. We used non-negative least-squares fitting algorithm to retrieve sO₂ values from the distorted spectra. Figure 3(b) depicted the relationship among pre-set sO₂, relative errors in the sampled OCT intensity, and sO₂ estimation error. We empirically assumed that sO₂ estimation error of less than 5 percentage points (pp) is acceptable, as is to the left of the highlighted white solid line. Using the relationship depicted in Fig. 3(a) and the average spectroscopic vis-OCT image rSNR of 3-5 dB, we found that the statistical fitting approach shifted the estimation error region from the right dashed box to the left one, which roughly indicates a 3 pp reduction in sO₂ estimation error.

Though seemingly small, such 3 pp reduction in error is significant in clinical sO₂ measurements. First, spectroscopic OCT have even lower rSNR (2-3 dB drop) than the one shown in Fig. 2(b) due to fewer detector elements contributed to each narrow-band images [36]. Second, a 3 pp reduction in the measurement error is critical in disease management, as pathophysiological alterations in retina circulation sO₂ can be very small in early stages of blinding diseases such as diabetic retinopathy [8, 27, 38]. Reducing measurement error will lead to higher confidence level in identifying a statistical significant alteration, leading to higher diagnostic sensitivity in human studies.

3.2 Vis-OCT oximetry

Figure 4 shows the vis-OCT imaging results from four volunteers. Images in each row were obtained from different individuals’ eyes. Table 1 summarizes the basic information of the volunteers and the condition of the eyes.

 

Fig. 4 Vis-OCT imaging results with different rSNRs from four healthy volunteers (A to D). Column i: vis-OCT en face fundus images. Column ii: vis-OCT B-scan images from the corresponding locations highlighted by the dashed lines in column i. Column iii: OCT intensity spectrum taken from the identified central retinal arteries (A₁, black curve) and their corresponding least-squares fitted results (red line). Column iv: OCT intensity spectrum taken from the identified central retinal veins (V₁, black curve) and their corresponding least-squares fitted results (red line).

Download Full Size | PDF

Table 1. Summary of the imaged volunteers.

View Table | View all tables in this article

Column i of Fig. 4 are the en face vis-OCT fundus images. The FOV of the fundus images covers 6.4°, positioned within close proximity to the optic nerve head (ONH) to visualize the first- and second- order branches of the central retinal arteries and veins. At the current stage, the FOV selection around the ONH is intended to test sO₂ quantification at different vascular branching levels.

Vis-OCT images from the volunteers A and B have the highest qualities, registering peak rSNR of 12 dB and 11 dB within the blood vessels, respectively. The image from the volunteer C has moderate peak rSNR of 9.3 dB. Volunteer D’s image has the lowest peak rSNR of 8.4 dB, possibly due to extensive attenuation in the probing beam caused by age and minor cataract. In A and B, vis-OCT visualized arterioles and venules as small as 25 µm in diameter. In addition, the stripe patterns found in A and B correspond to nerve fiber layers (NFL) [33]. Two vis-OCT B-scan images, each averaged from 2 adjacent frames, are shown for each volunteer. The locations of the B-scan images are highlighted by the dashed-lines in column i. Appendix Fig. 6 provides a detailed labeling of the 13 identified layers from the vis-OCT B-scan A-ii.

We manually selected rectangular ROI’s within major artery and vein branches, which were indicated by the yellow dashed-boxes in column i. For each ROI, we followed the steps outlined in the methods section to retrieve true OCT intensity spectra. The obtained spectra were plotted in column iii and iv of Fig. 4 for identified first order arteries (A₁) and veins (V₁), respectively. Least-squares fitting on each individual spectrum yielded sO₂ of the correspondent retinal vessel segment. These fitted curves were plotted on top of the corresponding experiment spectra along with the measured sO₂ values.

In additional to extracting retinal sO₂ from selected ROI’s, we further verified that our statistical fitting method (SF) can reduce OCT oximetry measurement error against arithmetic mean approach (AM). For each vessel branches identified in Fig. 4, we randomly adjusted the ROI position and size by a small amount within the vessel segment, and then re-calculated sO₂ using both SF and AM. We repeated the process ten times and summarized the results in Table 2. The human experimental results suggested similar findings to our simulation results in section 3.1. First, sO₂ measurement error generally increased as rSNR worsened. Second, SF sO₂ results yield smaller errors than the corresponding AM values at all tested rSNR levels. Third, the largest measurement reduction for SF method appeared at lower rSNR values. Specifically, when rSNR fell below 8.7 dB, SF gave us a 5 pp reduction in the error.

Table 2. Reduction on sO₂ measurement error (Δσ) using statistical fitting (SF) against arithmetic mean (AM). sO₂ values are mean ± S.D of 10 measurements.

View Table | View all tables in this article

Our measurement showed that first order arteries had an sO₂ value at least 20 pp higher than the corresponding first order veins (p-value < 0.01). For the volunteers B and C, where multiple orders of retina vessel branches could be observed, we further plotted blood sO₂ against the branch orders. The results shown in Fig. 5 suggest that the oxygen content decreases as arteries bifurcate and travel downstream. Meanwhile, central retinal vein tends to have higher sO₂ values when it is closer to the ONH.

 

Fig. 5 Variations of sO₂ in major retinal artery and vein branches. The vessel branches are taken from ROI’s identified in the corresponding vis-OCT fundus images in Fig. 4. The arterial and venous sO₂ difference is statistically significant (**: p < 0.01).

Download Full Size | PDF

Though we could not find reference to this finding in previous literature on spatial distribution of retinal sO₂ [39, 40], these observations are consistent with Murray’s law [41]. Murray’s law established a predictive model between transmural oxygen gradient and perfusion pressure gradient, which suggested that higher branch-level vessels should have greater oxygen unloading rate, or lower sO₂ values. Furthermore, due to the proximity to the supplying retinal arteries and lower photoreceptor density, the region proximity to the OHN may demonstrate higher regional oxygen tension, and thus higher venular sO₂ as unutilized oxygen molecules diffuse back into the vein [40]. Nevertheless, the observed variations are small (4-5 pp), and are likely to be neglected if the measurements were not precise enough.

4. Discussion

In this work, we showed vis-OCT results from four healthy volunteers of both genders with an age span from 20s to 60s. We retrieved sO₂ in major central retinal veins and arteries following spectroscopic analysis. To increase the accuracy of sO₂ estimation when the SNR was low, we developed a statistical-fitting based sampling strategy. We verified the sampling strategy both theoretically and experimentally. In addition, both numerical simulation and human experimental data showed that the new approach reduces sO₂ estimation error by about 3 percentage points. However, there are still some caveats associated with the statistical model, as detailed in the following discussion.

First, we neglected speckle variations in our OCT intensity statistical model as we focused on low SNR situations. However, speckle variance is one of the major noise sources in higher quality OCT images. Speckle arises as the interference of diffusively reflected light. Unfortunately, the diffusive nature of the speckle mandates that its properties are difficult to describe using a simple formula, as it is affected by a combination of the illumination light and the microscopic structure of the sample. Schmitt et al explored the characteristics of speckle in OCT when imaging biological tissues [42]. They found that the OCT intensity from single polarized quasi-monochromatic illumination follows a pdf of exponential distribution. For unpolarized light, the incoherent addition of the two orthogonal components leads to a first order negative exponential distribution, which is related to Rayleigh distribution [42]. In our study, we obtained a second-order negative exponential distribution, which is also related to Rayleigh distribution. As a result, it is possible that the effects of speckle might be absorbed into our formula, due to the fact that we used a higher order model. More recently, Mahmud et al reported in a review article on OCT speckle variance that OCT intensity follows Gaussian distribution in solid tissues, while it changes to Rayleigh distribution in fluids [43]. It is worth noting that both these two distributions are the special cases of Rice distribution ($\epsilon \ll v$ leads to Gaussian distribution; and $v\ll \epsilon$ leads to Rayleigh distribution). The former situations are likely to be satisfied in high SNR images within solid tissue, and the later in fluids with high speckle variance, respectively. In such sense, our model might serve as a generalized description for OCT intensity variances.

Second, we noticed that our statistical model fails when blood flow reaches fringe washout velocity. In this situation, OCT intensity decreases significantly within the blood vessels. We observed a left-shifted peak and elevated tails when doing statistical fittings on regions with high blood flow. We can avoid this detrimental effect by limiting the ROIs to where the Doppler angle is large. As Doppler effect is most pronounced when the actual flow is parallel to the probing beam, limiting the Doppler angle close to 90° can minimize the perceived velocity by OCT detection so as not to reach the fringe washout threshold [44]. In addition, we consistently monitored the fitting process and rejected any fittings with R² less than 0.95. In case of rejected fittings, the ROI is moved to a different appropriate region where the established model held. Using the approach, we were able to extract and calculate sO₂ from various branches of central retinal arteries and veins. We also observed lower sO₂ with increasing orders of vascular branching. The observation is in accordance with Murray’s law on oxygen exchange in capillaries [41].

It is worth noting that our measured sO₂ of retinal veins are at the higher range compared to previously reported values [39, 40, 45, 46]. All these results were obtained by the less accurate multiple wavelength diffusive fundus imaging, which has intrinsic vulnerability due to lack of axial resolution and strong influences from vessels diameter and RPE pigmentation variations [22]. There are also a few other possible reasons that we cannot completely rule out currently. First, our sampling ROI’s are close to the root of ONH. It is possible that the higher oxygen tension around ONH leads to elevated sO₂ within veins [40]. However, due to the lack of a gold standard in retinal oxygenation measurements, it is premature to draw a definitive conclusion on this discrepancy. Second, there is a possibility that the visible light illumination itself during the vis-OCT alignment and acquisition may cause elevated venous sO₂ due to neural-vascular coupling. Although we believe that vis-OCT will unlikely to cause any detectable hemodynamic variations during the 5-s data acquisition period as compared with the light energy and illumination duration in flicker studies [47, 48], we will conduct further investigations. A dual-band OCT system, using both visible and NIR light [28] will be an ideal tool to examine this issue. Nevertheless, vis-OCT oximetry is a principally new technology and we still need extensive studies to establish its accuracy and reliability. In addition, an automated algorithm to assign ROI to retinal vessels, preferably based on image segmentation, is desirable to facilitate clinical applications in the future.

In summary we have demonstrated the feasibility of vis-OCT oximetry in humans. Using the statistical-fitting approach, we improved the accuracy of sO₂ extraction when the SNR is low, which can frequently be encountered in human imaging. Though data from only four healthy volunteers are presented here, we have imaged over 30 volunteers and are actively recruiting more to establish a normative database of retinal sO₂ values. Further improvement of the vis-OCT system will be focused on higher stability and better patient comfort.

Appendix

 

Fig. 6 A labeled vis-OCT B-scan showing imaged retinal layers. ILM: inner limiting membrane; NFL: nerve fiber layer; GCL: ganglion cell layer; IPL: inner plexiform layer; INL: inner nuclear layer; OPL: outer plexiform layer; ONL: outer nuclear layer; ELM: external limiting membrane; IS: photoreceptor inner segment; OS: photoreceptor outer segment; RPE: retinal pigment epithelium; BM: bruch's membrane.

Download Full Size | PDF

Funding

National Institute of Health (NIH) (R24EY022883, DP3DK108248, and R01EY026078); National Science Foundation (NSF) (CBET-1055379); H. F. Zhang has financial interests in Opticent Health.

References and links

1. D.-Y. Yu and S. J. Cringle, “Retinal degeneration and local oxygen metabolism,” Exp. Eye Res. 80(6), 745–751 (2005). [CrossRef]   [PubMed]  

2. R. SimóC. Hernández, and European Consortium for the Early Treatment of Diabetic Retinopathy (EUROCONDOR), “Neurodegeneration is an early event in diabetic retinopathy: therapeutic implications,” Br. J. Ophthalmol. 96(10), 1285–1290 (2012). [CrossRef]   [PubMed]  

3. T. Kurihara, P. D. Westenskow, M. L. Gantner, Y. Usui, A. Schultz, S. Bravo, E. Aguilar, C. Wittgrove, M. Sh. Friedlander, L. P. Paris, E. Chew, G. Siuzdak, and M. Friedlander, “Hypoxia-induced metabolic stress in retinal pigment epithelial cells is sufficient to induce photoreceptor degeneration,” eLife 5, e14319 (2016). [CrossRef]   [PubMed]  

4. S. H. Hardarson and E. Stefánsson, “Retinal oxygen saturation is altered in diabetic retinopathy,” Br. J. Ophthalmol. 96(4), 560–563 (2012). [CrossRef]   [PubMed]  

5. C. M. Jørgensen, S. H. Hardarson, and T. Bek, “The oxygen saturation in retinal vessels from diabetic patients depends on the severity and type of vision-threatening retinopathy,” Acta Ophthalmol. 92(1), 34–39 (2014). [CrossRef]   [PubMed]  

6. C. Jørgensen and T. Bek, “Increasing oxygen saturation in larger retinal vessels after photocoagulation for diabetic retinopathy,” Invest. Ophthalmol. Vis. Sci. 55(8), 5365–5369 (2014). [CrossRef]   [PubMed]  

7. M. Hammer, W. Vilser, T. Riemer, A. Mandecka, D. Schweitzer, U. Kühn, J. Dawczynski, F. Liemt, and J. Strobel, “Diabetic patients with retinopathy show increased retinal venous oxygen saturation,” Graefes Arch. Clin. Exp. Ophthalmol. 247(8), 1025–1030 (2009). [CrossRef]   [PubMed]  

8. J. C. Lau and R. A. Linsenmeier, “Increased intraretinal PO₂ in short-term diabetic rats,” Diabetes 63(12), 4338–4342 (2014). [CrossRef]   [PubMed]  

9. Y. Wang, B. A. Bower, J. A. Izatt, O. Tan, and D. Huang, “In vivo total retinal blood flow measurement by Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. 12(4), 041215 (2007). [CrossRef]   [PubMed]  

10. W. Soliman, M. Vinten, B. Sander, K. A. E.-N. Soliman, S. Yehya, M. S. A. Rahman, and M. Larsen, “Optical coherence tomography and vessel diameter changes after intravitreal bevacizumab in diabetic macular oedema,” Acta Ophthalmol. 86(4), 365–371 (2008). [CrossRef]   [PubMed]  

11. L. Pedersen, M. Grunkin, B. Ersbøll, K. Madsen, M. Larsen, N. Christoffersen, and U. Skands, “Quantitative measurement of changes in retinal vessel diameter in ocular fundus images,” Pattern Recognit. Lett. 21(13–14), 1215–1223 (2000). [CrossRef]  

12. E. Stefansson, M. B. Landers 3rd, and M. L. Wolbarsht, “Increased retinal oxygen supply following pan-retinal photocoagulation and vitrectomy and lensectomy,” Trans. Am. Ophthalmol. Soc. 79, 307–334 (1981). [PubMed]  

13. N. M. Holekamp, Y.-B. Shui, and D. Beebe, “Lower intraocular oxygen tension in diabetic patients: possible contribution to decreased incidence of nuclear sclerotic cataract,” Am. J. Ophthalmol. 141(6), 1027–1032 (2006). [CrossRef]   [PubMed]  

14. R. A. Linsenmeier and R. D. Braun, “Oxygen distribution and consumption in the cat retina during normoxia and hypoxemia,” J. Gen. Physiol. 99(2), 177–197 (1992). [CrossRef]   [PubMed]  

15. O. S. Finikova, A. Y. Lebedev, A. Aprelev, T. Troxler, F. Gao, C. Garnacho, S. Muro, R. M. Hochstrasser, and S. A. Vinogradov, “Oxygen microscopy by two-photon-excited phosphorescence,” ChemPhysChem 9(12), 1673–1679 (2008). [CrossRef]   [PubMed]  

16. R. D. Shonat and A. C. Kight, “Oxygen tension imaging in the mouse retina,” Ann. Biomed. Eng. 31(9), 1084–1096 (2003). [CrossRef]   [PubMed]  

17. B. A. Berkowitz, C. McDonald, Y. Ito, P. S. Tofts, Z. Latif, and J. Gross, “Measuring the human retinal oxygenation response to a hyperoxic challenge using MRI: eliminating blinking artifacts and demonstrating proof of concept,” Magn. Reson. Med. 46(2), 412–416 (2001). [CrossRef]   [PubMed]  

18. Y. Ito and B. A. Berkowitz, “MR studies of retinal oxygenation,” Vision Res. 41(10-11), 1307–1311 (2001). [CrossRef]   [PubMed]  

19. J. V. Kristjansdottir, S. H. Hardarson, G. H. Halldorsson, R. A. Karlsson, T. S. Eliasdottir, and E. Stefánsson, “Retinal oximetry with a scanning laser ophthalmoscope,” Invest. Ophthalmol. Vis. Sci. 55(5), 3120–3126 (2014). [CrossRef]   [PubMed]  

20. D. J. Mordant, I. Al-Abboud, G. Muyo, A. Gorman, A. R. Harvey, and A. I. McNaught, “Oxygen saturation measurements of the retinal vasculature in treated asymmetrical primary open-angle glaucoma using hyperspectral imaging,” Eye (Lond.) 28(10), 1190–1200 (2014). [CrossRef]   [PubMed]  

21. H. Li, W. Liu, B. Dong, J. V. Kaluzny, A. A. Fawzi, and H. F. Zhang, “Snapshot hyperspectral retinal imaging using compact spectral resolving detector array,” J. Biophotonics 2016, 00053 (2016). [CrossRef]  

22. W. Liu, S. Jiao, and H. F. Zhang, “Accuracy of retinal oximetry: a Monte Carlo investigation,” J. Biomed. Opt. 18(6), 066003 (2013). [CrossRef]   [PubMed]  

23. X. Shu, W. Liu, and H. F. Zhang, “Monte Carlo investigation on quantifying the retinal pigment epithelium melanin concentration by photoacoustic ophthalmoscopy,” J. Biomed. Opt. 20(10), 106005 (2015). [CrossRef]   [PubMed]  

24. W. Song, Q. Wei, W. Liu, T. Liu, J. Yi, N. Sheibani, A. A. Fawzi, R. A. Linsenmeier, S. Jiao, and H. F. Zhang, “A combined method to quantify the retinal metabolic rate of oxygen using photoacoustic ophthalmoscopy and optical coherence tomography,” Sci. Rep. 4, 6525 (2014). [CrossRef]   [PubMed]  

25. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]   [PubMed]  

26. S. P. Chong, C. W. Merkle, C. Leahy, H. Radhakrishnan, and V. J. Srinivasan, “Quantitative microvascular hemoglobin mapping using visible light spectroscopic Optical Coherence Tomography,” Biomed. Opt. Express 6(4), 1429–1450 (2015). [CrossRef]   [PubMed]  

27. J. Yi, W. Liu, S. Chen, V. Backman, N. Sheibani, C. M. Sorenson, A. A. Fawzi, R. A. Linsenmeier, and H. F. Zhang, “Visible light optical coherence tomography measures retinal oxygen metabolic response to systemic oxygenation,” Light Sci. Appl. 4(9), e334 (2015). [CrossRef]   [PubMed]  

28. S. Chen, X. Shu, J. Yi, A. Fawzi, and H. F. Zhang, “Dual-band optical coherence tomography using a single supercontinuum laser source,” J. Biomed. Opt. 21(6), 066013 (2016). [CrossRef]   [PubMed]  

29. N. Bosschaart, G. J. Edelman, M. C. G. Aalders, T. G. van Leeuwen, and D. J. Faber, “A literature review and novel theoretical approach on the optical properties of whole blood,” Lasers Med. Sci. 29(2), 453–479 (2014). [CrossRef]   [PubMed]  

30. M. Adhi and J. S. Duker, “Optical coherence tomography--current and future applications,” Curr. Opin. Ophthalmol. 24(3), 213–221 (2013). [CrossRef]   [PubMed]  

31. J. Yi, Q. Wei, W. Liu, V. Backman, and H. F. Zhang, “Visible-light optical coherence tomography for retinal oximetry,” Opt. Lett. 38(11), 1796–1798 (2013). [CrossRef]   [PubMed]  

32. B. T. Soetikno, J. Yi, R. Shah, W. Liu, P. Purta, H. F. Zhang, and A. A. Fawzi, “Inner retinal oxygen metabolism in the 50/10 oxygen-induced retinopathy model,” Sci. Rep. 5, 16752 (2015). [CrossRef]   [PubMed]  

33. J. Yi, S. Chen, X. Shu, A. A. Fawzi, and H. F. Zhang, “Human retinal imaging using visible-light optical coherence tomography guided by scanning laser ophthalmoscopy,” Biomed. Opt. Express 6(10), 3701–3713 (2015). [CrossRef]   [PubMed]  

34. F. C. Delori, R. H. Webb, D. H. Sliney, and American National Standards Institute, “Maximum permissible exposures for ocular safety (ANSI 2000), with emphasis on ophthalmic devices,” J. Opt. Soc. Am. A 24(5), 1250–1265 (2007). [CrossRef]   [PubMed]  

35. S. Chen, J. Yi, W. Liu, V. Backman, and H. F. Zhang, “Monte Carlo investigation of optical coherence tomography retinal oximetry,” IEEE Trans. Biomed. Eng. 62(9), 2308–2315 (2015). [CrossRef]   [PubMed]  

36. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef]   [PubMed]  

37. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). [CrossRef]   [PubMed]  

38. W. Liu, S. Wang, J. Yi, K. Zhang, S. Chen, C. M. Sorenson, N. Sheibani, and H. F. Zhang, “Increased inner retinal oxygen metabolism precedes microvascular alterations in rodent model with Type 1 diabetes,” Invest. Ophthalmol. Vis. Sci.in press.

39. D. J. Mordant, I. Al-Abboud, G. Muyo, A. Gorman, A. Sallam, P. Ritchie, A. R. Harvey, and A. I. McNaught, “Spectral imaging of the retina,” Eye (Lond.) 25(3), 309–320 (2011). [CrossRef]   [PubMed]  

40. A. M. Shahidi, S. R. Patel, J. G. Flanagan, and C. Hudson, “Regional variation in human retinal vessel oxygen saturation,” Exp. Eye Res. 113, 143–147 (2013). [CrossRef]   [PubMed]  

41. C. D. Murray, “The physiological principle of minimum work II. oxygen exchange in capillaries,” Proc. Natl. Acad. Sci. U.S.A. 12(5), 299–304 (1926). [CrossRef]   [PubMed]  

42. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]   [PubMed]  

43. M. S. Mahmud, D. W. Cadotte, B. Vuong, C. Sun, T. W. H. Luk, A. Mariampillai, and V. X. D. Yang, “Review of speckle and phase variance optical coherence tomography to visualize microvascular networks,” J. Biomed. Opt. 18(5), 050901 (2013). [CrossRef]   [PubMed]  

44. H. C. Hendargo, R. P. McNabb, A.-H. Dhalla, N. Shepherd, and J. A. Izatt, “Doppler velocity detection limitations in spectrometer-based versus swept-source optical coherence tomography,” Biomed. Opt. Express 2(8), 2175–2188 (2011). [CrossRef]   [PubMed]  

45. D. Schweitzer, M. Hammer, J. Kraft, E. Thamm, E. Königsdörffer, and J. Strobel, “In vivo measurement of the oxygen saturation of retinal vessels in healthy volunteers,” IEEE Trans. Biomed. Eng. 46(12), 1454–1465 (1999). [CrossRef]   [PubMed]  

46. S. Palkovits, M. Lasta, R. Told, D. Schmidl, A. Boltz, K. J. Napora, R. M. Werkmeister, A. Popa-Cherecheanu, G. Garhöfer, and L. Schmetterer, “Retinal oxygen metabolism during normoxia and hyperoxia in healthy subjects,” Invest. Ophthalmol. Vis. Sci. 55(8), 4707–4713 (2014). [CrossRef]   [PubMed]  

47. M. Hammer, W. Vilser, T. Riemer, F. Liemt, S. Jentsch, J. Dawczynski, and D. Schweitzer, “Retinal Venous Oxygen Saturation Increases by Flicker Light Stimulation,” Invest. Ophthalmol. Vis. Sci. 52(1), 274–277 (2011). [CrossRef]   [PubMed]  

48. K. Polak, L. Schmetterer, and C. E. Riva, “Influence of flicker frequency on flicker-induced changes of retinal vessel diameter,” Invest. Ophthalmol. Vis. Sci. 43(8), 2721–2726 (2002). [PubMed]