1. Introduction

Oxygen is the only molecule as the primary biological oxidant, and it is also an important substance to maintain cell activity [1]. Since oxygen cannot be stored in tissues, a constant and adequate supply must be guaranteed [2]. As one of the most metabolically active tissues, the retina consumes oxygen faster than the other tissues [3]. This phenomenon means that subtle changes in oxygen kinetics, such as oxygen delivery ($\mathrm {DO_2}$), oxygen metabolism ($\mathrm {MO_2}$), and oxygen extraction fraction (OEF), probably cause retinal dysfunction. In clinics, measuring retinal oxygen kinetics is of paramount importance for understanding metabolism and pathophysiology in ophthalmic diseases, diabetes [4,5], hypertension [5], and Alzheimer’s disease [6] as well. In previous studies, retinal oxygen saturation ($\mathrm {SO_2}$) is widely used to estimate retinal oxygen kinetics [7–9]. However, it is not sufficient to use $\mathrm {SO_2}$ alone since it only indicates the percentages of oxyhemoglobin in the blood. For oxygen kinetics analysis, not only functional information, such as $\mathrm {SO_2}$, and blood flow velocity (BFV), but also structural information, such as lumen diameter (LD), are required [10]. Although there is a crucial need for retinal oxygen kinetics analysis in clinics, no current technologies are available to fuse above mentioned structural-functional imaging and provide all these parameters simultaneously or quasi-simultaneously in a single system. In addition, the dynamic measurement technology with high spatial resolution can help researchers to study retinal diseases more comprehensively, but the existing technology is unable to achieve high spatial resolution and dynamic measurement simultaneously.

Previous studies have proven that separate measurements of retinal oxygen concentration ($\mathrm {C_{O2}}$), BFV, and retinal vessel diameter are feasible. For $\mathrm {C_{O2}}$ measurement, the oxygen tension measured by microelectrodes is the gold standard in clinics [11], but the invasive procedure is challenging in practice. As an alternative, phosphorescence imaging is another way for oxygen tension measurement [10,12]. Unfortunately, contrast agents used in phosphorescence and fluorescent imaging still have potential risks. In recent years, $\mathrm {SO_2}$ measured by dual-wavelength spectroscopic reflectometry has become the preferred metric for $\mathrm {C_{O2}}$ measurement because of its advantages in conducting noninvasive and wide field imaging. In addition, BFV and retinal vessel diameters are also available using existing techniques such as ultrasound Doppler imaging, laser Doppler flowmetry, optical coherence tomography (OCT), and photoacoustic Doppler velocimetry [13–16]. Based on these independent techniques, researchers have attempted to perform oxygen kinetics analysis. Shahidi et al. [17] achieved noninvasive retinal oxygen kinetics analysis in healthy subjects and sickle cell retinopathy patients based on $\mathrm {SO_2}$ and Doppler OCT. Aref et al. [18] used the same method for oxygen kinetics analysis in glaucomatous visual field loss patients. Fondi et al. [14] used a custom-built bi-directional Doppler OCT with dual-wavelength spectroscopic reflectometry for retinal OEF analysis and discovered significant differences between diabetic patients and healthy control. Unfortunately, the mismatched field of view (FOV) in dual-wavelength imaging and OCT technologies limits the final FOV of oxygen kinetics analysis [19]. More importantly, there are risks when multiple devices are used to measure SO2, BFV, and vessel diameter, independently. Firstly, different optical designs will introduce more problems such as phase diversity, defocus, or magnification to the imaging results. Secondly, different devices have different FOV and resolution, which greatly increases the difficulty of registration and image processing algorithms. Finally, switching between different devices will not only bring discomfort to patients but also increase operational complexity and time. Yi et al. reported a $\mathrm {SO_2}$ measurement technology based on the visible light OCT [20] and further demonstrated the accurate and robust measurement of retinal oxygen metabolism in rat eyes [21]. However, the factors such as FOV (20$^{\circ }$), sampling points (256 $\times$ 256 pixels), discomfort with exposure to visible light, and high cost, confine the popularity of the technology. Furthermore, limited by temporal resolution, existing research cannot achieve dynamical oxygen kinetics analysis. Recently, Hosseinaee et al. [22] showed the possibility of photoacoustic imaging technology for retinal $\mathrm {SO_2}$ measurement but have not yet carried out oxygen kinetics analysis.

In addition to $\mathrm {C_{O2}}$ and BFV, retinal vessel diameter is also one of the important parameters in oxygen kinetic analysis. However, whether the diameter of retinal vessels denotes the external diameter (ED) or the LD is confused. Previous studies have shown that retinal vascular metrics including ED, LD, wall thickness (WT) and wall-to-lumen ratio (WLR) are valid biomarkers of cerebrovascular and cardiovascular risk [23]. Vessel diameter can be easily measured by fundus imaging, but measuring WT, WLR, and other vascular metrics is challenging. In 2007, Harazny et al. [24] proposed a novel method to measure LD, WT, and WLR based on scanning laser Doppler flowmetry. This method measured the ED in the reflection image and the LD in the perfusion image. Based on this method, Baleanu et al. [25] evaluated WLR in patients with cerebrovascular injury, the results suggested that WLR may reflect different pathophysiological changes or different cerebrovascular disease stages of the retinal arterioles. Rizzoni et al. [26] further carried out statistical analysis for normal and hypertension groups, proving the effectiveness of this method. Although this method offered noninvasive and convenient measurement of vascular metrics, it has fallen out of favor because of its small FOV (an arteriole with a size between 80 and 140 $\mu$m), and limited spatial resolution (256 points $\times$ 64 lines $\times$ 128 lines). At present, optics scanning laser ophthalmoscopy or scanning laser Doppler flowmetry is the main technique for vascular metrics analysis [23]. However, these techniques are based on complex technology and procedures, and lack additional functional information for oxygen kinetics analysis, which restrict their use in daily clinical routine.

