Evaluating the impact of vitreomacular adhesion on anti-VEGF therapy for retinal vein occlusion using machine learning.

Vitreomacular adhesion (VMA) represents a prognostic biomarker in the management of exudative macular disease using anti-vascular endothelial growth factor (VEGF) agents. However, manual evaluation of VMA in 3D optical coherence tomography (OCT) is laborious and data on its impact on therapy of retinal vein occlusion (RVO) are limited. The aim of this study was to (1) develop a fully automated segmentation algorithm for the posterior vitreous boundary and (2) to study the effect of VMA on anti-VEGF therapy for RVO. A combined machine learning/graph cut segmentation algorithm for the posterior vitreous boundary was designed and evaluated. 391 patients with central/branch RVO under standardized ranibizumab treatment for 6/12 months were included in a systematic post-hoc analysis. VMA (70%) was automatically differentiated from non-VMA (30%) using the developed method combined with unsupervised clustering. In this proof-of-principle study, eyes with VMA showed larger BCVA gains than non-VMA eyes (BRVO: 15 ± 12 vs. 11 ± 11 letters, p = 0.02; CRVO: 18 ± 14 vs. 9 ± 13 letters, p < 0.01) and received a similar number of retreatments. However, this association diminished after adjustment for baseline BCVA, also when using more fine-grained VMA classes. Our study illustrates that machine learning represents a promising path to assess imaging biomarkers in OCT.

In this paper, we use fully automated computational segmentation and classification of SD-OCT scans to evaluate VMA in patients receiving ranibizumab for macular oedema secondary to branch and central RVO in two multicentre prospective clinical trials, and compare visual acuity outcomes between patients with and without VMA.

Patients and Methods
This study is a post-hoc analysis of a comprehensive clinical trial image library at the Vienna Reading Center (VRC), Vienna, Austria. For central RVO, we included SD-OCT data of patients enrolled in CRYSTAL, a phase IIIb, multicentre trial assessing the efficacy and safety of an individualized visual acuity (VA) dosing regimen of ranibizumab 0.5 mg driven by stabilization criteria 6 . For branch RVO, we included SD-OCT data of patients enrolled in BRIGHTER, a phase IIIb, multicentre study assessing the efficacy and safety of an individualized stabilization criteria-driven pro re nata (PRN) dosing regimen of ranibizumab 0.5 mg alone or in combination with laser versus laser photocoagulation 7 . Patients randomized to the laser only arm of BRIGHTER were excluded from the current analysis because of previously demonstrated inferior VA outcomes 7 . Our study was conducted in compliance with the Declaration of Helsinki. Approval was obtained by the institutional review board or ethics committee at each site participating in the multicentre trial. Furthermore, approval for the present post-hoc analysis was obtained by the Ethics Committee at the Medical University of Vienna. All participants provided written informed consent before inclusion into the CRYSTAL and BRIGHTER trials. The multicentre trials are registered with clinicaltrials.gov (identifiers NCT01599650 and NCT01535261).

Inclusion and exclusion criteria. Detailed inclusion and exclusion criteria of the CRYSTAL and
BRIGHTER studies have been reported previously 6,7 . In brief, patients with macular oedema secondary to branch or central RVO and a visual acuity between 20/40 and 20/400 (Snellen equivalent) were included. Key exclusion criteria were the use of intraocular anti-VEGF treatment ≤3 months, intraocular corticosteroid use ≤3 months, and macular laser photocoagulation ≤4 months before the screening visit. Because the VMA segmentation algorithm developed for this work (see below) is currently only available for Heidelberg Spectralis OCT images, only data of patients scanned with this particular OCT machine were included into this post-hoc analysis. Treatment schedule; functional and structural assessments. All patients analysed in this study received intravitreal 0.5 mg ranibizumab monthly until stable VA was observed for 3 consecutive monthly visits (implying a minimum of 3 consecutive treatments per protocol). The investigators then monitored the patients for VA and disease activity on a monthly basis. Monthly ranibizumab injections were reinitiated if monitoring indicated a loss of VA resulting from disease activity, and were continued until VA stabilization. In BRIGHTER, patients randomized to the ranibizumab with laser arm received macular laser based on investigator discretion as reported previously 7 .
Best-corrected visual acuity (BCVA) was measured by certified, masked examiners at each monthly visit using early treatment diabetic retinopathy study (ETDRS) charts. SD-OCT was acquired at each visit by certified, masked operators. Using Heidelberg Spectralis OCT instruments, a volume scanning pattern of 512 A-scans × 49 B-scans covering 20° × 20° centred on the fovea was acquired. Eye tracking and automated real-time averaging were activated at 29 frames. Raw data was exported and uploaded to the VRC database.
Automated segmentation of the posterior vitreous boundary. A custom designed, fully automated segmentation algorithm was used to delineate the posterior vitreous boundary in the SD-OCT volume scans. The software method is described in detail in the section "Description of segmentation algorithm". Preliminary versions of the method have been reported previously 18,19 , as well as an extension to retinal fluid and layer segmentation 20 . In brief, a machine learning algorithm was trained to classify voxels belonging to the posterior vitreous boundary in the OCT volumes. The posterior vitreous interface was regarded as a continuous surface and a 3D graph cut segmentation algorithm was used to fit an optimal surface through the voxels classified as posterior vitreous boundary with highest probability. The segmentation algorithm also delineates the internal limiting membrane (ILM). The automated segmentation method was applied on all baseline SD-OCT volume scans to delineate the posterior vitreous boundary and the ILM. The distance between the two segmentation surfaces was computed at each A-scan and presented as ILM-vitreous distance map (Fig. 1).

