Three-dimensional structural and angiographic evaluation of foveal ischemia in diabetic retinopathy: method and validation

: Optical coherence tomography angiography (OCTA) allows us to noninvasively investigate foveal ischemia, a key feature of diabetic retinopathy (DR). However, the sizes of the foveal avascular zone (FAZ) have a significant variation in a normal population, preventing the objective assessment of pathological enlargement of FAZ due to capillary dropout. Based on the relationship between FAZ and ganglion cell complex (GCC) thickness in normal eyes, we defined a theoretical baseline FAZ (tbFAZ) on structural OCT and measured 2D and 3D vessel density in its vicinity on the simultaneously acquired OCTA in normal and diabetic eyes. We found that the structure-based tbFAZ was a reliable reference to identify foveal ischemia and that the 3D vessel density demonstrated ischemia more effectively than the 2D method. The proposed 3D para-FAZ vessel density correlates well with DR severity and potentially is a useful diagnostic biomarker, especially in the early stages of DR.

In this study, we propose a method of 3D foveal ischemia quantification that addresses each of these challenges and evaluate its utility in clinical setting. We defined the theoretical baseline FAZ (tbFAZ) using the relationship between the ganglion cell complex (GCC) thickness and the FAZ in healthy eyes as the basis for where the FAZ should be with the assumption that the laminar retinal structure remains constant, especially in the early stages of DR. We used this tbFAZ to define a para-FAZ region to evaluate the vessel density. Then, we calculated the 3D vessel density of this region, and compared it to the 2D method and the FAZ-based metrics in the evaluation of clinical DR.

Subjects
Healthy and diabetic participants were recruited from the Casey Eye Institute, Oregon Health & Science University (OHSU). The study's protocol was approved by the Institutional Review Board of OHSU and complied with the Declaration of Helsinki and the Healthy Insurance Portability and Accountability Act of 1996. The participants underwent full ophthalmic examinations, and were assigned into three predetermined groups based on photographic gradings based on the Early Treatment Diabetic Retinopathy Study (ETDRS) severity levels [25][26][27]. The groups were diabetes mellitus (DM) without retinopathy; mild to moderate nonproliferative DR (NPDR); and severe NPDR or proliferative DR (PDR). Eyes of healthy participants served as control. Exclusion criteria included presence of glaucoma, prior intraocular surgery, and significant ocular pathologies such as uveitis and retinal vascular occlusions, etc. We excluded scans with signal strength index (SSI) less than 55. Axial length was obtained for the correction of retinal magnification in all participants with the IOL Master (Carl Zeiss AG, Oberkochen, Germany).
A total of 102 subjects-22 healthy controls; 24 with DM without retinopathy; 25 with mild to moderate NPDR; and 31 with severe NPDR or PDR-were included in this study. Diabetic macular edema (DME) [28][29][30], as defined by the presence of cystic intraretinal fluid on the macular OCT scans, was present in 11 eyes with mild to moderate NPDR and 16 eyes with significant NPDR and PDR ( Table 1). The severity of DME were categorized by whether the cystic fluid occurred in GCC (significant DME), or restricted only at outer retina (mild DME).