Here, we propose a noninvasive retinal oxygen kinetics imaging and analysis (ROKIA) technology with high spatial resolution and dynamic imaging by integrating dual-wavelength imaging with laser speckle contrast imaging (LSCI) technology in a common optical path. Based on dual-wavelength imaging and LSCI technologies, ROKIA technology can analyze the $\mathrm {SO_2}$ and BFV, independently. By fusing these technologies, ROKIA technology can measure the ED, LD of blood vessels in a large FOV (45$^{\circ }$), which is essential for retinal blood flow(RBF) calculation and oxygen kinetics analysis. In addition, WT and WLR are also available in ROKIA technology, which are effective biomarkers of cerebrovascular disease. Finally, a mathematical model is developed to translate $\mathrm {SO_2}$, RBF, and LD metrics into retinal oxygen kinetics metrics. This is the first demonstration of the vascular metrics (including ED, LD, WT, and WLR), blood flow metrics (including BFV and RBF), oxygenation metrics (including $\mathrm {SO_2}$, $\mathrm {C_{O2}}$, $\mathrm {DO_2}$, $\mathrm {MO_2}$, and OEF) via a single system, to the best of our knowledge. We further demonstrate the significant differences between non-proliferative diabetic retinopathy (NPDR) patients and healthy controls, as proof of concept.

2. Materials and methods

In this section, the experimental device is introduced first. Then, measuring principles of vascular metrics, blood flow metrics, and oxygenation metrics are introduced, followed by a detailed ROKIA workflow, and finally, the subjects are described.

2.1 Optical implementation

The instrumentation of ROKIA is shown in Fig. 1. For $\mathrm {SO_2}$ analysis, central wavelengths at 548 $nm$ and 605 $nm$ LEDs were selected (Huashang Light & Electronics Inc.). The bandwidth was 20 $nm$ (FWHM). For retinal BFV analysis, the wavelength and bandwidth of the laser source (SPL TR85, OSRAM GmbH) were 850 $nm$ and 3 $nm$ (FWHM), respectively. Different wavelengths were controlled by the light source driver. In order to achieve fast light source switching, the light source driver was controlled by electronic switches through a computer. The illumination is delivered using a custom-made fiber bundle (Chunhui Inc., Nanjing, China). This bundle simplifies optical design by avoiding the need for multiple dichroic beam splitters and bulk optics. At the input end, each fiber bundle is 2.5 mm in diameter, consisting of 2000 multimode fibers. Each fiber has a core diameter of 50 $\mu$m, a 64-66 $\mu$ m cladding diameter and a numerical aperture of 0.62. Each wavelength of light is coupled to a bundle of fibers (2000 multimode fibers) via a coupling lens. At the output end, the fiber bundles of all these three wavelengths are uniformly distributed in a ring shape (6000 multimode fibers), so that the light of all wavelengths can be uniformly radiated on the retina. Successively, the light beam irradiated the retina through the annular fiber bundle, converging lens, reflecting mirror, hollow mirror, and objective lens. The accessible emission limit of retinal exposure was measured to be $2.87 \times 10^{-4}$ $W$, which was less than the maximum limit of $8.034 \times 10^{-4}$ $W$ in the International Standard ISO 15004-2:2007. The unweighted corneal and lenticular infrared radiation irradiance of LEDs were also satisfied with ISO 15004-2:2007 which were $1.5\times 10^{-5}$ $mW/cm^{2}$ and $1.9\times 10^{-5}$ $mW/cm^{2}$, respectively. The reflected light beam passed through the objective lens, hollow mirror, and focus lens successively. Finally, a $12$ bit CMOS camera (M3ST507M-H, DO3THINK, Shenzhen, Guangdong, China) with $2730$ $\times$ $2048$ pixels and $3.45$ $\mathrm {\mu } \mathrm {m}$ pixel size was adopted to acquire images. The magnification of the acquisition system was $0.63$, while the f-number was $7.36$. The FOV is $45^{\circ }$. Each pixel in retina images is $6.59\mu m$. The optical design was performed using Zemax (Radiant Zemax LLC, Redmond, Washington) and the optical system was constructed with customized lenses. For more details on the system and calibrating methods please refer to our previous publication [27].

 

Fig. 1. Instrumentation of ROKIA.

Download Full Size | PDF

2.2 Correction of retinal motion artifacts

