Inference of Capillary Nonperfusion Progression on Widefield OCT Angiography in Diabetic Retinopathy

Abstract

Purpose: The purpose of this study was to explore the spatial patterns of the nonperfusion areas (NPAs) on widefield optical coherence tomography angiography (OCTA) images in diabetic retinopathy (DR) and to investigate their associations with NPA progression and DR severity.

Methods: We prospectively enrolled 201 eyes from 158 patients with DR. Widefield images were obtained using a swept-source OCTA device (Xephilio OCT-S1), followed by the creation of 20-mm (1614 pixels) en face images. Nonperfusion squares (NPSs) were defined as 10 × 10-pixel squares without retinal vessels. Eyes with high-dimensional spatial data were mapped onto a two-dimensional space using the uniform manifold approximation and projection algorithm and divided by clustering. The patterns of NPA distribution were statistically compared between clusters.

Results: All eyes were mapped onto a two-dimensional space and divided into six clusters based on the similarity of NPA distribution. Eyes in clusters 1 and 2 had minimal and small NPAs, respectively. Eyes in clusters 3 and 4 exhibited NPAs in the temporal and inferotemporal regions, respectively. Eyes in cluster 5 displayed NPAs in both superonasal and inferonasal areas. The unique NPA distributions in each cluster encouraged us to propose eight possible pathways of NPA progression. DR severity was not equal between clusters (P < 0.001), for example, 8 (15.7%) of 51 eyes and 15 (65.2%) of 23 eyes had PDR in clusters 1 and 5, respectively.

Conclusions: Dimensionality reduction and subsequent clustering based on the NPA distribution on widefield OCTA enabled the inference of possible NPA progression in DR.

Diabetic retinopathy (DR) is a major cause of severe vision loss in working age populations worldwide.¹ Hyperglycemia initiates and promotes active microvascular lesions, for example, vascular hyperpermeability and angiogenesis, and the progressive degeneration of neurovascular unit, for example, nonperfusion areas (NPAs) and retinal neurodegeneration.² In the era of anti-VEGF therapy, diabetic macular edema (DME) and proliferative diabetic retinopathy (PDR) are clinically treatable in many cases, whereas NPAs and diabetic macular ischemia (DMI) exacerbate irreversible visual impairment in patients with diabetes.^3–10 This has prompted medical researchers to hypothesize “diabetic retinal disease” with respect to the worsening of vascular and neural function.¹¹ 

Quantitative and qualitative profiling of diabetic capillary nonperfusion is in progress. Automatically quantified NPAs on (ultra)wide field fluorescein angiography (FA) or optical coherence tomography angiography (OCTA) images are associated with DR severity.^4,12,13 Qualitative analyses documented unique distribution of NPAs, which may be explained by the three-dimensional characteristics of retinal vessels in each location.^14–16 Diabetic NPAs develop throughout the retinas, but predominantly in the midperiphery.¹⁷ Major arterioles segment all three capillary layers in the extramacular regions, leading to vulnerability to NPA progression. In contrast, overlapping perfusion by multiple arterioles and a seamless network of deep capillaries might maintain the perfusion status in the macula.^18,19 Radial peripapillary capillaries (RPCs) and imbalanced perfusion pressure between the nasal and temporal quadrants might contribute to the capillary nonperfusion in the temporal subfield.²⁰ However, the patterns and processes of NPA progression remain to be systematically understood in DR. 

Recent advances in data science have prompted biologists and clinicians to analyze the trajectory inference of cellular differentiation in development and cancer progression in individuals, respectively.²¹ High-dimensional datasets of mRNA levels are mapped onto two-dimensional spaces using dimensionality reduction algorithms, for example, principal component analysis and uniform manifold approximation and projection (UMAP), based on the similarity of gene expression patterns.²⁰ Subsequent clustering and pseudotime analyses in clinical datasets reveal the patterns and progression profiles in several diseases.^22,23 Despite its usability, such methodological approaches have never been applied to the analyses of DR progression. 

In this study, we aim to analyze the patterns of NPA distribution using widefield OCTA images via UMAP and clustering, infer the pathways of its progression statistically, and investigate its association with the international DR severity grades. 

Methods

Participants

In this prospective, observational, cross-sectional case series, we enrolled patients with DR who were examined at the Department of Ophthalmology in Kyoto University Hospital. This study was conducted under the approval of the Kyoto University Graduate School and Faculty of Medicine Ethics Committee and in adherence to the tenets of the Declaration of Helsinki. Written informed consent was obtained from each participant before the study. 