Automated classification of vitreomacular adhesion.
To automatically classify the ILM-vitreous distance maps of each patient into VMA and non-VMA, the histogram of the ILM-vitreous distances at each A-scan of the volume was computed. The resulting histograms were automatically clustered by unsupervised machine learning using the k-means algorithm 21 . To evaluate the robustness of the clustering, the number of clusters was systematically varied from 2 to 4. The clusters were interpreted clinically by expert evaluation of the ILM-vitreous distance map representing the cluster centre (Fig. 2). Statistical analysis. This is an exploratory post-hoc analysis of a clinical trial database. Hence, no formal hypotheses or adjustments for multiple comparisons were made, and p-values are to be viewed as hypothesis-generating. The main outcome measure of this study was the mean change in BCVA at the primary endpoint (Month 6 for BRIGHTER, Month 12 for CRYSTAL). Additional outcome measures were the change in reading-centre determined central retinal thickness (CRT) as well as the number of administered ranibizumab injections at the primary endpoint. The means of the outcome variables were compared between eyes with and without VMA using analysis of variance (ANOVA). Moreover, linear regression analysis was performed to assess the predictive impact of VMA on BCVA change, including the BCVA at baseline as an independent co-variate. The bold line in the 3D visualization corresponds to a 90° wedge being cut out of the volume; note that the line to the right of the image corresponds to the right half of the provided B-scan, while the part pointing towards the observer is perpendicular to the B-scan slice shown on the left. The linear regression analysis was performed separately for 2 clusters (VMA presence vs. VMA absence), as well as 3 and 4 clusters (VMA absence vs. different degrees of VMA). The formal significance level was set at p = 0.05.

Results
The main outcomes and patient characteristics of CRYSTAL and BRIGHTER have been published previously 6,7 . All 812 patients randomized in CRYSTAL (n = 357) and BRIGHTER (n = 455) were screened for the current analysis. Patients receiving macular laser photocoagulation only in BRIGHTER (n = 97) were excluded. Furthermore, patients imaged by OCT devices other than Spectralis OCT were excluded, resulting in 391 patients for analysis (198 in CRYSTAL;193 in BRIGHTER).