The LSCI technique is highly sensitive to motion artifacts. Motion artifacts not only reduce the spatial resolution of perfusion images but also cause errors in blood velocity measurement. For brain microvasculature imaging, researchers have proposed several LSCI registration algorithms [28–30] that effectively improve the quality of perfusion images. However, only low-power near-infrared light is acceptable in retinal imaging, resulting in a further decline in the quality of perfusion images. In addition, the complex and variable eye movements (such as nystagmus, saccade, etc.) in retinal LSCI imaging greatly increase the difficulty of registration. In this section, we proposed a stabilized laser speckle angiography method for retinal motion artifacts suppression. The pipeline of the stabilized laser speckle angiography is shown in Fig. 2. Firstly, the movement of optic disc is used to divide speckle pattern images (SPI) into three states: blinking, weak nystagmus and violent nystagmus. Since the intensity of the optic disc area in SPI is significantly higher than that of the vascular and tissue areas, weak and violent nystagmus can be effectively distinguished by projecting the maximum intensity of the SPI along the X-axis and Y-axis, as shown by the red arrow in Fig. 3. The frames between the red arrow are a group of weak nystagmus. Based on Fig. 3(a) and (b), the intensity curves of X-axis and Y-axis can be obtained by summing in pixels axes. And, we multiplied the X-axis intensity curve and Y-axis intensity curve. Then, the algorithm based on wavelet transform is used to detect all eye movement mutation points. All the mutation points are arranged according to the intensity from large to small and recorded as $\vec {P}$, and the index and intensity of $\vec {P}$ are recorded as ${\vec {P}}_{ind}$ and ${\vec {P}}_{int}$, respectively. A higher ${\vec {P}}_{int}$ means stronger eye movement at this mutation point. Finally, select the first N (N=10) points in $\vec {P}$ and sort them according to ${\vec {P}}_{ind}(N)$ from small to large to obtain N weak nystagmus groups (SPI marked by different colors in Fig. 2). The number of frames in each group is uncertain which depends on the subject. For each weak nystagmus group, the intra-group enhancement is performed to obtain a perfusion image (including rigid registration, laser speckle contrast (LSC) calculation and mean intensity projection). Finally, the inter-group enhancement is performed to obtain the final perfusion image. The standard steps of the registration algorithm include feature extraction, rigid registration, and non-rigid transformation. Generally, non-rigid registration algorithms have higher registration accuracy. However, limited by the low SNR of SPI images, the non-rigid registration algorithm may cause more errors, so the rigid registration algorithm (only translation and rotation) is used in the intra-group enhancement steps. In this work, the iterative closest point algorithm is used in the intra-group enhancement step, and the enhanced correlation coefficient (ECC) algorithm [31] is used in the inter-group step. Fig. 4 shows retina perfusion images obtained by different algorithms. Fig. 4(a) shows a typical retina perfusion image with motion artifacts. Fig. 4(b) shows the perfusion image that directly uses registration algorithms. Figure 4(b) is obtained by the following steps: first, SPI data are registered, and then the perfusion images are calculated using the LSC algorithm. Fig. 4(C) shows the perfusion image obtained by the algorithm proposed in this section. It can be seen that the algorithm proposed in this section effectively suppresses motion artifacts and provides high -quality retinal perfusion images, which provide important guarantees for subsequent analysis.

 

Fig. 2. The pipeline of stabilized laser speckle angiography.

Download Full Size | PDF

 

Fig. 3. The maximum intensity projection of SPI along X-axis and Y-axis. (a) X-axis; (b) Y-axis.

Download Full Size | PDF

 

Fig. 4. Retina perfusion images obtained by different algorithms. (a) A typical retina perfusion image with motion artifact; (b) The perfusion image that directly uses registration algorithms; (c) The perfusion image obtained by the method proposed in this section.

Download Full Size | PDF

2.3 Measuring principle of vascular metrics

This paper uses perfusion images to measure LD and 548 nm images to calculate ED. Since oxygenated hemoglobin and deoxygenated hemoglobin have the almost same absorption at the wavelength of 548 nm, retinal arteries and veins show similar gray values in 548 nm images. It is worth noting that the diameter of the blood vessel obtained in the 548 nm image is the ED (including the thickness of the blood vessel wall). The laser speckle imaging technology takes the moving blood cells as the measurement object, and the obtained blood vessel segmentation result can be regarded as LD. Therefore, the measurement of WT and WLR can be achieved by calculating the segmentation results of the blood vessels in the two modalities, as shown in Eq. (1), (2).