Individuals with DR for whom widefield OCTA images of sufficient quality were obtained were included in the study. The exclusion criteria were the presence of media opacities interfering with visual acuity (VA) measurement or image acquisition, other chorioretinal diseases, a history of vitrectomy, cataract surgery within 3 months, prior ocular steroid treatment, and prior anti-VEGF treatment; and an axial length < 22 mm or > 26 mm. We further excluded eyes with severe image artifacts on en face OCTA images or weak signal strength (< 5). 

Image Acquisition and Processing

After comprehensive ophthalmic examinations, the best-corrected decimal VA was measured and converted into logarithm of the minimum angle of resolution (logMAR). The axial length was measured using partial coherence interferometry (IOL Master; Carl Zeiss Meditec, Inc.). Ultra-widefield color fundus photographs were taken using Optos 200Tx (Optos PLC), on which fundus findings, for example, neovascularization, were assessed, and graded for DR severity according to the International Clinical Disease Severity Scale for Diabetic Retinopathy.²⁴ Two retinal specialists evaluated the fundus photographs, and any disagreements were discussed until they reached an agreement. 

Two en face, widefield OCTA images centered on the upper and lower quadrants were acquired with a scanning area of nominal 20 (height) × 23 (width) mm (1614 × 1856 pixels) using a swept-source OCTA device (Xephilio OCT-S1 device; Canon), as described previously.²⁵ After the creation of a montage image of superficial and deep layers (from the internal limiting membrane to the outer plexiform layer) using the built-in manufacturer's software, the 1614-pixel (nominal 20-mm) diameter circle centered on the fovea in the images were evaluated, as described previously.²⁵ We did not select the lamellar NPAs in the superficial or deep layers, because the automatic segmentation was not accurate in the extramacular areas in the default setting of this software. The images of the left eye were inverted horizontally to position the nasal subfield on the right-hand side. The circle with a diameter of 250 pixels located 375 pixels to the nasal side and 25 pixels to the superior side of the fovea was not evaluated in order to exclude various OCTA signal levels of the optic disc in each location in each eye. 

