Artificial intelligence as a gateway to scientific discovery: Uncovering features in retinal fundus images

In this work, we describe a novel methodological pipeline that leverages deep learning to reveal retinal features not currently part of diagnostician’s toolkit within the standard of care. This proposed pipeline is composed of three phases. In Phase 1, we train a CNN model to classify retinal fundus images based on a given binary label. Here, we have used patient sex (male vs. female) as a case study, though the methodology can be extended to other traits and diagnostic labels. The set of fundoscopic images used to train and evaluate this model (the CNN Development Set) is described in detail in Section 3.2. In Sections 3.3–3.6 we specify the CNN architecture as well as the details of the training and evaluation procedure. In Phase 2, we utilize CNN explainability tools to generate some insight as to what regions, anatomical landmarks, retinal features and signs visible in the fundoscopic images contributed to the successful classification of the trained model. These explainability tools are described in detail in Section 3.7. These observations are then used to enumerate a liberally construed list of exploratory hypotheses. These hypotheses are articulated as potential differences in various retinal parameters which were tested on the Exploration Set (see Section 3.2 for a detailed description). The various techniques to extract and measure the retinal parameters have been described in Section 3.8. Once the exploratory hypotheses are statistically evaluated based on the Exploration Set, a select set of candidate hypotheses, i.e., the ones that yielded significant differences in Phase 2 were then passed on to Phase 3 for a rigorous statistical verification on an independent set of fundoscopic images–the Verification Set (described in detail in Section 3.2). Phase 3 statistically re-evaluates and attempts to verify whether positive results obtained in Phase 2 are replicated when tested on the Verification Set. The details of this final analysis is described in Section 3.9.

Our methodology relies on forming three mutually-disjoint datasets: the CNN Development Set (which is further partitioned into training, validation and test sets, as commonly done in neural network training), the Exploration Set, and the Verification Set.

We used two different sources of retinal fundus images in this study. The first, ODIR [31] is a publicly available fundus dataset containing 7,000 annotated images from 3,500 patients with information regarding their age, sex, and pathological condition of each eye. In this dataset, samples are labeled for diabetic retinopathy, glaucoma, cataract, hypertension, AMD, myopia, and other diseases or abnormalities. Images are collected from different hospitals and medical centers and therefore captured by a variety of fundus cameras with different resolutions and angles of view. To prevent the model from using potential sex differences in the prevalence of diseases, we excluded all images with any eye disease or other abnormalities, as well as images with low quality, as identified by a trained ophthalmologist (Dr. Ozturan). The resulting dataset consists of 3146 images from 1991 patients, all of which contributed to the CNN development set. The summary statistics are reported in Table 1. The second source for fundus images was the retinal imaging database of Vancouver General Hospital Ophthalmic Imaging Department, named “VCH-source” from here on. It comprises a total of 2100 images from 1167 patients. These images are selected from individuals with healthy eyes, assessed by Dr. Ozturan to exclude the eyes with signs of ocular disease or abnormality. 1600 images (of 800 patients) of the VCH source was labeled as the “VCH phase-1” set, and contributed to the CNN development set. Table 2 represents the summary statistics of the VCH phase-1 dataset. In sum, a total of 4746 images (ODIR + VCH-phase 1) were placed in the CNN development set and thus used to train and assess the CNN model. Both ODIR and VCH phase-1 portions of the development set were randomly partitioned into training, validation, and test sets with approximately 70%, 15%, and 15% of the aggregate sets, respectively. Care was taken to make sure all images of any given patient remained in the same partition. All images were cropped to get a square image with equal height and width by detecting the circular contour of the fundus images and placing it at the center of the square.

The remaining 500 images from the VCH source were used to populate the Exploration set (100 images, phase-2) and the Verification set (400 images, phase-3). In both sets, exactly half of the male and female images are taken from the right eye. Age did not differ between male and female patients in both sets (p $>$ .2). The summary statistics are reported in Table 3.

The study was approved by the UBC Clinical Research Ethics Board and Vancouver Coastal Health Research Institute and the requirement of consent was waived (H21-02013).

The model architecture used in this study is Inception-Resnet-v2 [32]. The model’s parameters were initialized by the pretrained model weights on the ImageNet dataset [33]. Since the pretrained model has 1000 output classes consistent with the ImageNet classification contest, the model’s classifier module, which is the last fully connected layer, was replaced with a new randomly initialized fully connected layer containing 1536 inputs (the number of output features of InceptionResNet-v2 model) and two outputs corresponding to male and female classes.