The pipeline of vascular metrics measurement is shown in Fig. 5. Firstly, 548 nm retinal images and SPI were collected, during which subjects fixate on the same target to minimize registration requirements. For SPI data, the stabilized laser speckle angiography algorithm was used to calculate the blood perfusion image. Next, the multi-modal registration algorithm was used to register the image [32]. After registration, the 548 nm image and retinal perfusion image were segmented. Retinal vessel segmentation is a challenging task due to the complex structure of retinal vessels, the inconspicuous contrast between vessels and background, and the uneven intensity distribution caused by tissue reflections and scattering. In recent years, deep learning has shown remarkable success in medical image analysis. Therefore, in this work, the U-Net based deep learning technology is used for segmentation [33]. The architecture of the original U-Net consists of a dual path (including a contracting path and expansive path), with 23 convolutional layers in total. We followed the standard unchanged architecture in the previous research [33], except that the input size is 512x512. Our deep learning segmentation pipeline consists of a pre-training stage with the public high-resolution fundus (HRF) dataset [34] and fine-tuning stage with our private multispectral fundus vessel images dataset. To be specific, we initially extracted the green channel of fundus images in the HRF dataset and resized the images into 512x512. At the pre-training stage, 30 image pairs (fundus image and vessel map) were utilized to pre-train U-Net and the average accuracy of the test set (the remaining 15 fundus images) was 94.3816% $\pm$ 0.7144%. During the fine-tuning stage, the pre-trained network was utilized as initialized weights and 10 multispectral fundus image pairs (with 548 nm wavelength) labeled by an ophthalmologist were employed to finetune the U-Net model. Finally, the average accuracy on the test set, including additional 20 multispectral fundus images, was 99.3983% $\pm$ 0.9238%. In addition, to enhance the information on vascular structure and reduce the difficulty of segmentation, we preprocessed the retinal reflection (548 nm) images and perfusion images, including the Contrast Limited Adaptive Histogram Equalization (CLAHE) [35], Gamma correction, and Hessian matrix vascular enhancement. After segmentation, the 548 nm image and perfusion image were recorded as ${\rm Ves}_{548}$ and $\rm {Ves}_{lsc}$ respectively. Then, we binarized ${\rm Ves}_{548}$ and $\rm {Ves}_{lsc}$ and extracted centerlines. The traditional hessian based on centerline extraction is used in this work [36]. We traversed each point of centerlines to calculate the vascular angle and used the Gaussian function to fit the cross-section curve of the vessel, using the full width at half maximum (FWHM) as the vessel diameter.

 

Fig. 5. The pipeline of vascular metrics measurement. Abbreviations: SLSA:stabilized laser speckle angiography; PI,perfusion image; ED,external diameter; LD, lumen diameter.

Download Full Size | PDF

Figure 6 shows the comparison of vascular cross-sections in 548 nm images and in perfusion images. Figure 6(a) shows the 548 nm image, the yellow solid line box and the blue solid line box mark the representative venous and arterial regions, respectively; Fig. 6(b) shows an enlarged view of the venous region in the 548 nm image and perfusion image, in which the red line marks the vein in the 548 nm image, and the green line marks the vein in the perfusion image; Fig. 6(c) shows the cross-sections of retinal veins, the red curve is the cross-section of the vein in the 548 nm image, and the green curve is the cross-section of the vein in the perfusion image; Fig. 6(d) shows the perfusion image, the yellow dotted box and the blue dotted box mark the representative venous and arterial regions, respectively; Fig. 6(e) shows an enlarged view of the arterial region in the 548 nm image and perfusion image, in which the red line marks the artery in the 548 nm image, and the green line marks the artery in the perfusion image; Fig. 6(f) shows the cross-section of retinal arteries, the red curve is the cross-section of the artery in the 548 nm image, and the green curve is the cross-section of the artery in the perfusion image. From 6(c) and (f), it can be seen that the cross-section curve in the 548 nm image is wider than in the perfusion image, that is, the diameter measured in the 548 nm retinal image is ED and the diameter measured in the perfusion image is LD.

 

Fig. 6. Comparison of vascular cross-sections in 548 nm images and perfusion images. (a) 548 nm image, the yellow solid line box and the blue solid line box mark the representative venous and arterial regions, respectively; (b) Enlarged view of the venous region in the 548 nm image and perfusion image, in which the red line marks the vein in the 548 nm image, and the green line marks the vein in the perfusion image; (c) The cross-sections of retinal veins, the red curve is the cross-section of the vein in the 548 nm image, and the green curve is the cross-section of the vein in the perfusion image;(d) Perfusion image, the yellow dotted box and the blue dotted box mark the representative venous and arterial regions, respectively; (e) Enlarged view of the arterial region in the 548 nm image and perfusion image, in which the red line marks the artery in the 548 nm image, and the green line marks the artery in the perfusion image; (f) The cross-sections of retinal artery, the red curve is the cross-section of the artery in the 548 nm image, and the green curve is the cross-section of the artery in the perfusion image.

Download Full Size | PDF

2.4 Measuring principle of blood flow metrics

As a method for dynamic BFV evaluation, LSCI has the advantages of noninvasive, and high temporal resolution which has been used in retinal hemodynamic analysis [37,38]. The LSC calculation is usually performed in a sliding $H\times W\times F$ spatio-temporal window, where $H$ and $W$ represent the height and width of the SPI, and $F$ is the number of consecutive frames in time. BFV can be estimated by calculating the ratio of mean value to standard deviation in the spatio-temporal window [39,40], as Eq. (3) shown.

2.5 Measuring principle of oxygenation metrics

2.5.1 Retinal oxygen saturation and concentration analysis

The $\mathrm {SO_2}$ is measured by analyzing the absorption differences between oxygenated hemoglobin and deoxygenated hemoglobin at different wavelengths [43]. In spectral imaging, optical density (OD) is used to describe the capacity of spectral absorption, which can be quantified by Lambert Beer law, as Eq. (7) shown

On the basis of $\mathrm {SO_2}$, $\mathrm {C_{O2}}$ can be calculated according to Eq. (9)

2.5.2 Modelling of oxygen kinetics