Image acquisition
Volumetric macular scans (304 × 2 × 304 A-lines covering a 3 × 3 mm 2 area) centered at the fovea were obtained with a commercial spectral-domain OCT system (Avanti RTVue-XR, Optovue, Inc.). This system operates at an axial scan rate of 70 kHz and its spectrum is centered at 840 nm with a full-width half maximum bandwidth of 45 nm. It has an optical resolution in retinal tissue of 5 µm and 15 µm on axial and transverse directions, respectively, and a digital sampling interval 15 × 15 × 3 µm 3 /voxel. The split-spectrum amplitudedecorrelation angiography (SSADA) algorithm [31] generated the OCTA. PR-OCTA [32,33] and irb-BMS algorithms [34][35][36] further processed the images on MATLAB 2018b (MathWorks, Natick, MA, USA) platform.
Segmentations of the vitreous and inner limiting membrane (ILM) interface, inner plexiform layer (IPL) and inner nuclear layer (INL) interface, outer plexiform layer (OPL) and outer nuclear layer (ONL) interface as well as ellipsoid zone (EZ) were automatically performed by a guided bidirectional graph search algorithm [37]. As reported previously [37], the accuracy of the automatic segmentation algorithm 3.04 ± 0.48 µm for the vitreous/ILM interface, 6.20 ± 7.45 µm for the IPL/INL interface and 3.32 ± 1.76 µm for the EZ. In our study population, 5% of the B-scans were manually corrected for the minor segmentation errors by an expert grader. It took 3s to review each B-scan and 20s to manual correct if there were any segmentation errors. The inner retina is the region between ILM and OPL, and the ganglion cell complex (GCC) (Fig. 1(B)) is defined as the combination of nerve fiber layer (NFL), ganglion cell layer (GCL), and IPL [38]. The maximum flow projection was used to generate inner retinal angiogram ( Fig. 1(A)). Two independent graders manually delineated the angiographic FAZ boundaries in healthy controls. Inter-grader agreement was studied. The obtained angiographic FAZ area was highly correlated between the two graders, with the mean ± SD difference of −0.008 ± 0.011 mm 2 , and an intra-class correlation efficient (ICC) of 0.992 (95% CI, 0.961 -0.997). Jaccard coefficient was calculated as 0.94 ± 0.03 between the angiographic FAZs obtained by the two graders, with a false-positive error of 0.018 ± 0.012 and a false-negative error of 0.045 ± 0.028. These analyses showed a high inter-grader repeatability in angiographic FAZ delineation for healthy controls.

Theoretical baseline FAZ
Recently, Chui et al demonstrated that the inner retina had a relatively constant thickness at the boundaries of the FAZ in healthy participants [21], indicating the FAZ could be potentially obtained by inner retinal thickness. However, inner retinal thickness relies on the accurate segmentation of the OPL outer boundary, which is confounded near the fovea by the angle-dependent optical reflectivity of the Henle's fiber layer (HFL) [39]. Alternatively, Kulikov et al found the agreement between the contour of a constant GCC thickness and FAZ boundary in healthy controls [40]. These findings suggest that GCC thickness may better estimate the tbFAZ without capillary damage caused by DR since it is less prone to segmentation errors than OPL outer boundary.
We determined the GCC thickness threshold to estimate the tbFAZ boundaries on healthy controls. We measured the GCC thicknesses at the boundaries of angiographic FAZ manually delineated on OCTA ( Fig. 1(B)). The normal distribution of GCC thickness at the angiographic FAZ boundary had a maximum frequency bin of 31 μm (Fig. 2(A)) and we proposed that thickness as the criteria to define the boundaries of the theoretical undamaged FAZ of all subjects in this study. We then tested the agreement between the angiographic FAZ from OCTA and the tbFAZ determined using GCC in healthy controls. Jaccard coefficient, defined as the ratio of the intersection to the union of angiographic FAZ (A) and tbFAZ (T) (Eq. (1)), was calculated as 0.82 ± 0.059, with a false-positive error of 0.099 ± 0.068 and a false-negative error of 0.098 ± 0.078. Bland-Altman plot showed that the difference (mean ± SD) between the tbFAZ and angiographic FAZ on area was −0.001 ± 0.029 mm 2 ( Fig. 2(B)). These analyses confirmed good agreement between tbFAZ defined by GCC thickness and angiographic FAZ in healthy controls.
However, DME, characterized by intraretinal fluid and thickening of the contour, introduces potential errors in the definition of tbFAZ boundary. In mild DME cases, the cystic fluid occurs mainly in outer retina without affecting the GCC, the tbFAZ can be obtained as normal subjects. In significant DME cases, large cystic fluid pockets deform the GCC, requiring extra consideration and the corresponding strategies to obtain the tbFAZ. For eyes with significant DME (n = 13, 12.7% of total participants), the intraretinal fluid cysts were automatically detected from structural OCT by a deep-learning assisted segmentation algorithm [41,42]. We then determined the extent of affected area by projecting the automatically detected intraretinal fluid cysts on to inner retinal angiograms (Fig. 3). An edema severity scale was defined as the percentage of the area of tbFAZ with cystic fluid to the area of total tbFAZ. If the edema severity scale was smaller than 75% (n = 7, 6.9% of total participants) (Figs. 3(A1)-3(A3)), the foveal region was considered partially affected and the tbFAZ was obtained by using the unaffected portion with GCC thickness of 31 µm and fitting the rest to a circle. The final tbFAZ was then defined as the intersection of the fitted circle and the angiographic FAZ. If the edema severity scale was greater than 75% or the tbFAZ could not be defined by GCC thickness of 31 µm (n = 6, 5.9% of total participants) (Figs. 3(B1)-(B3)), a maximum inscribed circle of the angiographic FAZ was obtained as the tbFAZ. Moreover, to avoid the underestimation or overestimation of the tbFAZ, we calculated distribution of the angiographic FAZ size (97% CI: 0.198 mm 2 -0.288 mm 2 ) in the healthy control group. If the above described procedure produced a tbFAZ beyond this range, we replaced it with a circular region with equivalent area of the limit centered at the foveal center.

