Highlights
- •
Ultra-widefield imaging enables monitoring of retinal vascularization rate.
- •
Slow retinal vascularization is an early predictor of treatment-requiring retinopathy of prematurity.
- •
Low-dose bevacizumab does not slow physiological retinal vascularization.
PURPOSE
To assess the rate of retinal vascularization derived from ultra-widefield (UWF) imaging-based retinopathy of prematurity (ROP) screening as predictor of type 1 ROP and characterize the effect of anti–vascular endothelial growth factor (anti-VEGF) therapy on vascularization rate.
DESIGN
Retrospective, consecutive cohort study.
SUBJECTS
The study included 132 eyes of 76 premature infants with a mean gestational age (GA) of 26.0 (±2.0 SD) weeks and birthweight (BW) of 815 (±264) g, who underwent longitudinal UWF imaging for ROP screening, at a level 3 neonatal unit in Oxford, United Kingdom.
METHODS
The extent of retinal vascularization on each UWF image was measured as the ratio between “disc-to-temporal vascular front” and “disc-to-fovea” distance along a straight line bisecting the vascular arcades. Measurements from ≥3 time points plotted against post-menstrual age (PMA) enabled calculation of temporal vascularization rate (TVR) for each eye. Using TVR, GA, and BW as predictors, a machine learning model was created to classify eyes as either group AB (no ROP and type 2 ROP) or group C (type 1 ROP). The model was validated in a withheld cohort of 32 eyes (19 infants), of which 8 eyes (5 infants) required treatment. TVR in 37 eyes (20 infants) was compared before and after ultra-low-dose (0.16 mg) intravitreal bevacizumab treatment.
MAIN OUTCOME MEASURES
The rate of retinal vascularization was determined.
RESULTS
Slower retinal vascularization correlated with increasing ROP severity, with TVR being 29% slower in group C eyes (n=50) than group AB eyes (n=33 no ROP and n=49 type 2 ROP) ( P = .04). Our model correctly predicted ROP outcomes of 30/32 eyes, achieving a balanced accuracy of 95.8%. No significant change in TVR was found before and after bevacizumab treatment with mean posttreatment imaging follow-up of 7.7 (±7.9) weeks ( P = .60 right eyes, P = .71 left eyes).
CONCLUSIONS
UWF imaging-based ROP screening enables quantification of retinal vascularization rate, which can provide early prediction of type 1 ROP independent of BW and GA. Rate of physiological retinal vascularization does not appear to be significantly affected by ultra-low-dose anti-VEGF treatment, which has significant implications for the development of peripheral avascular retina and timing of anti-VEGF intervention to prevent disease progression in high-risk infants.
R etinopathy of prematurity (ROP) is a leading cause of childhood blindness worldwide. It is characterized by aberrant retinal vascular development that progresses through 2 main phases: (1) delayed retinal vascularization in phase 1 results in peripheral avascular retina (PAR), which stimulates release of angiogenic factors, followed by (2) a vasoproliferative phase 2 associated with extraretinal neovascularization and subsequent tractional retinal detachment. ROP causes visual impairment in more than 30 000 preterm babies globally, with rising incidence over the past decades owing to improved survival of extreme preterm infants (eg, those born as early as 22 weeks).
Timely detection of ROP enables treatment by anti–vascular endothelial growth factor (anti-VEGF) therapy or laser photocoagulation to prevent sight loss. Low gestational age (GA <31 weeks) and birthweight (BW ≤1500 g) are major risk factors for developing ROP, thus constitute the inclusion criteria for ROP screening. Other risk factors include oxygen supplementation, slow postnatal weight gain, sepsis, and respiratory distress. , Although many of these factors are interlinked, their individual predictive power for ROP occurrence remains poorly defined.
The current ROP grading system (based primarily on zones I-III, stages 0-5, and normal to plus disease spectrum) is applied during ROP screening examinations using either binocular indirect ophthalmoscopy (BIO) or digital retinal imaging. BIO is a technically demanding skill in premature infants and is mainly focused on detecting features of vasoproliferation, such as stage 3 (elevated neovascular ridge) and plus disease (vessel dilation and tortuosity), which pertains to the ROP treatment threshold. However, indirect ophthalmoscopy is less suited to assessing delayed retinal vascularization, which can only be gleaned approximately from 33 to 34 weeks postmenstrual age (PMA) at which stage vessels reach the edge of zone I or zone II.
Wide-field retinal imaging could enable accurate assessment and longitudinal monitoring of the rate of retinal vascularization from birth. Contact retinal imaging such as the 130° RetCam 3 Wide-field Digital Imaging System (Natus Medical Inc., USA) provides limited ability to visualize the posterior pole and far peripheral retina in a single imaging field. Noncontact 200° ultra-widefield (UWF) retinal imaging (eg, Optomap; Optos plc., UK) could potentially overcome this limitation to enable consistent and accurate quantification of the rate of retinal vascularization over time. We have previously demonstrated the utility of UWF imaging for ROP assessment with the infants being held to the camera using the “flying baby” technique. We have more recently transitioned from BIO examination to UWF imaging as the default method of ROP screening at our level 3 neonatal unit using an Optos California imaging system mounted on a mobile trolley with onboard battery.
In this study, we aimed to devise a practical method for measuring the rate of retinal vascularization in a retrospective analysis of serial UWF images obtained from routine ROP screening. Using retinal vascularization rate, GA, and BW as predictors, we aimed to create a machine learning–based gradient-boosting model for predicting ROP requiring treatment (type 1 ROP), which demonstrates high accuracy in a separate validation data set.
The presence of PAR has been reported in a high proportion (91%-100%) of infants who have received intravitreal anti-VEGF therapy for ROP. , However, the causal relationship between anti-VEGF treatment and PAR remains unclear because PAR may be associated with ROP itself. Nonetheless, in clinical practice, concerns about potential vascular arrest following anti-VEGF treatment may lead clinicians to hold off anti-VEGF intervention in borderline treatment-requiring cases with posterior disease, for example, posterior zone II stage 2 with plus disease (UK ROP guideline ). As part of this study, we aimed to quantify and compare the rate of retinal vascularization before and after intravitreal anti-VEGF treatment for type 1 ROP.
METHODS
DATA COLLECTION AND STUDY ELIGIBILITY CRITERIA
This study was conducted as an internal retrospective clinical audit approved by the Oxford University Hospitals Integrated Governance System (Oxford University Hospitals clinical audit approval no. 9080). Because of the retrospective and anonymous nature of this audit, informed consent was waived by the ethics committee.
The ROP database at the Oxford University Hospitals NHS Foundation Trust (a level 3 neonatal care unit, John Radcliffe Hospital, Oxford, UK) was searched retrospectively from May 2019 to February 2024 to identify all infants that underwent ROP screening in accordance with UK guidelines with known outcomes and who had at least three 200° UWF retinal images (Optos California) captured at least 1 week apart. The ROP outcomes were classified into 3 clinical categories:
- •
Group A (no ROP): eyes that were recorded as having full or zone 3 retinal vascularization without developing any ROP and discharged from ROP screening;
- •
Group B (ROP not requiring treatment): eyes that had complete regression of ROP without reaching treatment threshold (ie, type 2 ROP);
- •
Group C (ROP requiring treatment): eyes with type 1 ROP that required treatment for ROP (either anti-VEGF therapy or laser).
Infants without definitive ROP outcomes (eg, infants that were transferred to other units or passed away) and those with other significant ocular comorbidity were excluded from the analysis. Eyes that did not have serial UWF imaging of sufficient quality to measure the extent of temporal retinal vascularization were also excluded.
The data set was used for the following 2 major analyses:
Analysis 1 : Retinal vascularization rates as a predictor of type 1 ROP. To develop our machine learning model for predicting type 1 ROP, we excluded retinal images obtained after any ROP treatment or after 40 weeks’ PMA to avoid bias, as fewer scans beyond 40 weeks’ PMA would be associated with infants who do not develop any ROP. Infants that started ROP screening from May 2019 up to the end of March 2023 were used to train the ROP prediction model. Infants not included in model training, who started screening from April 2023 to February 2024, were used as an independent test set to validate the predictive model.
Analysis 2 : Effects of anti-VEGF therapy on retinal vascularization. To determine the effect of anti-VEGF treatment on the rate of retinal vascularization, we included all eyes with type 1 ROP with at least 2 UWF images captured before and 2 images after intravitreal injection of bevacizumab (0.16 mg in 0.025 mL). This ultra-low-dose bevacizumab is the clinical standard at Oxford University Hospitals NHS Foundation Trust. For this analysis, UWF images captured after 40 weeks’ PMA were included to assess the long-term effects of anti-VEGF therapy on retinal vascularization.
Alongside the UWF retinal images, the following information were also collected for each infant: gestational age in weeks at birth (GA), birthweight (BW), ROP outcome based on International Classification of ROP, PMA at each ROP screening time point, and PMA at the time of anti-VEGF injection (if applicable). PMA was defined as gestational age plus postnatal age.
UWF RETINAL IMAGE ANALYSIS
All UWF retinal images used for analysis were captured using the same camera (Optos California). For each UWF retinal image, the advancement of temporal vascular front was measured as the ratio between the distance from the optic disc center to the temporal vascular front and that from the optic disc center to the fovea along a straight line that bisects the superior and inferior venous arcades. Identifying the fovea in preterm infants can be challenging because it is not as well defined owing to delayed or altered development in prematurity and ROP. ,
For the purpose of consistent measurements across longitudinal UWF images, we defined the fovea as the point that bisects the widest distance between the superior and inferior venous arcades. Retinal veins were chosen instead of arteries, as the arteries tend to become more tortuous in pre-plus or plus disease. The red-free filter was used to provide improved contrast between vessels and the surrounding tissue to aid identification of vascular front and venous arcades.
To validate this approach, we performed fovea localization using the same method on confocal scanning laser ophthalmoscopy (cSLO) images accompanying Heidelberg Flex optical coherence tomography (OCT [Spectralis; Heidelberg Engineering, Germany]) in 13 eyes of 10 premature infants ( Figure 1 ). Our protocol for Flex OCT was as previously described, with 768 A-scans and 73 B-scans over 55 × 30°. The location of “presumed” fovea on the cSLO image was compared with the “actual” anatomical fovea as identified on the corresponding OCT. This showed a mean discrepancy of only 0.4 mm (SD = 0.36) between the presumed and actual fovea locations (Supplementary Table S1), thus validating our method for consistently identifying the fovea on UWF images.
The disc-to-fovea distance (D-F) was chosen as the reference unit of distance measurement so that subsequent disc-to-temporal vascular front distance measurements could be expressed as a ratio to D-F. This is consistent with the methodology used in previous studies. , In addition, De Silva and colleagues demonstrated that longitudinal D-F measurements did not change significantly over 9 weeks in preterm and full-term infants. This normalization helps to eliminate the effects of small changes in image magnification or angle and axial length growth on calliper measurements taken from images of the same eye at different time points (Supplementary Figure S1, A). Using this method, temporal vascularization was observed to progress in a linear fashion in the untreated eyes of premature infants (Supplementary Figure S1, B).
It should be noted that although all eyes grow in size between 32 and 52 weeks’ gestation (approximately 0.15 mm/wk, as described by Cook and colleagues ), the difference in the rate of axial length increase between normal and stage 1-3 ROP eyes is minimal and, importantly, all eyes follow linear growth trajectories over the period in question. The other practical benefit of using D-F distance as the reference unit is to allow inference of vascularization to the edge of zone I, because the diameter of zone I is defined as twice the distance from the fovea to the disc center (ie, 2 × D-F) (Supplementary Figure S1, B). The location of the temporal vascular front was measured until the posterior end of the vessel. However, if a vascular ridge was present, the location of the vascular front was defined as the center (or peak) of the ridge.
All distance measurements were made using the calliper tool within OptosAdvance (Optos plc.) with in-built peripheral retinal curvature correction that assumes the axial length to be 24 mm. Therefore, we validated the use of the calliper tool in premature infants (who have axial lengths generally less than 24 mm) by looking at the disc-to-temporal vascular front distance / disc-to-fovea distance (DT/DF) ratios in untreated eyes where the vessels grew more peripheral and found no deviation in TVR that would have been expected if distance measurements became more distorted toward the periphery (Supplementary Figure S1, B).
These measurements were reviewed twice by investigator E.C. and reviewed by a second independent investigator (R.P.) to assess interrater correlation and achieve a consensus. Measurements of the disc-to-temporal vascular front were regressed against PMA (in weeks) at the time of imaging and the gradient of the line of best fit taken as the temporal vascularization rate (TVR) for an individual eye.
STATISTICAL ANALYSIS
Statistical analysis was conducted using the programming software R (v4.2.1). Linear mixed effect modeling was used to account for repeated measures using the lme4 package. Model performance and diagnostics were carried out using the performance package, and effect size analysis for the linear mixed effects model was performed using the simr package.
Temporal vascularization rate (TVR) for each eye was calculated using a linear regression slope analysis of the advancement of temporal vascular front vs PMA at imaging. As part of our preliminary work, the disc-to-nasal vascular front was also measured in a similar manner to disc-to-temporal vascular front in eyes that had UWF images that captured nasal peripheral retina on more than 1 occasion (using D-F distance as the common unit). The advancement of nasal vascular front was found to correlate strongly with temporal vascular front with a correlation coefficient of 0.91 (Supplementary Figure S2). However, there were insufficient longitudinal data to calculate the nasal vascularization rate (NVR) for most eyes because of the technical challenge of consistently imaging the peripheral nasal retina. Therefore, all subsequent analyses used TVR as a surrogate measure of the overall rate of retinal vascularization for each eye.
To compare vascular growth rates between eyes that did not require treatment and those that did require treatment, group A and B eyes were combined to represent all those with no ROP or type 2 ROP (not requiring treatment). Group C included infants with type 1 ROP for which treatment is recommended per ICROP3 guidelines. To detect differences in the rate of temporal retinal vascularization as a surrogate marker of overall vascularization rates (slopes) between these 2 groups, a linear mixed effect model with random slopes and random intercepts was fitted using both eyes and patient ID as the nesting variables to account for similarities between right and left eyes in each patient and the repeated measures nature of the longitudinal data.
The fixed effect independent predictor variables investigated were the PMA, GA, BW, and ROP groups. The R package glmulti was used to compare all possible combinations of predictors along with their interactions. Model comparisons were performed and evaluated using Akaike information criterion and Bayes information criterion to arrive at the best model to describe the independent variable of temporal vessel extent. The linear mixed effects analysis omnibus test was followed by a Tukey post hoc comparison of slopes.
Model assumptions were verified using visual normal distribution checks, q-q plot analysis, symmetry of histograms, and assessment of heteroscedasticity (Supplementary Figure S3). Where model assumptions were violated, such as in the comparison of TVR pre and post anti-VEGF treatment, robust linear mixed modeling was used. Robust models mitigate the effects of the outliers by applying a trimmed weighting to extreme outlier patients. An alternative Bayesian framework was also investigated; however the robust model fit was found to be superior, likely as a result of the presence of extreme outliers and the setting of vague priors because of a lack of existing knowledge in this novel domain.
Monte Carlo simulation adjusting effect sizes was implemented to determine the minimum effect size that could be reasonably detected at an alpha level of 0.05 and 90% power when compared to a null model.
DEVELOPMENT OF ROP CLASSIFICATION PREDICTION MODEL
Gradient boosting is an ensemble machine learning technique that generates a large number of decision trees sequentially, each learning from the last, to accurately classify according to the observed predictors. We developed a model using gradient-boosting machine learning with the predictors, BW, GA, and TVR (slope coefficient of the temporal vascularization extent vs PMA), to predict whether an infant belongs to one of 2 groups: group AB (no treatment) or group C (ROP requiring treatment). We compared 3 models for prediction: one with TVR alone, one with BW and GA, and one model with all 3 predictors.
We chose to use a gradient boosting rather than random forest model because of the presence of nested data: the right and left eyes of each patient were included in the classification model. We used the package GPBoost, which combines linear mixed effects modeling and tree boosting to train the fixed effects and assign patient ID as the random effect to account for nested eyes. , GPBoost uses cross validation to tune the hyperparameters on the training data set.
Precision-recall curve (PRC) with area under the PRC (AUPRC) were calculated in addition to the receiver operating characteristic (ROC) curves and AUROC to assess each model from the training set. Both ROC and PR curves use sensitivity but the second axis differs: ROC uses a false positive rate (1 – specificity) whereas PRC uses precision defined as how many true positives out of all that have been predicted as positives. Although AUROC is more familiar to most, the AUPRC is a more suitable measure where there is class imbalance, particularly where there are far more negative cases than positive. , The moderate imbalance between the nontreatment group (n=48) and the ROP requiring treatment group (n=28) in our data set justifies the use of AUPRC over AUROC.
Bernoulli probit likelihood distribution was used for this binary classification training model, and the optimizer employed was Nesterov-accelerated gradient descent. The gradient-boosting predictive model was then tested using a fully independent validation set and assessed using balanced accuracy.
To establish whether TVR provides predictive information above that of BW and GA, we used a variable inflation factor analysis from the ‘car’ package that assessed variable importance and the presence of collinearity. This assesses multicollinearity between the predictors (GA, BW, and TVR) by performing pairwise comparisons of every combination of predictors and analyzing each correlation coefficient obtained as a result.
For full code or information see Github link: github.com/amanasj/ROP-prediction-from-vascularisation-rates
RESULTS
DEMOGRAPHICS OF PATIENT GROUPS
A total of 1204 infants underwent ROP screening during May 2019 to February 2024. Of these, 95 infants were included in the imaging analyses. Finally, 1109 infants were excluded mainly because of gradual transition from BIO to imaging-based ROP screening, transfer to other units, and insufficient number of UWF images.
Data from 76 infants (132 eyes) that started ROP screening from May 2019 to March 2023 were included in the training data set, and so used to investigate the rate of retinal vessel growth in infants with and without ROP requiring treatment and develop a predictive model. Note that the number of eyes to infants is not at 2:1 ratio due to exclusion of eyes with insufficient quality of UWF images for the delineation of temporal vascular front. Their mean GA was 26.0 (SD ±2.0) weeks and mean BW 815 (±264) g. Of the 76 infants, 22 infants (33 eyes) who did not develop any ROP were assigned to group A; 28 infants (49 eyes) who had type 2 ROP were assigned to group B; and 28 infants (50 eyes) with type 1 ROP requiring treatment were assigned to group C. Representative UWF images for groups A, B. and C are shown in Figure 2 , A.
Patients in group A and group B were subsequently merged into a combined group AB, which represents all infants who did not require ROP treatment. Group AB consisted of 82 eyes from 48 infants. There is a discrepancy of 2 infants in the combined group, as these infants had one eye in group A and the other eye in group B. No infant in group AB received any anti-VEGF treatment or laser during the study.
The predictive model was tested using a fully independent test set containing data from 19 infants (32 eyes) that underwent ROP screening from April 2023 to February 2024 with all inclusion criteria fulfilled. Their mean GA was 26.4 (SD ±1.5) weeks and mean BW 803 (±240) g. Among these, 14 infants (24 eyes) did not require any treatment for ROP (group AB) whereas 5 infants (8 eyes) required treatment for ROP (group C). The demographics of training and test data sets are summarized in Table 1 .
| Characteristic | Training Set (n=76 Infants) | Test Set (n = 19 Infants) |
|---|---|---|
| GA, wk, mean ± SD | 26.0 ± 2.0 | 26.4 ± 1.5 |
| BW, g, mean ± SD | 815 ± 264 | 803 ± 240 |
| Infants with type 1 ROP, n (%) | 28 (37) | 5 (26) |
ANALYSIS 1
Slow rate of advancement of temporal retinal vascular front is associated with ROP
Across the training cohort, we found a trend for the TVR to decrease with increasing severity of ROP: the mean (±SEM) TVR for group A eyes was 0.15 (±0.017) disc-to-fovea distance per week (D-F/wk), group B (ROP not requiring treatment) was 0.13 (±0.013) D-F/wk, and group C (ROP requiring treatment) 0.10 (±0.013) D-F/wk ( Figure 2 , B). Omnibus analysis of variance revealed significant differences in TVR ( P = .049); however, post hoc testing with Tukey correction revealed no significant pairwise differences in TVR between group A and B ( P = .60), between group B and C ( P = .26), or between group A and C ( P = .059).
Given the clinical importance of distinguishing treatment-requiring vs non–treatment-requiring ROP and the somewhat subjective distinction between no ROP and stage 1 ROP, we merged group A and B eyes into a combined group AB (ie, all eyes that did not require any therapeutic intervention) and compared this group with group C ( Figure 2 , C). In this case, the TVR of group C was significantly slower (by 29%) than the TVR of group AB ( P = .04, TVR of group AB = 0.14±0.01 D-F/wk, group C = 0.10±0.014 D-F/wk).
Vascularization rate as a predictor of type 1 ROP
Given the significant difference in TVR between infants with and without ROP treatment, we hypothesized that TVR, calculated from at least 2 successive ROP screening visits, could be used as an independent predictor of type 1 ROP, over and above decisions based on GA and BW alone. A machine learning predictive model was developed using the training data set, comparing 3 models incorporating (1) TVR only, (2) BW+GA, or (3) TVR+BW+GA as predictors ( Figure 3 , A). The models were then tested using a fully independent test set as detailed above.