The oxygen kinetics metrics include $\mathrm {DO_2}$, $\mathrm {MO_2}$, and OEF [48]. The $\mathrm {DO_2}$ is the product of blood flow and arterial oxygen content. The $\mathrm {MO_2}$ represents the oxygen metabolism which is the product of blood flow and arteriovenous oxygen difference. Finally, the OEF, defined as the ratio of $\mathrm {MO_2}$ to $\mathrm {DO_2}$, quantifies how much the retinal tissue extracts oxygen for metabolism. Combined with the concept proposed by Blair et al. [48], the expression is modified in ROKIA. The $\mathrm {DO_2}$ (unit: $mlO_2/dl$) is calculated as Eq. (10) shown

With high temporal resolution (80 fps), the pulsatility of these kinetic metrics in the cardiac cycle is available. The pulsatility is mainly caused by RBF, and subtle changes in $\mathrm {SO_2}$ and vascular diameters during cardiac cycles are ignored in this study. Figure 7 shows a representative pulsatility derived from a healthy subject to explain the process of metrics calculation. The pulsatility (Fig. 7) can be divided into two phases, systole (A to B) and diastole (B to C). By analyzing the pulsatility ($P$), systolic oxygenation rate ($(\sum _{A_x}^{B_x}P(x))/(B_x - A_x$)), diastolic oxygenation rate ($(\sum _{B_x}^{C_x}P(x))/(C_x - B_x$)), amplitude ($B_y - (A_y + C_y)/2$), and the mean value of oxygen kinetics metrics were analyzed.

 

Fig. 7. A representative pulsatility derived from a healthy subject. A to B represents the systole phase and B to C represents the diastole phase.

Download Full Size | PDF

2.6 ROKIA workflow

The workflow of ROKIA is shown in Fig. 8. Dual-wavelength images and SPI were collected, firstly. Based on these raw data, ED, LD, BFV, RBF, $\mathrm {SO_2}$, $\mathrm {C_{O2}}$, and pulsatility maps were calculated. Then, these maps were fused for vascular metrics, blood flow metrics, and oxygenation metrics analysis.

 

Fig. 8. Workflow of the ROKIA pipeline. Abbreviations: ED,external diameter ; LD, lumen diameter; $\mathrm {SO_2}$, retinal oxygen saturation; $\mathrm {C_{O2}}$, retinal oxygen concentration; BFV, retinal blood flow velocity; RBF, retinal blood flow; WT, wall thickness; WLR, wall-to-lumen ratio; A, artery; V, vein; $\mathrm {DO_2}$, oxygen delivery; $\mathrm {MO_2}$, oxygen metabolism; OEF, oxygen extraction fraction. A representative video is provided to demonstrate the dynamic velocity measurement results, as shown in Visualization 1.

Download Full Size | PDF

We adopt dual-wavelength images at 548 $nm$ and 605 $nm$ for $\mathrm {SO_2}$ calculation. Before that, ED was calculated from 548 $nm$ image. ODs could be calculated by the Lambert Beer law (Eq. (7)). The mean gray value of the intravascular pixel value was $I$, and the mean gray value of extravascular pixels was $I_0$. To reduce the effects of light scattering, for each pixel along central lines, two points that 15 pixels outside the central line were selected and the average intensity of these two points was used as the reflection signal before hemoglobin absorption, $I_0$. The connecting line of these two outside points was perpendicular to the central line. Then, the $\mathrm {SO_2}$ model was established according to Eq. (8). Finally, the $\mathrm {C_{O2}}$ could be calculated by combining $\mathrm {c_{HB}}$, $\gamma$, and $\mathrm {SO_2}$, as Eq. (9) shown.

It is worth noting that the precision of $\mathrm {SO_2}$ based on dual-wavelength imaging technology is susceptible to interference from uneven illumination, tissue absorption deviation, or tissue scattering, especially in small blood vessels. Therefore, researchers usually choose thicker arteries and veins(>50 $\mu$m) within optical disc-centered, 1.5 to 3 disc diameters areas for analysis [49]. This research followed the similar scheme. As shown in Fig. 8, the inner and outer diameter of the annular ROI was 1.5 and 3 times of disc diameter. Within this annular ROI, four main arterial vessel segments and four main venous vessel segments $\vec {V_{ROI}}$ (superior temporal, inferior temporal, superior nasal, and inferior nasal) were selected. All the metrics were calculated and averaged independently for each vessel segment.

For LSCI, 400 SPIs were recorded over 5 seconds, with a frame rate of 80 fps. After registration, the LSC algorithm was used for LSC images calculation. In this study, the spatial LSCI algorithm was used, and the size of spatio-temporal window was $3\times 3\times 1$ [50]. After LSC calculation, a total of 400 LSC images were obtained. For each LSC image, velocity maps could be calculated independently according to Eq. (4), (5). In order to improve the signal-to-noise ratio, we averaged these 400 velocity maps as the BFV map (As shown in Fig. 8). After that, LD was calculated from the BFV map, and then BRF map was further analyzed according to Eq. (6). In addition to the BFV and RBF maps, we also analyzed the pulsatilities. For vessels in the annular ROI (As shown in Fig. 8), each velocity map was calculated and averaged independently for each vessel segment. For each velocity map, we calculated and averaged velocity values for each vessel segment with the ROI map (Fig. 8) as a reference. After the above step, each vessel map will get 8 mean flow velocity values. We arranged these values according to the temporal axis to obtain 8 pulsatilities.