Automated segmentation and classification of vitreomacular adhesion. Automated segmen-
tation of the posterior vitreous boundary and generation of ILM-vitreous distance maps was successful in all patients. Representative segmentation results and ILM-vitreous distance maps of study eyes are provided in Fig. 1. Unsupervised clustering of the ILM-vitreous distance histograms revealed clinically distinctive classes of vitreomacular interface configurations (Fig. 2). According to the number of selected clusters (2)(3)(4), the k-means algorithm grouped non-VMA eyes; VMA with shallow vitreomacular separation and VMA with steeper adhesions Linear regression models were fitted to evaluate the interaction between baseline BCVA, VMA presence and BCVA change to the primary endpoint. In BRVO (using 2 clusters), the model fit for BCVA change was R = 0.30, p = 0.03. BCVA at baseline was a significant predictor (p = 0.03) whereas VMA presence was not (p = 0.28). For 3 clusters, the model fit was R = 0.28, p = 0.05; p = 0.01 for baseline BCVA and p = 0.74 for VMA class. For 4 clusters, the model fit was R = 0.35, p < 0.01; p < 0.01 for baseline BCVA and p = 0.05 for VMA class. Similarly, in CRVO, for 2 clusters, the model fit was R = 0.36, p = 0.01, with BCVA at baseline being a significant predictor (p = 0.02) as opposed to VMA presence (p = 0.10). Using 3 clusters, the model fit was R = 0.30, p = 0.05; p = 0.02    Tables 1 and 2. Baseline CRT did not differ among eyes with or without VMA in CRVO, while eyes without VMA showed significantly lower baseline CRT compared to eyes with VMA in BRVO.
Correspondingly, eyes with VMA demonstrated significantly larger CRT reductions at the primary endpoint in BRVO. The mean number of ranibizumab injections received up to the primary endpoint did not differ significantly between the vitreomacular interface clusters (Tables 1 and 2).