We utilized a transfer learning approach followed by a fine-tuning step to train the network using the training subset of the CNN development set. During the first 20 epochs, the network’s weights were frozen while the new classifier layer was learning the task. This is a common technique to prevent the gradient calculated based on the initial random weights of the classifier layer from changing the network’s parameters in a direction that is not meaningful and not necessarily aligned with the task. At the conclusion of the first 20 epochs, as the classifier has learned the task to some extent, we unfroze the network’s weights and allowed them to change during the subsequent 100 epochs. Hyperparameters were tuned based on the validation performance and by trying various combinations. A summary of the hyperparameters used for training and evaluating the model can be found in Table 4.

To further improve the performance of the model, we took advantage of data augmentation and image transforms. During the training process all images were rotated by a random amount chosen uniformly from -10 to +10 degrees to prevent the network from memorizing image-label pairs. Furthermore, we utilized a novel idea, not used before in similar studies to the best of our knowledge, that can be applied to fundus image datasets specifically because of their nature: The left and right retinas are anatomical mirror-images, approximately symmetrical along vertical axis, leading to a large image-level change (left vs. right) that is not related or informative to the sex classification task. We removed this image variance across the data set by horizontally flipping all images taken from the right-eye so they appear like a left-eye fundus image. This “horizontal flipping” transform removes part of the image variance across the dataset that is known a priori to be irrelevant to the task and improves the model performance for sex classification. Arguably this is because horizontal flipping allows the model to expect the same anatomical parts of the retina in nearly same locations of the input image (i.e., optic disc on the left side and fovea on the right side) and, in turn, learns the features more efficiently.

The validation subsets of the CNN development dataset was used for evaluating the Model during the training process and tuning the hyperparameters. At every training epoch, various performance metrics were calculated and recorded on both training and validation sets: area under the receiver operating characteristic curve (AUC) score, accuracy, hit rate, false alarm rate, and binary cross-entropy (BCE) loss. Once the training course is over, the epoch at which the highest validation AUC occurred is selected to report validation metrics. In addition, the Model’s weights from the same epoch were saved as “the best model’s weights”. Then, the best Model was reloaded to obtain the performance metrics separately on the two unseen test sets obtained from ODIR and VCH phase-1 to assess generalizability performance.

We used non-parametric bootstrapping to evaluate the significance of the results. We generated B = 1000 bootstrap replicates of the test sets to obtain the confidence interval for each performance metric. The chance level calculated from the ratio of male images in each test set was then compared to the AUC confidence intervals. p-values were calculated based on the percentile rank of chance-level (50%) performance in the bootstrap AUC distribution.

In order to maintain a benchmark performance for comparison, three models are trained independently using the identical procedure described in Section 3.4, namely, “Normal Model”, “Flip Model”, and “Random Model”. The Flip Model is trained with the horizontal flip of right-eye images as explained earlier in Section 3.5, while this transform is not used for training the Normal Model. Random Model was trained on the same datasets with the only difference that male and female labels were randomly shuffled in advance to training. Performance metrics on the test sets are reported for the three models.

We used the Grad-CAM [19] technique to generate saliency maps. In this method, input images were first normalized by the average and SD values of each color channel calculated based on the ImageNet dataset (see 5), the data on which the pretrained models are trained. Each image was then fed to the network to complete the forward path needed for calculating the gradient during the backward process, and the predicted label was saved. The model’s output was one-hot coded, i.e., the output corresponding to the predicted class was set to one and the other class to zero. Next, the Grad-CAM saliency map was generated by back-propagating the gradient of the predicted label to the last convolutional layer. In an independent process, the Guided Backpropagation map was also calculated by the deconvolutional network created as an inverse of the trained model (see Ref. 18 for details) and again setting the predicted class’s output to one and the other class’s output to zero. Finally, these two matrices were multiplied pixel-wise to form the Guided Grad-CAM saliency maps.

To visualize the preferred stimuli for male and female classes, we use the regularized activation maximization method [20]. This technique allows the user to input a sample image and provides a transformed version that maximizes the response at some layer of the network. Here, we input both noise images as well as male and female fundoscopic images from the test sets. Similarly with the saliency map method, input images are normalized before going through the optimization process, and the model’s parameters are fixed. Stochastic Gradient Descent (SGD) was used as the method to maximize the activation of the output class of interest. We tried different combinations of hyperparameters in the implementation of Ref. 34 and selected the values that led to qualitatively more interpretable results. The hyperparameters used to generate the final results are summarized in Table 5.