Para-FAZ vessel density parameters
To evaluate the foveal ischemia in DR, we analyzed a para-FAZ volume including regions between ILM and 15 µm above EZ in axial direction and annulus regions of 600 µm distance from the tbFAZ boundary in the transverse direction (Fig. 4). Para-FAZ vessel density (volume %) was characterized by a 3D vessel density (3D-PFVD, %) measurement, defined as the percentage of vascular perfused voxels to the total voxels in the analytic volume. Para-FAZ vessel density (area %) was characterized by a 2D vessel density (2D-PFVD, %) measurement in the maximum projection of the corresponding region. For patients with DME, volumetric retinal cystic fluid were excluded when calculating the 3D-PFVD.

Statistical analysis
All spatial parameters, including the angiographic FAZ area and perimeter were calculated and adjusted by individual axial length. Shapiro-Wilk was used to test normality of each data set of angiographic FAZ area, perimeter, 2D-PFVD and 3D-PFVD. Since not all groups of angiographic FAZ area and perimeter met the requirements for normality, we compared the groups using the nonparametric Kruskal-Wallis H-test. In contrast, each study group of 2D-PFVD and 3D-PFVD satisfied the conditions for the one-way analysis of variance (ANOVA) i.e. normality and equality of variances. The post-hoc tests with Bonferroni correction were used for multiple comparisons. Spearman rank-order correlation analysis tested the correlation of these four parameters with the increasing level of retinopathy severity and bestcorrected visual acuity (VA), as determined by ETDRS vision scores. A P value of less than Bonferroni adjusted threshold (0.05/6) was considered statistically significant. The sensitivity of detecting patients with DM from healthy control individuals, patients with retinopathy from those without retinopathy, and patients with severe NPDR or PDR from patients with less severe retinopathy were calculated for these four parameters. Statistical procedures were performed using SPSS 23.0 (IBM Corporation, Chicago, IL).

Results
The distributions of the obtained tbFAZ area and perimeter for the four study groups were shown in Fig. 5. The mean ± SD values for tbFAZ area were 0.244 ± 0.104 mm 2 in healthy control group, 0.229 ± 0.087 mm 2 in DM without DR group, 0.217 ± 0.085 mm 2 in mild to moderate NPDR group, and 0.197 ± 0.094 mm 2 in severe DR group. Kruskal-Wallis H-test showed no significant differences in the tbFAZ area and perimeter between the four study groups. The results indicated that the detected tbFAZ did not change with DR progression and validated its effectiveness in estimating the FAZ without the effects of DR. In contrast, the angiographic FAZ area and perimeter increased with DR severity, and were significantly larger in the eyes with DR compared to the eyes without DR. In particular, the angiographic FAZ area in the severe DR group (0.333 ± 0.167 mm 2 ) was significantly greater than the angiographic FAZ area in DM without DR group (0.227 ± 0.083 mm 2 ), but not the healthy control group (0.243 ± 0.102 mm 2 ) and mild to moderate NPDR group (0.308 ± 0.122 mm 2 ) (Fig. 6(A)). The angiographic FAZ perimeter distinguished the groups slightly better, which was significantly greater in the severe DR group (2.915 ± 1.522 mm) when compared to the healthy control (1.973 ± 0.440 mm) and DM without DR (1.981 ± 0.430 mm) groups (Fig. 6(B)), but not the mild to moderate NPDR group (2.489 ± 0.748 mm). Fig. 6. Boxplots of A: angiographic FAZ area, B: angiographic FAZ perimeter, C: 2D-PFVD and D: 3D-PFVD for the four study groups. By comparing to the Bonferroni corrected threshold (0.05/6), the significant differences between groups were annotated by asterisks.
The correlation coefficients of angiographic FAZ area, angiographic FAZ perimeter and 3D-PFVD and 2D-PFVD parameters with increasing level of retinopathy severity, bestcorrected VA and the patient's age were shown in Table 2. The sensitivities of detecting eyes with 95% specificity of healthy controls were given in Table 3. Overall, correlations with participant age were the weakest. The 3D-PFVD had the strongest correlation (rho: 0.544, P<0.01) with best-corrected VA and the strongest correlation (rho: −0.792, P<0.01) with DR severity. The 3D-PFVD had the same sensitivity of 81.3% in differentiating DM eyes compared with 2D-PFVD, however the 3D-PFVD had the highest sensitivity of 92.9% in differentiating DR eyes, and sensitivity of 96.7% in differentiating severe DR eyes, indicating an excellent diagnostic accuracy.