Discussion
In this study, we used fully automated computational analysis based on machine learning to characterize the condition of the vitreomacular interface in a large dataset of SD-OCT images. The proposed method demonstrated that fully automated characterization of retinal morphology is feasible in a large-scale randomized clinical trial setting, where manual analysis of SD-OCT images may not be practicable. Our results suggest that patients with VMA may show larger functional gains with individualized anti-VEGF therapy in macular oedema secondary to central and branch RVO, although this effect was not visible any longer when the analysis was adjusted for baseline BCVA. In the era of personalized medicine and the big-data situation entailed by modern OCT imaging, automated analysis of ophthalmic imaging is a promising step with the potential to transform clinical care in the future 22 .
The role of the vitreous in anti-VEGF therapy for exudative macular disease has attracted scientific interest in the recent years. While ample data exists in the setting of neovascular age-related macular degeneration, scientific studies are scarce in macular oedema of other origins. In anti-VEGF therapy of RVO, which was investigated here, Terao et al. reported superior functional and anatomical response to bevacizumab in eyes with VMA 11 . However, their study included only branch RVO and was limited in patient number. On the other hand, a recent study by Singh et al. reported no influence of VMA on outcomes of intravitreal therapy for RVO 16 . Our results show a positive impact of VMA on outcomes of anti-VEGF therapy (albeit not confirmed when adjusting for baseline BCVA) in a much larger number of patients, including those with central RVO, using data of a prospective randomized controlled trial. Other studies have investigated the role of VMA in anti-VEGF therapy of diabetic macular oedema and of macular oedema in uveitis, with similar results 15,23 .
The potential reasons for the interaction between VMA and the response to intravitreal therapy are currently poorly understood. Since the vitreous acts as a reservoir for anti-VEGF agents injected into the eye, pharmacokinetic mechanisms influencing drug half-life or diffusion in senile vitreous liquefaction or PVD may be responsible for some of the observed phenomena 12 . Although the distribution of drug molecules inside the vitreous cavity is governed by passive diffusion, the interaction between the vitreous gel and charged large molecules, as well as vitreous fluid currents may play a role in drug distribution and clearance within the eye 24,25 . Potentially, egress of VEGF molecules from the retina is altered by the presence of aqueous instead of vitreous at the vitreomacular interface 26 . Laboratory studies have demonstrated a faster clearance of anti-VEGF agents from vitrectomised eyes, which would suggest a shorter half-life of ranibizumab also in eyes with complete PVD and/or senile vitreous liquefaction 27 . A relatively prolonged clearance of anti-VEGF molecules from eyes with VMA could explain the beneficial functional and anatomical outcomes in the presence of VMA in this and other studies. However, in neovascular age-related macular degeneration, the presence of a posterior vitreous detachment (PVD) also seemed to be associated with reduced needs for retreatment, which clearly contradicts this hypothesis 12 . Recently, investigators have detected that PVD affects the intraocular cytokine levels as measured in the aqueous, which may present another contributing factor to the role of the vitreomacular interface configuration in anti-VEGF therapy 28 . Furthermore, the presence of a PVD was associated with improved oxygenation of the retina, which may likely alters VEGF expression in the eye 29 . The finding that eyes with VMA exhibited poorer mean baseline BCVA levels (at least in CRVO) may be a sign of increased VEGF load in eyes that have not yet experienced PVD. Our study adds clinical evidence supporting a role of VMA in anti-VEGF treatment of RVO, however, further experimental studies are required to shed light onto the contributing pathomechanism. Ongoing developments in advanced in-vivo imaging continue to change management paradigms in clinical practice on a frequent basis 30 . Particularly in OCT, the ever-increasing resolution and scanning speed allow to capture massive amounts of highly detailed morphologic data with each single OCT exam. Unfortunately, innovation in analysing retinal imaging data is significantly lacking behind the rapid development cycles for imaging hardware. This dissociation between hardware capabilities and limitations conferred by manual review of hundreds of images has given rise to substantial interest in automated computational analysis methods 31,32 . Validated software algorithms have the potential to rapidly and reliably process data, capture the morphology of disease in a meaningful manner, and provide clinically relevant outcome measures to the ophthalmologist.
In our study, we use a custom-developed, fully automated segmentation algorithm to detect the posterior vitreous boundary in relation to the ILM in clinical SD-OCT data. The volumetric segmentation enables examination of the vitreomacular interface configuration as a ILM-vitreous distance map as well as three-dimensional inspection (Fig. 1), which could provide a clinically useful interpretation of the posterior vitreous, e.g. prior to vitreoretinal surgery. In addition, we employ unsupervised clustering, a data-driven, hypothesis-independent machine learning technique, to automatically group individual study eyes into clinically meaningful morphological subclasses. Our clustering approach is not limited by arbitrary thresholding or other a-priori constraints. Depending on the number of selected clusters, the unsupervised k-means classifier identified clinically distinctive classes of different degrees of VMA, as well as the non-VMA class. Regardless of the number of clusters, the non-VMA class consistently demonstrated lowest BCVA gains. However, in the cohort of patients with CRVO, eyes without VMA also showed higher baseline BCVA scores, confounding these differences in BCVA gains 6 .
This study is mainly limited in its retrospective nature. Since the vitreous segmentation algorithm is currently only available for Spectralis OCT images, only 55% of the total eligible study population could be used for analysis. Potentially, the exclusion of patients imaged with other SD-OCT devices may have caused selection bias. However, the inclusion of prospectively collected, randomized controlled trial data confers many advantages, including a standardized SD-OCT acquisition protocol, treatment regimen and BCVA measurement. Our algorithm is currently being extended to also segment images of other OCT machines. As the algorithm is based on a 2D approach, feature discontinuities across different B-scans lead to a slightly reduced smoothness in the final 3D segmentation results. A feature representation in 3D would result in a smoother segmentation, however, this benefit is mitigated by the highly anisotropic voxel size of the Spectralis instrument (~3.8 × 11.1 × 118 µm). Furthermore, the non-VMA class detected by our clustering approach may contain both patients with completely attached posterior vitreous (where the ILM-vitreous distance is 0) and with complete PVD (where the vitreous boundary is antepositioned beyond the OCT scanning range and the distance is therefore unknown). In the age group of patients enrolled in the CRYSTAL and BRIGHTER studies, the prevalence of completely attached  posterior vitreous is reported to be low (14% in patients aged 60-69 years); thus, eyes with complete PVD are likely representing the majority of the non-VMA cluster 33 . Future work will address the use of Fourier-space features to investigate alternative means of introducing a degree of rotation invariant features 34,35 . Furthermore, deep learning approaches, for instance convolutional neural networks (CNN), have been shown as a promising way forward for image segmentation in retinal OCT data 36 .
In conclusion, we have demonstrated the feasibility of fully automated computational diagnosis of VMA in a large-scale clinical SD-OCT dataset. Data-driven unsupervised clustering of the segmentation results yielded clinically meaningful subclasses of non-VMA as well as different degrees of VMA. In macular oedema secondary to branch and central RVO, eyes with VMA afforded larger functional benefits from anti-VEGF therapy, although this association diminished after adjustment for baseline BCVA. Automated computational analysis of ophthalmic imaging data may help to characterize morphologic changes on OCT effectively, and may facilitate clinical management in a field of high-throughput patient care in the future.
Description of the segmentation algorithm. In order to fully automatically segment the posterior vitreous boundary, we designed an image processing pipeline consisting of four major steps: (1) feature extraction, (2) 5 class region prediction, (3) 3D graph theoretic segmentation and (4) auto-context refinement.
We treat the vitreomacular segmentation problem as a voxel-wise region probability prediction problem and define the retinal layer surfaces as the boundaries between these regions. The regions used in this algorithm are from exterior to interior the vitreous, the posterior vitreous boundary, the space between the posterior vitreous boundary and the inner limiting membrane, the retina (i.e. the space between ILM and the retinal pigment epithelium (RPE)) and the space exterior of the RPE.
The final surface segmentation is then defined as the set of four boundaries (i.e. the top and the bottom surface of the posterior vitreous boundary, the ILM and the RPE) that separate the interior from the exterior part of the OCT volume so that the corresponding regions lie above respectively below them.  . with mask[x, y] being 0 outside and 1 inside the circular patches.
By extracting sliding windows I of size 21 × 21 pixels and computing ϕ(I) for each window a map of the local orientations can be generated (Fig. 3).
Using the local image orientation, we can now extract small 21 × 21 masked image patches I and rotate them by ϕ(I). This rotated image patches are the used in a principal component analysis (PCA) to compute the first 20 eigenvectors of the patches. These eigenvectors (shown in Fig. 4) are then used as convolution kernels resulting in the image features used for the machine learning algorithm, as described in detail elsewhere 20 .
Voxelwise region prediction. After computing the voxel wise features as described above we trained a random forest classifier 38 on 337 B-scans of 88 macula-centred SD-OCT scans in which the posterior vitreous boundary was manually annotated by trained readers of the Vienna Reading Center as described previously 39 . Since only the exterior boundary of the vitreous cortex is of interest and was thus annotated we assumed the vitreous cortex region to be at least 3 voxels thick and excluded 20 voxels above the annotation from the training.
For each voxel in the training set, the random forest classifier is presented with the feature vector computed above and a label encoding which of the 5 classes the voxel belongs to (above vitreous, VMI, VMI-ILM, ILM-RPE Figure 8. Increase in segmentation accuracy due to auto-context approach. Red dots indicate the signed error of each A-scan in the test set (with random horizontal displacement for visualization purposes), white markers show quartile and median error values. and below RPE). Using this information, the classifier is able to predict the probability that a previously unseen voxel belongs to any of the 5 classes. Figure 5 illustrates an example of the probability map for the ILM-RPE class.
Surface extraction. In order to extract the surface segmentation from the generated probability maps we used a 3D graph theoretic approach described by Garvin et al. 40 . This graph cut approach requires for each voxel a cost term that represents the penalty or benefit of that voxel being below or above the final segmentation. Since the result of the classifier used in the previous step directly encodes this information, the generated probability maps can directly be used as region cost terms. For example, to segment the ILM the sum of the regions vitreous, vitreous cortex and between vitreous cortex and ILM is used as region cost in the graph theoretic approach. The dashed lines shown in Fig. 5 show the result of such a segmentation.
Iterative improvement. As shown in Fig. 5, the segmentation results (while mostly correct) contain severe segmentation errors. In order to improve the result, we apply a technique called auto-context in a similar manner as described previously 37,41 . We compute for each voxel the axial distance to all extracted surfaces and extract raw probability samples of certain classes in the close neighbourhood of the voxels (e.g. what is the probability that the voxel above me belongs to the ILM-RPE region?). These distances and samples are then merged with the image features described previously. Together they are used to train a second classifier that can leverage the information learned by the first classifier. This process can be repeated multiple times to iteratively improve the segmentation results.
The probability map and the surface segmentation after two such iterations can be seen in Fig. 6. Note that the large segmentation errors are mostly gone and the segmentation approaches the ground truth.
Validation of segmentation algorithm. The segmentation accuracy of the proposed method was tested by comparing the pixel distance between the ILM and the posterior vitreous boundary as determined by the algorithm against the manual annotation. For this purpose, 11 SD-OCT volumes (with a total of 26,189 annotated A-scans) were randomly chosen from the data set and excluded from the training set used in the steps described above. For each A-scan the difference between the automatic and the manual segmentation was calculated (Fig. 7).
As shown in Fig. 8, the iterative improvement method described above increases the segmentation accuracy significantly after only a few iterations. However, due to the higher computation costs and the diminishing increase in segmentation accuracy, only few iterations are practicable.