We used the PyTorch implementation by Ref. 34 for both Guided Grad-CAM and regularized feature visualization. The main challenge in visualizing CNN models is their complexity; therefore, a simpler model was used in this phase. A VGG16 [35] model was trained using the same procedure explained in Section 3.4 and achieved test AUC scores of 0.66 and 0.73 on ODIR and VCH phase-1 datasets, respectively. This model was used for interpretation purposes as it resulted in better and more understandable visualizations because of its simpler architecture compared to Inception-Resnet-v2.

In order to test the exploratory hypotheses, a wide variety of retinal parameters were measured. Thus, we needed to segment the main anatomical parts of the retina, specifically, the optic disc, retinal vasculature, and the peripapillary area. The segmentation results were then used to quantitatively measure the variables of interest and statistically test the hypothesized sex differences.

The optic disc and the fovea were segmented using the Sefexa software [36] under the supervision of an ophthalmologist (Dr. Ozturan). A sample segmentation output is depicted in Fig. 1. The optic disc mask covers the non-parametric oval-like shape of the optic disc, and the fovea mask is a small dot marking the approximate location of the fovea since a precise location cannot be determined in this imaging modality. To segment the vasculature, we trained the LadderNet model [37], a deep learning architecture specialized for image segmentation, using the DRIVE dataset. The implementation by Ref. 38 was used for training the model on the DRIVE dataset [39] and then applying the trained model to the datasets used in this study. See the supplementary materials for the details of the vessel segmentation procedure and the DRIVE dataset.

Based on the segmentation masks, multiple measurements were made to statistically test the hypotheses derived from CNN interpretation results. To characterize the optic disc (OD), we measured the area, sharpness of the edge, and brightness. The area covered by the vessels in different parts of the retina, as well as the number of nodes, number of branches, and the total length of branches in the vessel graph, were measured to characterize vasculature. Also, the radius of the foveal avascular zone (FAZ) was estimated. In order to account for the variation in the angular field of view, measurements were normalized by the distance between the fovea and the center of the optic disc, where appropriate. This is the same normalization method used by ophthalmologists; for instance, to assess the changes in optic disc size, they measure its diameter as a fraction of the distance between the fovea and optic disc.

In the definition of the measurements and variables, $G$ denotes the original fundus image with size $w\times h$ converted to gray-scale by using cvtColor function from OpenCV Python library. $G_{x,y}$ denotes the value of the pixel located at column x and row y of the gray-scale fundus image. In the binary masks (for the optic disc, vessels, and field of view (FOV)), each pixel can take the value of either zero (not annotated as the region of interest) or one (annotated as the region of interest). $OD_{x,y}$, $V_{x,y}$, and $FOV_{x,y}$ represent the value of pixel (x,y) of the optic disc mask, the vessel mask, and the FOV mask respectively.
Center of optic disc
The optic disc center is calculated as the average coordinates of pixels included in the optic disc mask:

Then the gray-scale image is passed through the Sobel operator to obtain the derivative of the image (denoted as $S$) and is averaged by the edge mask.

A two-sample t-test (male eyes vs. female eyes) has been performed for each of the retinal measurements. For the short-listed hypotheses (phase-3 ”verification”), the p-values were then adjusted [40] using the Benjamini-Hochberg method [41] to address the multiple comparisons problem.

In order to track the training progress, the network’s performance was measured real-time over the epochs on both training and validation sets. Three evaluation metrics, namely AUC, accuracy, and BCE loss, are plotted separately for the training (blue) and validation (red) sets along with the learning rate in Fig. 2. The vertical black line indicates the epoch at which the highest validation AUC was achieved, and the model was saved for reporting metrics on the unseen test sets. The colored horizontal lines show the value of each metric at the best epoch on the training and validation sets. The vertical gray line shows the epoch at which the model is unfrozen, and the whole network is trained. The learning rate also resets to its initial value at this point. By looking at the patterns of changes in the performance metrics, we ensured that the model is learning the task properly and that the training course is running favorably. As shown in Fig. 2, AUC and accuracy, target indications of classification performance, improve over the training epochs on both training and validation sets. Although, the loss shows an increasing behavior on the test set a few epoch after unfreezing the network, this pattern has been attributed to the nature of BCE loss (see Ref. 22 for details).

The three models were tested on the unseen ODIR and VCH phase-1 test partitions. To obtain the significance level of the results and calculate p-values, the performance metrics achieved by each model were compared to the chance level accuracy, which was 54.16% and 50% on ODIR and VCH phase-1 datasets, respectively (see Table 6).