Discussion
Foveal ischemia in DR is associated with vision loss and retinopathy severity [15][16][17]. In this study, we have demonstrated the feasibility of 3D vessel density in para-FAZ volume to evaluate the foveal ischemia in DR. This quantification algorithm involves three steps: (1) obtain tbFAZ based on GCC thickness; (2) construct para-FAZ volume according to the detected tbFAZ; and (3) calculate projection artifact-resolved 3D vessel density. The 3D para-FAZ vessel density measurement appeared to be the most sensitive indicator for identifying eyes with DR compared to FAZ area and perimeter obtained through the angiographic FAZ, and 2D para-FAZ vessel density.
The tbFAZ defined in this way showed a good agreement with angiographic FAZ in healthy control group, and showed no significant differences with increasing DR severity (Figs. 5(A) and 5(B)). This demonstrates that constant GCC thickness could be used to estimate the tbFAZ without the effects of DR. It is worth noting that the FAZ area and perimeter vary considerably in healthy control group, resulting in substantial overlap between the healthy control and diabetics. The para-FAZ volume from the tbFAZ boundary objectively measures the capillary dropout without being affected by the normal population variations in FAZ size.
Previously, G. Lynch et al have documented that the enlargement of FAZ area by manually delineation of tbFAZ and angiographic FAZ to predict worsening of DR, but its strict exclusion criteria on cases with DME and laborious delineation of tbFAZ on en face OCT reflectance images restricted its clinical implications, making it better suited to tracking dynamic changes than to assessing the current ischemic status [8]. Furthermore, many prior studies examining foveal ischemia [5,8,10,13,43,44] excluded patients with DME. Once again, given the frequency of DME, DR quantification metrics that exclude these eyes have limited clinical utility. In this study, we included eyes with a range of DME severity. Notably, when calculating the 3D vessel density in our study, the cystic fluid volumetric voxels were excluded from the para-FAZ volume. This makes underestimation of 3D VD less likely in DME eyes due to increase in the denominator volume.
The other strengths of this study are the advanced PR-OCTA and irb-BMS algorithms used to generate 3D binary maps of vascular voxels. These approaches reduce projection artifacts and remove the decorrelation background noise, improve the reliability of axial vascular perfusion evaluation and makes volumetric vessel density quantification feasible.
Limitations of this study include ignoring the probable thinning of GCC with progression of DR, which can affect both determination of the theoretical baseline FAZ and the para-FAZ volume. A modest number of healthy eyes were used to determine the GCC thickness of 31 µm. A larger reference population may improve the threshold value. Additionally, for patients with significant DME, although we have taken into account the underestimation or overestimation of the area of tbFAZ, further investigation is needed to improve the estimation of tbFAZ in eyes with significant DME. Furthermore, the four groups were not age or sex matched. The participants in DM without DR, mild to moderate NPDR, and severe DR were older than the control group. However, we did not observe a significant age effect on these measurements ( Table 2).

Conclusions
We have developed a novel OCTA-based metric that can measure foveal ischemia more effectively. Theoretically, 3D-PFVD can be more sensitive to the FAZ enlargement introduced by vascular pathologies, because it takes into account the normal population variation in FAZ size and the 3D structure of the foveal circulation. Our study proved that 3D-PFVD is better correlated with VA and DR severity, compared to previous known parameters related to foveal ischemia. Overall, this metric shows promise as a biomarker for clinical evaluation of DR.