Based on the above procedures, vascular metrics, blood flow metrics, and oxygenation metrics were further analyzed. In vascular metrics analysis, eight mean values of eight segments in $\vec {V_{ROI}}$ were calculated independently. So did the LD map. Eight WT and WLR values were obtained according to Eq. (1), (2). Finally, for each metric, we averaged four arteries values and four veins values respectively. Similarly, blood flow metrics were calculated based on BFV, RBF, and ROI maps. The oxygenation metrics estimation were based on $\mathrm {SO_2}$, $\mathrm {C_{O2}}$, RBF, and ROI maps.

2.7 Subjects

For the reproducibility evaluation, eleven healthy subjects (9 males and 2 females) with an average age of 34.4$\pm$6.5 years (mean $\pm$ SD) were analyzed. Each subject performed 2 consecutive imaging on the same eye (oculus dextrus) with a 2-minute interval for rest. As a proof of concept study, we then analyzed 9 NPDR eyes (oculus dextrus) with an average age of 54$\pm$5.20 (6 males and 3 females) and nine age-matched healthy controls (58.22$\pm$5.09, 6 males and 3 females). The study was approved by the institutional review board of Peking University Shenzhen Hospital and informed consent was obtained from the subjects before the experiment. All procedures adhered to the tenets of the Declaration of Helsinki.

3. Results

3.1 Retinal oxygen kinetics imaging

Figure 9 shows representative imaging results from a 26-year-old, healthy male volunteer. Figure 9(a) and (b) show retinal images at 548 $nm$ and 605 $nm$. For veins, both images show higher absorption and lower gray value, and arteries show higher absorption at 548 $nm$ and lower absorption at 605 $nm$. Based on the above phenomena, the measurement results of $\mathrm {SO_2}$ and $\mathrm {C_{O2}}$ maps are shown in Fig. 9(c) and (d), where the color red indicates arteries with a higher percentage of oxygenated hemoglobin, and the color green indicates veins with a lower percentage of oxygenated hemoglobin. To validate the oxygen mapping qualitatively, we provided a standard artery-vein mapping labeled by an experienced ophthalmologist on color fundus image in Fig. 9(j). The ophthalmologist was blinded to the $\mathrm {SO_2}$. The red ’A’ represents arteries, and the black ’V’ represents veins. The qualitative comparison between Fig. 9(c) and (j) exhibits good consistency in distinguishing arteries and veins. The BFV and RBF maps are shown in Fig. 9(e) and (f). Due to the periodic pulsation of the heart, the pulsatility is shown in Fig. 9(i). The ED and LD maps are illustrated in Fig. 9(g) and (h).

 

Fig. 9. Results of Retinal Oxygen Kinetics Imaging in a healthy subject. (a) The image at 548 nm. (b) The image at 605 nm. (c) The $\mathrm {SO_2}$ map. (d) The $\mathrm {C_{O2}}$ map. (e) The BFV image. (f) The RBF image. (g) The ED map. (h) The LD map. (i) The pulsatilities in arteries and veins. (j) The standard artery-vein mapping labelledby an experienced ophthalmologist on color fundus image. The red ’A’ represents arteries, and the ’V’ represents veins. Abbreviations: ED, external diameter ; LD, lumen diameter; $\mathrm {SO_2}$, retinal oxygen saturation; $\mathrm {C_{O2}}$, retinal oxygen concentration; BFV, retinal blood flow velocity; RBF, retinal blood flow; WT, wall thickness; WLR, wall-to-lumen ratio; A, artery; V, vein; CFI, color fundus image.

Download Full Size | PDF

3.2 Retinal oxygen kinetics metrics

Two parameters indicating reproducibility, the coefficient of variation (COV) and the intraclass correlation coefficient (ICC), were analyzed with the oxygen kinetics metrics. The oxygen kinetic metrics and reproducibility results are summarized in Table 1. The highest reproducibility is observed in the $\mathrm {\overline {OEF}}$ (COV: 4.62%; ICC: 0.932), and the lower reproducibility is observed in the $\mathrm {WT^{a}}$ (COV: 12.87%; ICC: 0.754).

Table 1. Reproducibility Result of Retinal Oxygen Kinetic Metrics for Healthy Group^a

View Table | View all tables in this article

The vascular metrics in ROKIA show similar results with previous studies [51–53]. Both venous ED ($\mathrm {ED^{v}}$, $137.46\pm 16.71 \mu m$) and LD ($\mathrm {LD^{v}}$, $114.46\pm 12.74 \mu m$) are larger than arteries ($\mathrm {ED^{a}}$, $116.19\pm 16.71 \mu m$;$\mathrm {LD^{a}}$, $87.19\pm 15.00 \mu m$), while the WTs of arteries ($\mathrm {WT^{a}}$, $14.50\pm 4.14 \mu m$) and veins ($\mathrm {WT^{v}}$, $11.71\pm 6.01 \mu m$) are close. In addition, the ROKIA independently analyzes the arterial and venous metrics, which provides more information for future studies.