Edges of the retinal vessels were automatically detected using the Canny Edge Detector plugin of ImageJ software (NIH, http://imagej.nih.gov/ij/), as described previously.²⁶ The image was divided into squares of 10 × 10 pixels. The pixels of vessel edges were automatically counted in each 10 × 10-pixel square using ImageJ. The squares with no pixels of vessel edges were defined as nonperfusion squares (NPSs). The squares were assigned coordinates with the fovea as the origin and binarized based on the NPS. In other words, each eye had 19,623 dimensional data of the NPAs. 

We defined the NPS ratio for each square as follows:  

\begin{eqnarray*} {\text{NPS ratio}} = \frac{{{\text{The number of eyes with NPSs}}}}{{{\text{The total number of eyes examined}}}}. \end{eqnarray*}

 

Dimensionality Reduction and Clustering

UMAP, a technique for nonlinear dimensionality reduction, was used to visualize eyes with high-dimensional data of NPSs in a two-dimensional space, based on the characteristics of the NPSs distribution.²⁷ The Euclidean distance between the 3 nearest neighbors for each point with a minimum distance of 0.1 was used to develop the manifold. All other parameters remained at the default values reported by McInnes et al.²⁷ Each eye was mapped onto a two-dimensional space according to the similarity of the NPSs distribution. Subsequently, we used the elbow method and silhouette analysis to determine the optimal number of clusters and performed clustering using the k-means algorithm (Supplementary Fig. S1). 

To compare the characteristics of NPSs distribution between individual clusters, the 20-mm circle was divided into 17 sectors (Supplementary Fig. S2): a central 1-mm area, 8 sectors (0–1.5, 1.5–3, 3–4.5, 4.5–6, 6–7.5, 7.5–9, 9–10.5, and 10.5–12 o'clock) of the inner ring (from a 0.5-mm radius to a 5-mm radius), and 8 sectors of the outer ring (from a 5-mm radius to a 10-mm radius). The percentages of NPSs to the total squares in each sector were calculated as follows:  

\begin{eqnarray*} && {\text{NPS percentage }}\left( {\rm{\% }} \right)\\ && = \frac{{{\text{The number of NPSs in a sector}}}}{{{\text{The number of all squares in a sector}}}}{\rm{ \ \times \ 100}}. \end{eqnarray*}

 

Inference of NPAs Progression

To infer the progression of NPAs, we performed two analyses: (1) statistical comparisons of NPSs between clusters, and (2) probabilistic calculation of possible transition from one cluster to another. 

Empirically, the NPAs progress gradually but not stepwise in DR, and eyes were mapped according to the similarity of NPA distribution. We therefore considered that the transition is possible between nearby clusters. Because typical NPAs progress irreversibly, we hypothesized that NPS ratios would sustain or increase in any retinal areas. First, we compared the NPSs ratios in each sector statistically, and considered that NPS ratios would not be different or would statistically increase as NPAs progress. We then repeated such analyses in all 17 sectors. If cluster X had higher NPS ratios in one or more sectors than cluster Y, we defined the possible NPA progression from cluster Y to cluster X. 

In the other analysis, based on a hypothesis that transitions of clusters would occur between neighboring clusters, we counted the number of all combinations of eyes between two clusters. A possible transition was defined as the combinations in which there was the same or greater numbers of NPSs across all sectors. The percentages of possible transition were calculated as follows:  

\begin{eqnarray*} && {\text{Percentage of possible transition }}\left( {\rm{\% }} \right) \\ && = \frac{{{\text{The number of possible transitions}}}}{{{\text{The total number of combinations}}}}{\rm{ \ \times \ 100}}. \end{eqnarray*}

 

Statistical Analyses

All values were expressed as the median and interquartile range (IQR). Statistical significance was set at P < 0.05. The Kruskal-Wallis test with a Bonferroni correction was performed to compare continuous variables. The Chi-square test was used to compare categorical variables. All statistical analyses were conducted using SPSS (version 24; IBM). 

Results

Multiple Patterns of NPA Distribution

We enrolled 201 eyes from 158 participants that met the inclusion and exclusion criteria in this study (Table 1). The NPS ratio of all cases was shown on the pseudocolored rate map image (Fig. 1A). The NPS ratios were lower along the vascular arcade and around the optic disc, whereas the ratios were higher in the extramacular areas and the temporal subfield. The UMAP algorithm mapped 201 eyes onto a 2-dimensional space based on the characteristics of the NPS distribution (Fig. 1B). We first investigated the NPS counts in each eye on the UMAP projection and found that eyes with minimal NPAs were found in the upper areas of the UMAP and corresponded to those in the initial stage of NPAs severity. Eyes were distributed toward the left, the right lower, or the lower side. Most eyes with a greater number of NPSs were found at the marginal zones. 

Table 1.

 

Patient Characteristics

 

Figure 1.

 

Eyes with DR on the two-dimensional space using UMAP based on the NPAs distribution. (A) The NPS ratio in each square on the pseudocolored rate map in all 201 eyes with DR. (B) Each eye is shown as a dot on the two-dimensional pseudocolored map using UMAP algorithm. (C) K-means clustering divides 201 eyes into 6 groups. (D) Boxplots show the differences of NPS percentages between clusters in all 201 eyes. *P < 0.05; †P < 0.01; ‡P < 0.001.

Subsequent clustering using the k-means algorithm divided 201 eyes into 6 groups (Figs. 1C, 1D, Supplementary Fig. S3). The upper, left, and right lower areas corresponded to clusters 1, 3, and 4, respectively. Eyes within the lower areas were further divided into three clusters: those in the middle left (cluster 2), lower left (cluster 6), and lower right (cluster 5) areas. The NPS ratios in eyes of each cluster were exhibited in pseudocolored rate maps; minimal NPAs in cluster 1 (Fig. 2A); a random distribution of mild NPAs in cluster 2 (Fig. 2B); greater NPAs in the temporal regions in cluster 3 (Fig. 2C); NPAs localized to the inferotemporal areas in cluster 4 (Fig. 2D); severe NPAs in both the superonasal and inferonasal areas in cluster 5 (Fig. 2E); and NPAs limited in the macula in cluster 6 (Fig. 2F). 

Figure 2.

 

NPS ratios in each square in eyes of each cluster. (A) The heatmap of cluster 1 showing minimal NPSs. (B) In eyes of cluster 2, lower NPS ratios are found throughout the retinas except the vascular arcade and optic disc. Higher NPS ratios in the temporal subfield in eyes of cluster 3 (C), in the inferotemporal subfield in eyes of cluster 4 (D), in both superonasal and inferonasal subfields in eyes of cluster 5 (E), and in the macula in eyes of cluster 6 (F).

Inference of the Progression of Capillary Nonperfusion

We used statistical and probabilistic methods to infer the progression of NPAs. After dividing the 20-mm circle OCTA image into 17 sectors, we compared the NPS percentages in each sector between clusters (Fig. 3, Supplementary Fig. S4). In the central sector, the NPS ratios were the lowest in eyes of cluster 1 and the highest in eyes of cluster 6, respectively (see Fig. 3A). The superior and inferior sectors of the outer ring had higher ratios in eyes of cluster 5 (see Fig. 3B, Supplementary Fig. S4). The NPS ratios in the temporal sectors were higher in eyes of cluster 3 than those of other clusters (see Fig. 3C). 

Figure 3.

 

Comparisons of the NPS percentages between clusters in each sector. Boxplots of each cluster in sector 1 (A), sector 11 (B), sector 14 (C), sector 15 (D), and sector 17 (E). *P < 0.05; †P < 0.01; ‡P < 0.001. (F) Possible NPA progression between clusters. Black dot = centroid of each cluster; and arrow = possible transition between clusters.

Considering the statistical differences in NPS ratios across all sectors between clusters, we identified 8 possible patterns of NPA progression: from clusters 1 to 2, 3, and 4; from clusters 2 to 3, 4, 5, and 6; and from clusters 4 to 5 (Fig. 3F). Cluster 2 was connected to the other 5 clusters, making it likely to be a central node. The probabilistic analyses demonstrated that the percentages of possible transition from clusters 1 to 2, 3, and 4 were higher, whereas those from clusters 2 to 3, 4, 5, and 6 were lower (Table 2). 

Table 2.

 

The Percentage of Possible Transitions in all 17 Sectors Between 2 Clusters

 

Association Between the Patterns of NPA Distribution and DR Severity

When examining the relationship between international DR severity grades and the presence of NPA, we found that eyes with mild nonproliferative diabetic retinopathy (NPDR) were mostly concentrated in the superior part of the UMAP projection (Fig. 4A). Conversely, eyes with moderate NPDR were scattered throughout. Notably, those with severe NPDR and PDR were situated near each other at the marginal zones. Eyes with retinal neovascularization (NVE), neovascularization of the disc (NVD), more than 20 intraretinal hemorrhages in each of the 4 quadrants, and definite venous beading (VB) in 2 or more quadrants were also located at the marginal zones (Figs. 4B–E). Eyes with prominent intraretinal microvascular abnormalities (IRMAs) in one or more quadrants were distributed randomly throughout the areas (Fig. 4F). 

Figure 4.

 

DR severity grades of each eye in the UMAP projection. Each dot corresponds to each eye with DR severity grade as indicated (A), in the presence of NVE (B), NVD (C), more than 20 intraretinal hemorrhages in each of 4 quadrants (hem; D), definite VB in 2 or more quadrants (E), and IRMAs (F).

Significant differences were observed between clusters for the following variables: sex, logMAR, DR severity, multiple intraretinal hemorrhages, IRMAs, NVE, NVD, and prior panretinal photocoagulation (PRP; P < 0.001; Table 3). 

Table 3.

 

Comparisons of Each Parameter Between Clusters

 

Discussion

In the current study, we used dimensionality reduction using UMAP and subsequent clustering for automatically detected NPAs on widefield OCTA images in DR. This approach elucidated several patterns of NPAs distribution in DR, namely; minimal type, mild type, temporal type, inferotemporal type, both superior and inferior types, and macular type. Comparisons of the NPA extents among clusters encouraged us to infer the unique progression of NPAs in this cross-sectional study, which should be confirmed by future longitudinal studies. Further investigation demonstrated that the patterns of NPA distribution are associated with DR severity to some extent, although the amounts of NPAs did not always correspond to DR severity. 

The NPS ratios were the lowest around the optic disc and along the superotemporal and inferotemporal vascular trunks in all 201 eyes. This observation may be explained by an additional capillary layer, the radial peripapillary capillaries, and higher perfusion pressure there.^16,28 Furthermore, NPAs were more frequently seen in the extramacular areas than in the macula. This may depend on the predisposition of vascular structure specific to the extramacular and macular regions.¹⁸ Because the temporal raphe is the most distal to the optic disc, perfusion pressure is the lowest and capillary nonperfusion might be more likely to progress in the temporal subfield.^20,29 

Eyes in cluster 1 had minimal NPAs, although some eyes suffered from PDR. We were able to evaluate vascular lesions from the macula to the midperipheral regions using widefield OCTA images, although peripheral images could not be obtained. Recent publications showed that some cases have peripheral lesions including NPAs on ultrawide field FA images.^4,12,13 Future studies should elucidate the status of capillary nonperfusion in the periphery on ultrawide field montage of OCTA images. Empirically, some eyes of young patients with PDR had no NPAs, which might not be compatible with the general concept that VEGF derived from NPAs develops neovascularization. 

Clusters 3 (temporal type) and 5 (both superior and inferior types) showed an exclusive relationship with respect to the NPA distribution. We therefore hypothesized that blood flow redirects toward the remaining patent vessels after the NPA extension in the specific areas. When NPAs extended in the temporal subfield, blood flow was maintained in the superior and inferior peripheral regions in eyes of cluster 3. We might explain that the reduced amounts of capillaries in the temporal subfield increase perfusion pressure and concomitant resistance to capillary obstruction in other retinal regions. In eyes of cluster 5, segmental perfusion in the superonasal and inferonasal subfield might contribute to NPA progression there. As a result, the perfusion pressure might increase in the temporal subfield, and the seamless network in the deep capillaries there might also make NPAs progression less likely.¹⁸ 

In eyes of cluster 4, the NPA in the inferotemporal region could be attributed to 2 factors: the trait during the development and the perfusion pressure. First, the retina is thinner and vascular density is lower in the inferior subfield than in the superior subfield in healthy adults.^30,31 This might lead to reduced vascular density or segmental perfusion of the extramacular regions in the inferior subfield. The second point is that lower perfusion pressure on the temporal side than on the nasal side makes it more likely for NPA to develop in the inferotemporal region. The comparative study may support the possible transition from clusters 4 to 5, although they were far from each other in the UMAP projection. 

The heatmaps of clusters 1 and 6 seemed similar to each other, although there were statistical differences in NPS percentages of several sectors. The main difference was the NPS counts in the central 1 mm, and we speculated that eyes in cluster 6 might correspond to those with DMI.³ The foveal avascular zone is the most distal to the optic disc and parafoveal capillaries are composed of only one or two layers, compared to three or more layers in other areas within the macula.^14,15,32 Additionally, deep capillaries may reduce in eyes with DME, and VEGF derived from retinal neurons might also decrease in eyes with diabetic neurodegeneration.^33–35 However, it should be noted that we cannot exclude the possibility of artifacts from the OCTA machine with an artificial intelligence (AI)-driven denoise algorithm or the false negative in the image processing of vessels with high density.^36,37 

The diversity in the NPA distribution suggests that multiple factors influence the progression of capillary nonperfusion. In this cross-sectional study, we could not guarantee the processes of NPA progression completely, although the statistical comparisons suggested that there are eight pathways of the NPA progression. The transition from cluster 1 to 2 might correspond to the initiation; at random development of minimal NPAs. Because clusters 3 and 4 were near cluster 1 on the UMAP projection, the unknown and specific factors or predispositions might promote the NPAs extension to the specific regions vigorously. After the initiation, four pathways for NPA progression might lead to the specific distribution as discussed above. However, the probabilistic analyses suggest less frequency of NPA progression from cluster 2 to clusters 3, 4, 5, and 6. The probabilities of the transitions were higher than those of random transitions (data not shown). However, future longitudinal studies should confirm the pathways for NPA progression and the actual probability for the transition between clusters per year. 

There are several limitations to this study. The inclusion and exclusion criteria were applied in this single-center study, which may result in selection bias. The OCTA machine used in this study delineates angiographic images from the macula to the midperiphery, but not to the periphery. The image quality may affect the automatic quantification. Because the longer time was required for image acquisition of wide-field OCTA, images of sufficient quality were not obtained in eyes with poor fixation and concomitantly there was another selection bias. The artifacts at OCTA image acquisition may modify the subsequent assessment.^38,39 We used images of superficial and deep layers in this study, and future analyses of layer-specific NPAs would provide more detailed and accurate data. Although image processing in this study had advantages in the NPA distribution, the vessel edges could not be detected completely, and the NPA amounts were approximately measured. The presence of NVE could potentially lead to an underestimation of the NPA. We selected a UMAP algorithm to perform dimensionality reduction, and future studies may compare the results from other algorithms. We tried to infer the NPA progression from the macula to the midperiphery using the cross-sectional study, although future longitudinal studies should be planned to confirm the characteristics of the NPA progression in the whole retinas and to elucidate whether the patterns of NPA progression explain the differences in DR severity or visual impairment.⁴⁰ 

In conclusion, the dimensionality reduction using the UMAP algorithm and subsequent clustering divided DR eyes into six clusters based on the NPA distribution in the widefield OCTA images. Statistical and probabilistic analyses allowed us to infer the pathways for the NPA progression. 

Acknowledgments

Funded by a Grant-in-Aid for Scientific Research of the Japan Society for the Promotion of Science (Grant Number: 23K09004). The funding organization had no role in the design or conduct of this research. 

Disclosure: M. Yoshida, None; T. Murakami, None; K. Kawai, None; K. Nishikawa, None; K. Ishihara, None; Y. Mori, None; A. Tsujikawa, None