On the unseen ODIR test partition, the AUC and accuracy achieved by the Normal Model were 0.645 and 0.625, respectively, both significantly greater than chance level ($p^{\prime}s<.001$). As anticipated, Flip Model performed slightly better than the Normal Model obtaining an AUC score of 0.676 ($p<.001$) and an accuracy of 0.646 ($p<.001$). The random model achieved 0.480 AUC ($p=.992$) and 0.542 ($p=.477$) accuracy, which are not significantly different from and very close to the chance level, as expected. The 95% confidence intervals (CIs) are calculated based on non-parametric bootstrapping with $B=1000$ bootstrap resamples of the ODIR test set with size $N=480$. All metrics, p-values, and confidence intervals obtained from ODIR by the three models are shown in Table 6.a. The results imply that the network classifies patient sex significantly higher than chance level at an $\alpha=.05$ confidence level.

Based on the unseen VCH phase-1 test partition, the Normal Model obtained 0.760 AUC ($p<.001$) and 0.692 accuracy ($p<.001$). The Flip Model achieved higher performance on this dataset as well by an AUC of 0.785 ($p<.001$) and accuracy of 0.721 ($p<.001$). The Random model again performed similar to chance level with 0.422 AUC ($p=.979$) and 0.496 accuracy ($p=.476$). $B=1000$ bootstrap resamples of the test set with $N=240$ sample size were used to calculate the confidence intervals. The AUC score, accuracy, p-values, and CI reported in Tabel 6.b indicates that the network has obtained a sex classification performance significantly higher than the chance level.

Based on the observations on the saliency map and feature visualization results, possible sex differences in the retina were described as general exploratory hypotheses over three main themes: (1) retinal vessels are more prominent in males than females, (2) optic disc is more prominent in males than females, (3) the peripapillary area is darker in males than females.

According to the two-sample t-test results, out of the 14 retinal parameters that were tested on the exploratory dataset, 9 were significantly different between males and females. The normalized area of optic disc is significantly larger in males. Female retinas show significantly brighter peripapillary areas compared to males. Regarding vessel properties, male eyes have higher vessel coverage in the superior temporal quadrant, inferior temporal quadrant, and macula. Also, the FAZ normalized radius is greater in males than in females. All three vessel graph properties, i.e., number of nodes, number of branches, and total length of branches, are significantly higher in males. On the other hand, we did not detect any significant sex differences in the normalized brightness of optic disc and the sharpness of optic disc edge. Also the vessel coverage in the superior nasal quadrant, inferior nasal quadrant, and macula did not show significant differences between the two groups. The average and standard deviation of all tested parameters for males and females, along with the p-values and effect sizes, are reported in Table 7.

The verification results appear in Table 8 in which the average and SD values for male and female groups, along with the effect sizes, p-values, and BH-adjusted p-values, are reported. Out of the nine parameters that resulted in significant sex differences on the exploration dataset, five again showed significant results in the expected direction on the verification dataset. Male fundus images showed more nodes and branches in the vessel graph as well as a greater total length of branches. Also, the area covered by the vessels in the superior temporal quadrant of the retina was higher in men. These results confirm that retinal vasculature is more prominent in males compared to females. However, no significant differences were observed in vessel coverage of the other quadrants of the retina and in the FAZ radius. Female fundi showed significantly brighter peripapillary areas aligned with what was observed in the CNN interpretation results. There was no significant difference in the normalized area of the optic disc.

DRIVE [39] is one of the most used benchmark datasets for retinal vessel segmentation, and numerous papers have published the performance of their automated vessel segmentation method on this dataset. DRIVE contains 40 samples which are originally separated into two sets of equal size of 20 as the training and test sets. Each sample is a pair of a fundus image and a black-and-white binary mask of the retinal vessels as shown Fig. S1.

Based on Ref. 38 implementation, images were cropped into multiple smaller and overlapping patches to enhance the dataset size as it has a limited number of images. The patch size was set to 64 and the stride to 16, which resulted in 150,000 image patches. Some sample training patches, along with their ground truth vessel masks, are shown in Fig. S2. The hyperparameters are summarized in Table S1.

Then, to take advantage of the entire DRIVE dataset and enhance the performance, the same parameters were used to train the model using both training and test sets. The resulting model was applied to the VCH dataset for obtaining final vessel segmentation masks. Also, because of the domain adaptation issues, the VCH images were first preprocessed before being fed to the model. The preprocessing procedure involved increasing contrast by a factor of 1.5, decreasing brightness by a factor of 0.6, and slightly blurring the images with a Gaussian blur filter. These steps were taken to make VCH images get closer to DRIVE images that the model was trained on, and results were assessed qualitatively by try and error to choose the best preprocessing hyperparameters.

After training, the model achieved a test AUC score of 0.987 and an F1 score of 0.853. Fig. S3 shows sample vessel segmentation results from the test set.