For $\mathrm {SO_{2}}$ analysis, the arteries $\mathrm {SO_{2}}$ ($91.38\pm 4.89 \%$), veins $\mathrm {SO_{2}}$ ($57.08\pm 8.89 \%$), and arteriovenous difference ($34.30\pm 5.08 \%$) are in accordance with previous researches [8,9]. For blood flow analysis, the average total arterial ($12.91\pm 3.82$ $\mu l/min$) and venous ($13.83\pm 2.37$ $\mu l/min$) blood flow are in close agreement with Murray’s law [14,54,55]. For oxygen kinetics analysis, the $\mathrm {\overline {DO_{2}}}$, $\mathrm {\overline {MO_{2}}}$, and $\mathrm {\overline {OEF}}$ are consistent with the previous researches [17,56,57]. All these metrics show well reproducibility, which means that ROKIA is a reliable technique that can be used for follow-up research.

3.3 Oxygen kinetic metrics changes in NPDR patients

Table 2 lists the results of a comparison study between NPDR and control groups. All the entire retinal oxygen kinetic metrics are compared between control and NPDR groups using unpaired t-tests. A P-value of <0.05 is considered statistically significant. No significant differences in age (P=0.15) and gender (P=1) are observed.

Table 2. Comparison of Retinal Oxygen Kinetic Metrics Between Control (N = 9) and NPDR (N = 9) Groups^a

View Table | View all tables in this article

Compared with control group, NPDR group has higher arterial WT ($\mathrm {WT^{a}}$), higher WLR ($\mathrm {WLR^{a}}$) and lower LD ($\mathrm {LD^{a}}$). But, there are no significant differences in arterial ED ($\mathrm {ED^{a}}$) and veinous metrics (including $\mathrm {ED^{v}}$, $\mathrm {LD^{v}}$, $\mathrm {WT^{v}}$, and $\mathrm {WLR^{v}}$). As shown in Table 2, arterial $\mathrm {SO_2}$ ($\mathrm {SO_2^{a}}$), $\mathrm {C_{O2}}$ ($\mathrm {C_{O2}^{a}}$) and venous $\mathrm {SO_2}$ ($\mathrm {SO_2^{v}}$), $\mathrm {C_{O2}}$ ($\mathrm {C_{O2}^{v}}$) are higher in the NPDR group than the control group. However, the arteriovenous difference of $\mathrm {SO_2}$ ($\mathrm {SO_2^{av}}$) and $\mathrm {C_{O2}}$ ($\mathrm {C_{O2}^{av}}$) are lower in NPDR group. The $\mathrm {BFV^{a}}$ and $\mathrm {BFV^{v}}$ show no significant differences but the $\mathrm {RBF^{a}}$ and $\mathrm {RBF^{v}}$ show a significant decrease in the NPDR group than the control group. Compared with the control group, the NPDR group has lower value in $\mathrm {DO_2^{sr}}$, $\mathrm {DO_2^{dr}}$, $\mathrm {\overline {DO_{2}}}$, $\mathrm {MO_2^{sr}}$, $\mathrm {MO_2^{dr}}$, $\mathrm {\overline {MO_2}}$, and $\mathrm {\overline {OEF}}$. Nevertheless, the $\mathrm {DO_2^{amp}}$ and $\mathrm {MO_2^{amp}}$ are higher in the NPDR group as compared to the control group.

4. Discussions

The feasibility of ROKIA is demonstrated in this study. Retinal structural and functional imaging in a single system can fully integrate structural and functional information and further analysis of $\mathrm {DO_{2}}$, $\mathrm {MO_{2}}$, and OEF. In addition, the measurement of LD, WT and WLR can also be preformed with the proposed single system. Based on the ROKIA classic metrics including $\mathrm {ED^{a}}$, $\mathrm {LD^{a}}$, $\mathrm {WT^{a}}$, $\mathrm {WLR^{a}}$, $\mathrm {ED^{v}}$, $\mathrm {LD^{v}}$, $\mathrm {WT^{v}}$, $\mathrm {WLR^{v}}$, $\mathrm {SO_2^{a}}$, $\mathrm {SO_2^{v}}$, $\mathrm {SO_2^{av}}$, $\mathrm {C_{O2}^{a}}$, $\mathrm {C_{O2}^{v}}$, $\mathrm {C_{O2}^{av}}$, $\mathrm {BFV^{a}}$, $\mathrm {BFV^{v}}$, $\mathrm {RBF^{a}}$, $\mathrm {RBF^{v}}$, $\mathrm {\overline {DO_{2}}}$, $\mathrm {\overline {MO_{2}}}$, $\mathrm {\overline {OEF}}$ are measured. Furthermore, additional six oxygen kinetics metrics are proposed based on the advantages of ROKIA technology with high temporal resolution. These 27 metrics systematically describe the functional and structural information. The reproducibility of these metrics is good, since the ICC of oxygen kinetics metrics are greater than 0.75, and the COVs are less than 16%.

