Pointwise visual field estimation from optical coherence tomography in glaucoma: a structure-function analysis using deep learning

Background/Aims: Standard Automated Perimetry (SAP) is the gold standard to monitor visual field (VF) loss in glaucoma management, but is prone to intra-subject variability. We developed and validated a deep learning (DL) regression model that estimates pointwise and overall VF loss from unsegmented optical coherence tomography (OCT) scans. Methods: Eight DL regression models were trained with various retinal imaging modalities: circumpapillary OCT at 3.5mm, 4.1mm, 4.7mm diameter, and scanning laser ophthalmoscopy (SLO) en face images to estimate mean deviation (MD) and 52 threshold values. This retrospective study used data from patients who underwent a complete glaucoma examination, including a reliable Humphrey Field Analyzer (HFA) 24-2 SITA Standard VF exam and a SPECTRALIS OCT scan using the Glaucoma Module Premium Edition. Results: A total of 1378 matched OCT-VF pairs of 496 patients (863 eyes) were included for training and evaluation of the DL models. Average sample MD was -7.53dB (from -33.8dB to +2.0dB). For 52 VF threshold values’ estimation, the circumpapillary OCT scan with the largest radius (4.7mm) achieved the best performance among all individual models (Pearson r=0.77, 95% CI=[0.72-0.82]). For MD, prediction averaging of OCT-trained models (3.5mm, 4.1mm, 4.7mm) resulted in a Pearson r of 0.78 [0.73-0.83] on the validation set and comparable performance on the test set (Pearson r=0.79 [0.75-0.82]). Conclusion: DL on unsegmented OCT scans accurately predicts pointwise and mean deviation of 242 VF in glaucoma patients. Automated VF from unsegmented OCT could be a solution for patients unable to produce reliable perimetry results. Introduction Glaucoma causes retinal ganglion cell (RGC) loss, resulting in structural and functional changes in the visual system. Standard automated perimetry (SAP) is the reference technique to follow functional visual field (VF) loss during glaucoma management.1,2 Current SAP devices such as the Humphrey® Field Analyzer (HFA) have high intra-subject variability and a lengthy examination time.3,4 Retinal nerve fiber layer (RNFL) thickness measurements in circumpapillary optical coherence tomography (OCT) scans are standard practice for the quantitative assessment of structural retinal damage and RGC loss during glaucoma progression.5 The structure-function relationship in glaucoma is an intensively studied topic.6–9 Functional damage only becomes noticeable with current SAP when extensive RGC loss has occurred.10 SAP detects progressive glaucomatous damage better than OCT when VF loss has already occurred.11 Studies that use OCT to estimate visual field report limited accuracy. Most models rely on assumptions such as log transformation to predict decibel (dB) VF values from retinal layer thickness values.12–14 Deep learning (DL) approaches overcome these constraints because they can model non-linear functions.15–19 Still, most current DL developments have limits because they (1) estimate global VF indices such as mean deviation (MD)15,16, or (2) use OCT-derived information such as RNFL thickness maps that require a priori retinal layer segmentation17–19. This study goes beyond state-of-the-art by using unsegmented SPECTRALIS® OCT scans to develop and evaluate DL models that estimate the visual field sensitivity threshold at each location (52 threshold values) and MD as measured by the HFA. Materials and Methods Data initially extracted comprised 1643 matched OCT-VF pairs corresponding to the 998 eyes of 542 patients who visited the University Hospitals Leuven's glaucoma clinic between 2015-2019. This work is part of the larger study on “automated glaucoma detection with deep learning” (study number S60649), approved by the Ethics Committee Research UZ/KU Leuven in November 2017. Informed consent was waived due to the retrospective nature and patient reidentification made impossible as the link between patient ID and study ID was deliberately removed. The research adhered to the tenets of the Declaration of Helsinki. Inclusion criteria were: (1) the availability of a SPECTRALIS® OCT (Heidelberg Engineering, Heidelberg, Germany) scan using the Glaucoma Module Premium Edition (GMPE), containing one scanning laser ophthalmoscopy (SLO) en face image, 24 radial scans and three circumpapillary rings, and (2) the results of an HFA3 exam with the strategy 24-2 SITA Standard (54 test points) obtained with the Humphrey® Field Analyzer (model 850 v1.3.1.2, Carl Zeiss Meditec, Dublin, CA, USA). We used the circumpapillary rings as these scans cover all nerve fibers that pass through the optic nerve, unlike the radial scans. The SLO image served as intra-study benchmark. The 1643 OCT-VF pairs of 542 patients were allocated to train, validation, and test set, accounting for 60/20/20 of the patients, respectively. We took care that all data from a single patient were stored in the same partition to avoid leakage. Subsequently, VF data of the validation and test data were filtered on standard HFA reliability indices, with false positives (FP), false negatives (FN), and fixation losses (FL) not exceeding 15%, 33%, and 20%, respectively.20 We did not filter unreliable visual field data in the training set but used it as data augmentation (noise). OCT data were extracted in RAW format using Heyex v6.12.1 software (Heidelberg Engineering, Heidelberg, Germany). The binary files were subsequently processed with heyexReader v0.1.3, a Python package for reading Heyex OCT files. The three circumpapillary RNFL rings (3.5mm, 4.1mm, and 4.7mm) and SLO were extracted as lossless image files with dimensions 768x496 and 1536x1536, respectively. We obtained VF data with HFA3 that were analyzed in PeriData v3.5.7 (PeriData Software GmbH, Hürth, Germany). Pointwise sensitivity threshold values were extracted from the individual patients' printouts using an optical character recognition (OCR) tool developed for this task. OCR output was manually verified on 10% of the data, matching perfectly with actual threshold values. These values were paired with global indices such as MD that were exported as a comma-separated value text file by PeriData. Two VF test points were discarded in all analyses, as these are located on the anatomical blind spot (see Figure 1), resulting in 52 threshold values to model. Eight DL models were trained using the 3.5mm, 4.1mm, 4.7mm circumpapillary rings and SLO images, for both MD and 52 threshold values. We compared the single models with ensemble (averaged) predictions. We selected the Xception21 architecture pretrained on ImageNet22. The convolutional neural network (CNN) was followed by a global average pooling and convolution operation, avoiding any fully connected layers to minimize overfitting. The final convolution operation had either one or 52 filters, depending on the target (MD or 52 threshold values), featuring a linear activation to allow for regression. Models were trained using mean squared error (MSE) loss, optimized using Adam23 with a starting learning rate of 1e-4. The latter was reduced to 75% of its value after ten epochs without improving the validation loss. Each epoch featured 300 training steps of batches containing four preprocessed and augmented images. Single-channel input images were either upscaled to 768x512 (circumpapillary scans) or downscaled to 512x512 (SLO), with intensity values rescaled between 0 and 1. Augmentation included horizontal flipping, elastic deformation, and cutout24. Model development and evaluation were performed using Keras25 v2.2.4, Tensorflow26 v1.12.0, in a Python 3.6.7 environment running on a server with six GTX 1080 Ti and two TITAN V graphics processing units. The coefficient of determination (R2), Pearson's r (r), and mean absolute error (MAE) were the metrics to evaluate model performances. R2 and MAE were computed using the scikit-learn library27, and r using the NumPy library28. The best model configuration (one for MD, one for threshold values) was selected on the highest R2 metric and evaluated with the independent test set. We computed a baseline MAE for validation and test sets by predicting the mean MD or threshold value. 95% confidence intervals (CI) were obtained through bootstrap sampling (5000 iterations). The model trained on SLO images offered an objective way to assess the added value of OCT (with complete retinal layer information) over en face retinal imaging. We visualized correlation per individual VF point and VF sector (Garway-Heath29) to verify which VF regions are better modeled than others. Finally, we compared the 90% CI of VF sensitivity threshold predictions against empirically established test-retest variability published by Artes et al.3 Results Study sample Study sample characteristics are presented in Table 1. Of note, the average MD was -7.53dB (ranging from -33.8dB to +2.0dB), which is expected considering that the data come from a tertiary glaucoma clinic. After filtering on standard reliability indices, 1378 OCT-VF pairs were eligible for model development and evaluation. The baseline MAE values for MD and threshold estimation were 7.15dB and 8.17dB for the validation set and slightly lower for the test set (6.26dB and 7.64dB, respectively, see last column of Table 2). MD estimation The CNN model with an Xception encoder explained up to 71% (95% CI: 0.62-0.78) of the variance in the validation set of 190 OCT-VF pairs, equivalent to an average correlation coefficient of 84%. The MAE of 3.33dB obtained using the 4.7mm scan reduced the baseline MAE by 53%. The performances of single models were similar between circumpapillary scans, but significantly outperformed the model trained with SLO images, with the latter explaining 46% (0.27-0.55) of the variance (r=0.68). The averaging of predictions of the four single models resulted in improved evaluation metrics: R2=0.74 (+0.03), r=0.87 (+0.02), MAE = 3.21 dB (-0.12). This ensemble model had similar results on the independent test set