As a representative oxygen metabolism disease, NPDR has complex evolution process of oxygen metabolism dysfunction. Unfortunately, the lack of monitoring technology limits the study of oxygen metabolism in NPDR patients. Previous studies have shown that patients with NPDR have higher $\mathrm {SO_2^{a}}$ and $\mathrm {SO_2^{v}}$, but lower $\mathrm {SO_2^{av}}$ [4,58,59]. The results in Table 2 demonstrate similar conclusions. The values of $\mathrm {SO_2^{a}}$ ($91.73\pm 2.15$%) and $\mathrm {SO_2^{v}}$ ($55.69\pm 3.91$%) in healthy controls are smaller than values of $\mathrm {SO_2^{a}}$ ($97.14\pm 6.09$%) and $\mathrm {SO_2^{v}}$ ($64.98\pm 7.07$%) in NPDR patients, and the value of $\mathrm {SO_2^{av}}$ in healthy controls ($36.10\pm 3.94$%) is greater than that ($32.16\pm 2.54$%) in NPDR patients. These results prove that ROKIA has reliable repeatability and sensitivity for classic $\mathrm {SO_2}$ analysis.

Unlike $\mathrm {SO_2}$, the results of RBF in diabetic retinopathy (DR) patients are controversial in previous studies. Some studies have found that the RBF of NPDR patients is significantly lower than that of healthy controls [59,60], which is similar to our result. However, some other studies have reported the increase of RBF in NPDR patients [61,62]. There are several explanations for the ambiguity. First, some long-term follow-up studies find that RBF of diabetic patients without DR increases gradually. However, when DR occurs, RBF decreases with the aggravation of DR [61,63,64]. Second, RBF is affected by several factors, such as blood glucose concentration, blood pressure, and heart rate among others [65]. Third, only a few NPDR patients have been investigated for demonstrating proof of the concept of ROKIA technology in this study. More patients will be studied in the future. In addition, there is no significant difference between $\mathrm {BFV^{a}}$ and $\mathrm {BFV^{v}}$, but significant difference exists between $\mathrm {RBF^{a}}$ and $\mathrm {RBF^{v}}$. This phenomenon indicates that vascular metrics need to be considered during the RBF calculation.

Besides above metrics, oxygen kinetics ($\mathrm {\overline {DO_{2}}}$, $\mathrm {\overline {MO_{2}}}$, and $\mathrm {\overline {OEF}}$) in healthy controls are also comparable with previous researches [56,57]. The results of $\mathrm {\overline {DO_{2}}}$ and $\mathrm {\overline {MO_{2}}}$ of NPDR patients in this study are lower than those in healthy controls. Recently, based on Doppler OCT and dual-wavelength oximetry, Rahimi et al. performed oxygen kinetics analysis and found a decrease in $\mathrm {\overline {DO_{2}}}$ in the NPDR group, which is consistent with the results of this study. The $\mathrm {\overline {OEF}}$ in healthy controls ($0.36\pm 0.04$) are significantly higher (P=0.033) than that in NPDR patients ($0.31\pm 0.04$). The lower $\mathrm {\overline {OEF}}$ in NPDR patients is identical with the decreased oxygen extract capacity in NPDR patients described in previous studies [14,66]. In addition, the $\mathrm {DO_{2}}$ and $\mathrm {MO_{2}}$ in systolic and diastolic phases are analyzed independently with the advantage of the high temporal resolution of ROKIA. These oxygen kinetics metrics have sufficient sensitivity for detecting the abnormal oxygen metabolism of NPDR in the early stage, which will be helpful for DR diagnosis and intervention.

Compared to structural imaging systems, ROKIA technology has stronger capability in functional imaging. Compared with dual-wavelength imaging device or LSCI device, ROKIA technology deeply integrates these two technologies and calculate unique parameters of this technology: ED, LD, WT, WLR, $\mathrm {DO_2}$, $\mathrm {MO_2}$ and OEF. However, there are some limitations in the present study. Firstly, the ROKIA technology achieves oxygen kinetics imaging with large FOV, high temporal resolution and high spatial resolution. But there are only 4 arteries and 4 veins in the area of 1.5-3x optic disc diameter are considered in the calculation of oxygen kinetics metrics in this study. Hence, the oxygen kinetics metrics in microcirculation needs to be further developed. Secondly, for $\mathrm {SO_2}$ analysis, the bandwidth of the light source in this work seems a bit large, which may reduce the precision. In future research, we will further reduce the bandwidth of LED light sources. Thirdly, for LD measurement, velocity near vessel walls will decrease, which may cause inaccurate boundary positioning. Finally, considering that the LSC is affected by many factors (such as velocity, wavelength, system design, etc.), the $\beta$ value needs to be set carefully for each device.

5. Conclusion

In this article, we propose a novel retinal oxygen kinetics imaging (ROKIA) technology. By integrating dual-wavelength imaging technology with LSCI technology, ROKIA achieves high temporal resolution, and structural-functional imaging. This is the first demonstration of the vascular metrics (including ED, LD, WT, and WLR), blood flow metrics (including BFV and RBF), oxygenation metrics (including $\mathrm {SO_2}$, $\mathrm {C_{O2}}$, $\mathrm {DO_2}$, $\mathrm {MO_2}$, and OEF) via a single system, to the best of our knowledge. We verified the 27 metrics proposed by ROKIA in healthy and NPDR groups. The results show that the oxygen kinetics metrics in NPDR and healthy groups are significantly different. In addition, significant differences among OEF are first observed, which may imply abnormal oxygen metabolism in the early stages of DR. ROKIA simplifies clinical workflow and comprehensively analyzes the oxygen kinetics which is essential for diagnosis, screening, and prognosis. At the same time, the scalability of ROKIA provides a fresh idea for future research.