Decoding systemic vascular health and hypertensive disorders in pregnancy through retinal imaging and Visionary AI

Study design, ethics and participants

This prospective study was approved by the Institutional Review Boards of CU and NYU and adhered to the tenets of the Declaration of Helsinki. Informed consent was obtained from every participant before enrollment. The study was designed to evaluate whether retinal vascular features captured during early pregnancy could predict subsequent HDPs. CU served as the development cohort, and NYU served as an independent external validation cohort.

CU development cohort

Pregnant individuals receiving prenatal care at NewYork-Presbyterian/CU Irving Medical Center between 2021 and 2025 were enrolled in the CU development cohort (n = 1,267). Eligible participants were at least 18 years of age, able to provide informed consent and could be enrolled through the second trimester of pregnancy.

Maternal demographic and clinical data were extracted from the electronic health record, including age, race, ethnicity, obstetric history, medical history, medication exposures, laboratory values and pregnancy-related diagnoses (that is, gestational diabetes mellitus and HDPs). Delivery outcomes were also recorded, including birth weight, small for gestational age, severe small for gestational age, stillbirth and Apgar scores.

Retinal images were obtained from each eye by trained research technicians using the Optos ultrawidefield Primary scanning laser ophthalmoscope, which captures 200° images covering approximately 80% of the retina. Images were acquired during first-trimester, second-trimester and third-trimester study visits when available. Because the primary objective was early prediction before clinical diagnosis, CU model development analyses used first-trimester images when available. For participants without first-trimester imaging, second-trimester images obtained before 20 weeks of gestation were used. History of ocular conditions and prior ocular surgery was captured and used to determine imaging-related exclusions as described below.

NYU external validation cohort

Pregnant individuals receiving prenatal care at the Mignone Women’s Health Collaborative at NYU Langone Health between 2025 and 2026 were enrolled in the NYU external validation cohort (n = 79). Eligible participants were at least 18 years of age and able to provide informed consent. Maternal demographic, clinical and pregnancy-outcome data were extracted from the electronic health record using the same general categories as in the CU cohort.

Retinal images were obtained from each eye by trained research technicians using the Optos ultrawidefield California scanning laser ophthalmoscope, which also captures 200° images covering approximately 80% of the retina. The NYU cohort differed from the CU cohort in clinical workflow, imaging-device configuration, race and ethnicity distribution and comorbidity profile. For the external validation analysis, only retinal images acquired during the first trimester, until 13 weeks of gestation were used. Thus, the NYU validation analysis provided an early-pregnancy external test of the CU-trained model under a different site, population, imaging workflow and device configuration.

Extended Data Tables 1 and 2 summarize demographic, clinical and pregnancy-outcome characteristics for the CU and NYU cohorts.

Clinical outcomes and case definitions

In both cohorts, HDPs were diagnosed by participants’ obstetric clinicians during routine prenatal care and extracted from the electronic health record. Preeclampsia was defined according to the 2013 American College of Obstetricians and Gynecologists criteria⁸. GHTN was defined as new-onset hypertension after 20 weeks of gestation without diagnostic criteria for preeclampsia. CHTN was defined as hypertension present before pregnancy or diagnosed before 20 weeks of gestation.

Preeclampsia cases were further classified by disease severity and timing of onset. SF was defined as preeclampsia accompanied by at least one severe clinical feature, including severe-range BP (systolic BP ≥ 160 mmHg or diastolic BP ≥ 110 mmHg on two occasions more than 4 h apart), thrombocytopenia, impaired liver function, renal insufficiency, pulmonary edema or new-onset neurologic symptoms⁸. NSF met diagnostic criteria for preeclampsia but did not meet criteria for severe features (systolic BP ≥ 140 mmHg or diastolic BP ≥ 90 mmHg on two occasions that were at least 4 h apart, accompanied by either proteinuria or mild signs of maternal organ dysfunction)⁸. EOPE was defined as symptoms at or before 34 weeks of gestation and LOPE was defined as symptoms after 34 weeks of gestation.

In the CU development cohort, 55 participants were diagnosed with preeclampsia (corresponding to 4.3% observed preeclampsia prevalence), including 38 with severe features and 17 without severe features. In total, 20 cases were classified as EOPE and 35 were classified as LOPE. All CU preeclampsia cases were independently reviewed by a maternal–fetal medicine specialist who were masked to retinal imaging results. One CU participant with preeclampsia was excluded from retinal-model analyses because of poor retinal image quality; therefore, retinal-model analyses included 54 CU preeclampsia cases where indicated.

In the NYU external validation cohort, 14 participants were diagnosed with preeclampsia, corresponding to an observed validation cohort preeclampsia prevalence of approximately 17.7%. Preeclampsia cases included nine with severe features and five without severe features; one case was classified as EOPE and 13 cases were classified as LOPE. All NYU preeclampsia cases were independently reviewed by a maternal–fetal medicine specialist who was masked to retinal imaging results. One NYU participant with preeclampsia was excluded from the retinal-model evaluation analyses because of poor retinal image quality; therefore, retinal-model analyses included 13 preeclampsia cases were indicated.

Birth outcomes were extracted from the electronic health record, including gestational age at delivery, birth weight, small for gestational age, severe small for gestational age, stillbirth and Apgar scores. Small for gestational age was defined as birth weight below the tenth percentile for gestational age and severe small for gestational age was defined as birth weight below the third percentile.

Retinal image acquisition and preprocessing

Ultrawidefield retinal images were acquired by trained research technicians at each study site. In the CU development cohort, images were acquired using the Optos ultrawidefield Primary scanning laser ophthalmoscope. In the NYU external validation cohort, images were acquired using the Optos ultrawidefield California scanning laser ophthalmoscope. Both systems capture 200° retinal images covering approximately 80% of the retina. After image-level preprocessing and quality control, the analytic retinal image dataset included 4,361 images in the CU and NYU cohorts.

For the CU development analyses, retinal images acquired during the first trimester were used when available. For participants without first-trimester imaging, second-trimester images obtained before 20 weeks of gestation were used. This approach was selected to maximize early-pregnancy sample size while preserving the presymptomatic prediction framework. In contrast, the NYU external validation analysis was restricted to retinal images acquired until 13 weeks of gestation, only during the first trimester.

Artifact mitigation began at image acquisition. During clinical imaging, imagers were instructed to acquire multiple images when eyelashes, eyelids, shadows, positioning differences or other localized artifacts were observed. This multi-image acquisition strategy allowed vascular regions obscured in one image to be recovered from another image of the same eye during VSI generation and merging, reducing sensitivity to localized artifacts and limiting reliance on manual exclusion of lower-quality images.

Before vessel segmentation, retinal images were standardized to support consistent downstream processing. Each image was then cropped using a standardized elliptical mask to remove nonretinal regions and peripheral artifacts outside the retinal field of view, including eyelashes, eyelids and acquisition-related obstructions. This preprocessing step was applied programmatically and uniformly across images before AI-based vessel segmentation. The elliptical mask was 3,904 × 3,008 pixels and captured approximately 50% of the image where the majority of the retinal vasculature is captured. Images were rescaled to a uniform zoom level and resized to 912 × 912 pixels.

Images with inadequate quality for vessel segmentation or downstream feature extraction were excluded according to prespecified quality-control procedures described below. History of ocular conditions or prior ocular surgery was also captured and used to guide image-level or participant-level exclusions where relevant.

AI-based vessel segmentation and VSI generation

Preprocessed ultrawidefield retinal images were converted into VSIs using a dedicated retinal vessel segmentation AI model developed for this study. VSIs provide binary representations of the retinal vascular tree and serve as the input for downstream vascular feature extraction.

The segmentation model was inspired by retinal vessel generative adversarial network (RV-GAN) architectures and used a U-Net-based generator designed to capture both large-scale vascular structure and fine vessel detail⁴⁹. The RV-GAN model was initially trained using publicly available 45° color fundus-image datasets paired with expert-annotated vessel segmentation maps^{50,51,52,53,54}. It was then fine-tuned on a curated set of 20 ultrawidefield retinal images from the study cohort with corresponding hand-annotated vessel maps. These ultrawidefield vessel maps were reviewed by a retinal specialist to ensure anatomical plausibility and segmentation quality.

The segmentation model generated vessel-probability maps, in which each pixel was assigned a probability of representing retinal vasculature. Probability maps were binarized using a threshold of 0.3 to preserve fine vascular detail and small disconnected components were removed during postprocessing. The resulting VSIs were used for downstream VSI merging and feature extraction.

When multiple images were available from the same eye, VSIs were merged to generate a more complete representation of the retinal vasculature and to recover regions obscured by localized artifacts in individual images (Extended Data Fig. 6). VSI merging was performed sequentially at the eye level. First, oriented fast rotated BRIEF (ORB) keypoint detection and descriptor matching were used to identify corresponding vascular structures across VSIs. A homography transformation was then estimated using random sample consensus to achieve coarse alignment while reducing the influence of mismatched keypoints. Next, a gradient-descent image-registration algorithm was used to identify the local transformation that maximized vessel overlap between the two VSIs being merged. This deformable registration step locally adjusted the images to account for nonlinear differences in image acquisition and retinal geometry. When more than two VSIs were available for the same eye, images were merged sequentially according to acquisition order.

After eye-level VSI generation, vascular features were extracted independently from the right and left eyes and when both were available averaged to generate participant-level feature values. One-eye imaging data were retained when image quality was sufficient for vessel segmentation and downstream feature extraction. In the CU cohort, 29 subjects only had data for one eye, while two subjects had data for only one eye in the NYU cohort.

Image and VSI quality control and analytic image or VSI exclusions

Image and VSI quality were evaluated using three complementary metrics designed to capture vascular completeness and common acquisition artifacts: vessel-density score, eyelash score and artifact score (Supplementary Fig. 12 and Supplementary Note 4). These metrics were used to assess whether retinal images and derived VSIs were suitable for downstream vascular feature extraction and model evaluation.

The vessel-density score was calculated from each VSI to quantify the spatial completeness of the detected retinal vasculature. Each VSI was divided into a 12 × 12 grid and the proportion of grid squares containing vessel pixels above a prespecified threshold was calculated. A grid square was considered vessel-containing if vessel pixels exceeded 1% of the square. This score captured incomplete vascular coverage that could arise from poor image quality, obstruction or segmentation failure.

The eyelash score was generated using an ordinal image-quality classifier trained to identify eyelash obstruction in fundus images. Images were assigned an ordinal score corresponding to no, mild, moderate or severe eyelash obstruction. This metric was used to quantify localized obscuration that could affect vessel visibility and downstream feature extraction.

The artifact score was generated using a binary classifier trained to identify acquisition-related artifacts, including machine artifacts, shadows, eyelids and other image-quality defects that could interfere with segmentation or vascular feature extraction. Images assigned high artifact burden were reviewed according to prespecified quality-control criteria.

Quality-control metrics were evaluated at the image, eye and participant levels.

When multiple images were available for the same eye, localized artifacts in individual images were mitigated through VSI merging, as described above. Images or participant-level vascular representations were excluded from retinal-model analyses only when image quality was insufficient for reliable vessel segmentation or downstream feature extraction, according to prespecified segmentation quality criteria. These criteria included sparse vessel segmentation, defined by a vessel-density score < 0.25, and additional segmentation quality failures, described in Supplementary Note 5. Participants were excluded only when adequate vascular representations could not be obtained from either eye. Using these criteria, four participants were excluded from the performance-optimized CU analysis, six from the stability-optimized CU analysis with expanded controls and 13 from the NYU cohort. Sensitivity analyses evaluating robustness to image-quality variation and single-eye availability are described under model evaluation procedures, with detailed results provided in Supplementary Note 5.

Retinal vascular feature generation

Retinal vascular features were generated from participant-level VSIs. Feature extraction was designed to quantify complementary dimensions of retinal microvascular architecture, including graph topology, vessel geometry, vascular complexity and spatial organization and hierarchical loop structure (nesting tree). Detailed feature definitions are provided in Supplementary Note 1.

VSIs were first converted into graph-based representations of the retinal vascular tree. Binary vessel maps were skeletonized and nodes were assigned to vessel bifurcations and terminal points. Edges represented vessel segments connecting adjacent nodes. From these graphs, we extracted features describing branching architecture and network connectivity, including numbers of nodes, bifurcation points, terminal points, connected components, internode distances and topological length. Topological length was calculated as a graph-based measure of hierarchical vessel depth, with both weighted and unweighted versions computed, and was evaluated separately in superior, inferior and whole-retina regions when applicable.

Geometric features were extracted to quantify vessel shape, caliber and curvature. These included vessel-length and vessel-density measures, vessel-thickness features and multiple tortuosity measures. Tortuosity features summarized deviations from a straight vessel path using complementary metrics, including sinuosity, inflection-based tortuosity, curvature-based tortuosity and tortuosity density.

To capture higher-order vascular complexity and spatial organization, we extracted box-counting, fractal and TDA features. Box-counting features quantified how retinal vessels occupied space across scales and dimensionality of high-dimensional box-counting outputs was reduced using principal component analysis before model training. TDA was used to quantify connected components and loops across multiple filtrations, including radial and distance-based filtrations. Persistence diagrams were converted into vectorized representations suitable for downstream predictive modeling, with principal component analysis used where appropriate to reduce dimensionality.

We also extracted nesting-tree features to quantify the hierarchical organization of vascular loops and bifurcations. These features summarized loop count, loop organization, asymmetry and redundancy of the vascular network, providing complementary information about network resilience and hierarchical structure.

For selected feature classes, distributional-shift features were calculated by comparing each participant’s empirical feature distribution to distributions observed in healthy controls using two-sample Kolmogorov–Smirnov statistics. These distributional features were calculated within the cross-validation framework to prevent information leakage; when a healthy control served as the held-out participant, that participant was excluded from the reference set used to construct training features.

Extracted features were grouped into 32 semantically related feature sets spanning graph topology, geometry, complexity and organization, nesting-tree structure and mixed feature categories. These feature sets served as inputs to independent base learners in the Visionary AI model architecture described below.

Visionary AI model architecture

Visionary AI was designed as a stacked ensemble framework that integrates multiple biologically defined retinal vascular feature sets into a participant-level risk prediction. The model architecture consisted of two levels: feature-set-specific base learners and a metalearner that aggregated base-learner predictions.

Each of the 32 retinal vascular feature sets described above was used to train independent base learners. For each feature set, we evaluated three model classes: logistic regression (LR), random forest (RF) and extreme gradient boosting (XGB). This generated an initial candidate space of 96 retinal vascular base learners. Each base learner was trained on a single semantically defined feature set, allowing model predictions to remain interpretable at the level of vascular feature sets.

Throughout the retinal-model analyses, we included two minimal obstetric-history adjustment variables: first pregnancy and prior preeclampsia history (Supplementary Table 7). These variables were included to account for key pregnancy-history context while preserving the primary retinal vascular focus of the model. No broader clinical risk model was included in the primary Visionary AI architecture (Supplementary Note 6).

Base-learner predictions were integrated using a stacked generalization framework. In the CU development cohort, LOOCV was used to generate out-of-fold predictions for each participant. For each fold, base learners were trained on all participants except the held-out participant and predictions for the held-out participant were stored as out-of-fold base-learner predictions. These out-of-fold predictions were then used as inputs to an XGB metalearner, which generated the final participant-level risk prediction.

RFE was used within the metalearner framework to select a subset of informative base learners and reduce redundancy across correlated vascular representations. This approach allowed Visionary AI to integrate complementary vascular signals while limiting reliance on an unnecessarily large ensemble. All model fitting, base-learner prediction generation, metalearner training and feature-selection steps were performed within the appropriate training folds to mitigate overfitting.

For participants with features available from both eyes, right-eye and left-eye feature values were averaged before model training to generate a participant-level retinal vascular representation. Intereye differences may arise from biological asymmetry, acquisition variability, eyelash or eyelid obstruction, localized artifacts or segmentation incompleteness. Bilateral averaging was, therefore, used to reduce sensitivity to any single-eye artifact or local segmentation error and to better capture systemic retinal vascular architecture. When only one eye was available after quality control, the available eye was used. Final model outputs were participant-level.

Model development and evaluation settings

Visionary AI was evaluated across a series of increasingly stringent model development and validation settings. These included a high-contrast CU analysis comparing preeclampsia cases to narrowly defined healthy controls, a performance-optimized CU population-wide analysis, a stability-optimized CU model designed to prioritize reproducibility across control definitions and an independent NYU external validation analysis performed without retraining or optimization using NYU outcome labels.

High-contrast CU analysis

The high-contrast CU analysis was designed to determine whether retinal vascular features captured early in pregnancy contained detectable signal associated with subsequent preeclampsia. In this setting, preeclampsia cases were compared to a narrowly defined healthy-control group selected to minimize medical, ocular, medication-related and pregnancy-related conditions that could independently influence retinal vascular structure (n = 136, 54 cases and 82 healthy controls).

Healthy controls were drawn from singleton pregnancies without documented preexisting maternal comorbidities, ocular conditions, relevant medication exposures or pregnancy complications. Exclusion criteria included cardiometabolic, neurologic, vascular, hematologic, endocrine, infectious or genetic diseases, prior ocular surgery or trauma, multifetal gestation, use of medications that could reflect or modify vascular risk, including antihypertensive agents, insulin or aspirin, and pregnancy complications, including hypertensive disorders, gestational diabetes, cholestasis, hyperemesis, stillbirth, multifetal gestation or smoking-related exposures. From the eligible healthy-control pool, a subset of 82 controls was selected to support class balance and approximate matching on key demographic and imaging variables.

Performance-optimized CU population-wide analysis

To evaluate model performance in a more clinically heterogeneous obstetric population, we performed a population-wide CU analysis in which preeclampsia cases were compared to broader nonpreeclampsia controls from the CU cohort (preeclampsia: n = 1,137, including 54 preeclampsia cases, 82 healthy controls and 1,001 population-wide controls; GHTN: n = 1,132, including 49 GHTN cases, 82 healthy controls and 1,001 population-wide controls; CHTN: n = 1,106, including 61 CHTN cases, 82 healthy controls and 963 population-wide controls). Unlike the high-contrast analysis, the population-wide control pool retained participants with common clinical comorbidities, including CHTN, diabetes and obesity, to better reflect the heterogeneity encountered in prenatal care.

Exclusions from the population-wide control pool were limited to prespecified conditions that could confound outcome definition or preclude valid retinal-model evaluation, including multifetal gestation, HELLP syndrome, other maternal hypertensive disorders, stillbirth and unavailable or low-quality retinal images. For both preeclampsia and GHTN, CHTN cases were not excluded from the control set. For GHTN, preeclampsia cases were excluded and not considered as controls. For CHTN, both GHTN and preeclampsia cases were excluded and not considered as controls. For the performance-optimized model, repeated sampled control sets were generated from the broader population-wide control pool and the selected healthy controls were included in each sampled evaluation set. Specifically, each population-wide ensemble model was trained and evaluated using the 54 preeclampsia cases, the 82 selected healthy controls and one of ten randomly sampled population-wide noncase sets, each containing at least 100 controls sampled without replacement from the broader population-wide control pool. Performance metrics were aggregated across repeated population-wide control samplings. Because these sampled evaluation sets did not reflect the full underlying CU preeclampsia prevalence, PPV and PR metrics calculated directly from the sampled sets were interpreted as sampled evaluation metrics rather than real-world screening estimates. Prevalence-adjusted PPV and NPV were, therefore, calculated separately as described below.

The performance-optimized model used the stacked ensemble architecture described above to integrate base learners trained on complementary retinal vascular feature sets. Base learners were constructed by training LR, RF and XGB classifiers on each feature set. Hyperparameters were tuned using LOOCV, with the best combination selected on the basis of F₁ score. The LR grid included penalty {L1, L2}, C {0.01, 0.1, 1.0, 10} and solver {saga}. The RF grid included number of estimators {100, 200}, max depth {1, 3, none}, minimum samples per leaf {1, 2}, minimum samples per split {2, 5} and max features {sqrt, 0.5, none}. The XGB grid included number of estimators {100, 250}, max depth {1, 2, 3}, learning rate {0.01, 0.1, 0.3}, subsample {0.8, 1.0} and scale-positive weight {2, 8, 12}. Base learners with AUC ≤ 0.5 were excluded from metalearner consideration.

Out-of-fold predictions from retained base learners were used to train an XGB metalearner. Metalearner hyperparameters were tuned using LOOCV to prioritize F₁ score. The metalearner grid included number of estimators {100, 250}, max depth {1, 2, 3}, learning rate {0.01, 0.1, 0.3}, subsample {0.8, 1.0}, scale-positive weight {2, 8, 12} and number of base learners selected by recursive feature elimination (RFE) {8, 16, 32}. RFE was used to select a subset of informative base learners and reduce redundancy among correlated retinal vascular representations. For GHTN, the learning rate was limited to 0.01 and the same base-learner parameters were selected for all metalearners on the basis of average performance.

Stability-optimized CU model

To reduce model complexity and prioritize reproducible retinal vascular signal, we developed a stability-optimized CU model (preeclampsia: n = 1,188, including 54 preeclampsia cases and 1,134 population-wide controls; GHTN: n = 1,138, including 49 GHTN cases and 1,089 population-wide controls; CHTN: n = 1,143, including 61 CHTN cases and 1,082 population-wide controls). For the stability-optimized model, control exclusion criteria were similar to those used for the performance-optimized model, with one modification intended to better reflect the clinical heterogeneity of routine obstetric care. Specifically, participants with HDPs (for example, GHTN) that were not diagnosed with preeclampsia, were retained in the control set. Thus, for GHTN and CHTN, preeclampsia was excluded from their control sets but GHTN was included in the preeclampsia and CHTN control set. This design allowed the stability-optimized preeclampsia model to be evaluated against a broader and more clinically realistic spectrum of nonpreeclampsia pregnancies. This model was designed to favor consistency across nonoverlapping population-wide control settings rather than maximal performance in a single site configuration.

For this analysis, repeated use of the same healthy controls across control settings was removed. Each healthy control was included only once, together with an expanded population-wide control pool, and the resulting controls were divided into five nonoverlapping population-wide control settings.

To construct the five population-wide groups, we partitioned eligible controls into five mutually exclusive sets—PW1 through PW5—with similar distributions of prespecified clinical characteristics and image-quality measures. We first identified all healthy controls, as defined above, and divided them into five equally sized groups. We then grouped the remaining controls according to shared profiles across the following characteristics: history of preeclampsia, GHTN or gestational diabetes, multigravidity, obesity, advanced maternal age, gestational anemia, IVF, cardiac disease, history of hypertension and eyelash and artifact ranks derived from the image-quality models described below. Participants with the same profile across these characteristics were initially assigned to the same stratum. Strata with fewer than five participants were merged with the most similar stratum on the basis of the Hamming distance between their binarized characteristic profiles. The controls were then allocated across PW1–PW5 on the basis of the resulting strata. This procedure produced more consistent distributions of key clinical and image-quality characteristics across the five population-wide control groups.

Candidate base learners and hyperparameter configurations were evaluated across these population-wide settings. Hyperparameter configurations were selected on the basis of low variability in AP and high minimum AUC across control settings, thereby prioritizing base learners with consistent performance across population definitions.

Specifically, the base-learner search space consisted of 96 candidate base learners representing 32 retinal vascular feature sets and three model types: LR, RF and XGB. LOOCV was used to train base learners within each of the five CU population-wide control groups. The LR grid included penalty {L1, L2}, C {0.01, 0.1, 1.0, 10} and solver {saga}. The RF grid included number of estimators {100, 200}, max depth {1, 3, none}, minimum samples per leaf {1, 2}, minimum samples per split {2, 5} and max features {sqrt, 0.5, none}. The XGB grid included number of estimators {100, 250}, max depth {1, 2, 3}, learning rate {0.01, 0.1, 0.3}, subsample {0.8, 1.0} and scale-positive weight {4, 5, 6}.

Performance of each base learner and hyperparameter combination was calculated within each population-wide control group using AUC and AP. Hyperparameter combinations with all zero feature importance or coefficients for at least one control group were excluded. For each base learner, hyperparameter combinations were then summarized across the five control groups by calculating the minimum AUC and the s.d. of AP. Hyperparameter combinations were filtered to retain those with AP s.d. below the 25th percentile for that base learner (Supplementary Fig. 1). Among the remaining combinations, the configuration with the highest minimum AUC was selected as the stability-optimized hyperparameter configuration. Base learners were retained for metalearner training only if the minimum AUC of the stability-selected configuration exceeded 0.5 (Extended Data Fig. 2).

During metalearner training, out-of-fold probabilities from retained base learners were used for hyperparameter tuning and RFE. The XGB metalearner grid included number of estimators {100, 250}, max depth {1, 2, 3}, learning rate {0.01, 0.1, 0.3}, subsample {0.8, 1.0}, scale-positive weight {2, 8, 12} and number of base learners selected by RFE {8, 16, 32}. The number of retained base learners was tuned during metalearner training and used to perform RFE, thereby limiting model complexity and removing redundant base learners.

For preeclampsia, the highest AUC achieved during metalearner hyperparameter tuning across all five population-wide control groups used eight selected base learners. Because different LOOCV folds could select different sets of eight base learners, final model construction required an additional deduplication step. For each population-wide control group, we considered the union of base learners selected across LOOCV folds and selected the final eight base learners on the basis of the average metalearner base-learner importance, while disallowing base learners trained on repeated feature sets to reduce collinearity and correlated signal amplification. For GHTN and CHTN, some population-wide control groups selected more than eight base learners; in these cases, up to the number of base learners selected during hyperparameter tuning were retained, depending on the number remaining after feature-set deduplication.

These steps ensured that the final stability-optimized model used for validation was substantially more constrained than the initial candidate search space and retained base learners reflected recurrent retinal vascular representations across the five population-wide control groups. The resulting stability-optimized CU model was used for external validation.

NYU external validation

The stability-optimized CU model was evaluated in the independent NYU validation cohort without retraining, refitting, feature reselection or hyperparameter optimization using NYU outcome labels. The CU-derived model architecture selected base learners, hyperparameters and learned retinal vascular representations were preserved. NYU validation was restricted to retinal images acquired up until 13 weeks of gestation, only during the first trimester.

Unlike the CU population-wide analyses, the NYU validation analysis did not rely on repeated or balanced and sampled control sets. Cases and controls in the NYU analytic validation cohort were evaluated together at their observed prevalence. Because the observed preeclampsia prevalence in the NYU validation cohort was higher than the annual institutional preeclampsia prevalence at NYU Langone Health, PPV and NPV were reported both at the observed validation cohort prevalence and after adjustment to the institutional annual prevalence estimate.

To reduce site-associated and device-associated feature-scale differences before applying the CU-trained model to NYU, we applied a prespecified robust feature distribution harmonization procedure to the NYU features. This preprocessing step was applied uniformly to all NYU validation samples and did not use individual-level NYU outcome labels or perform model fitting, feature selection, hyperparameter tuning or threshold optimization on the NYU data. CU reference distributions were constructed by repeatedly sampling CU controls without replacement and combining them with CU preeclampsia cases to match the aggregate case fraction of the NYU validation cohort. For each feature, the median and interquartile range (IQR) were calculated within and averaged across each CU reference set. NYU features were then robustly rescaled by centering the distribution (that is, subtracting the NYU median), dividing by the NYU IQR, multiplying by the corresponding CU reference IQR and shifting by the CU reference median. Medians and IQRs were used rather than means and s.d. to reduce sensitivity to outliers. Because this harmonization used only aggregate validation cohort composition and did not use individual-level NYU outcome labels or NYU labels for any model-fitting decision, we treat it as a feature distribution preprocessing step rather than as model retraining, threshold optimization or probability calibration.

Preeclampsia subtype analyses

Preeclampsia subtype analyses were performed as post hoc stratified evaluations of models trained on overall preeclampsia. The model was not retrained separately for SF, NSF, EOPE or LOPE. For each subtype, model predictions from the overall preeclampsia model were evaluated among participants belonging to the corresponding subtype group. These analyses were used to assess whether the overall preeclampsia model retained predictive signal across clinically heterogeneous preeclampsia presentations.

Differential-diagnosis analyses

To assess whether Visionary AI distinguished preeclampsia from related HDPs, differential-diagnosis analyses compared preeclampsia cases with participants diagnosed with other hypertensive disorders, including GHTN and CHTN. These analyses were intended to evaluate whether the preeclampsia-associated retinal vascular signal was separable from retinal vascular patterns associated with other hypertensive pregnancy phenotypes.

Subgroup analyses

Subgroup analyses were conducted to evaluate model behavior across demographic and clinical risk factor strata when sample sizes permitted. Subgroups included race and ethnicity categories, obesity status, CHTN and other clinically relevant risk factor groups available in the electronic health record. These analyses were interpreted descriptively because several strata contained limited numbers of preeclampsia cases.

Sensitivity and robustness analyses

Sensitivity and robustness analyses were conducted to evaluate whether model performance depended on image quality, eye availability or segmentation quality. Image-quality sensitivity analyses considered vessel-density, eyelash and artifact scores. Technical details for image-quality metrics, vessel segmentation and VSI generation are provided above and in Supplementary Note 4; model evaluation sensitivity analyses were interpreted as robustness checks rather than independent validation cohorts.

To evaluate whether model performance depended on bilateral feature averaging, we performed a single-eye sensitivity analysis using the final Visionary AI model. Rather than averaging features across both eyes, we selected one eye per participant and generated predictions using single-eye feature summaries. To assess whether single-eye performance was influenced by image quality, eye selection was based on the vessel-density score. Specifically, we evaluated model performance when selecting, for each participant, either the eye with the highest or lowest vessel-density score, eyelash score and artifact score. This analysis tested whether prediction was robust to single-eye selection and whether performance was disproportionately affected by eyes with lower vascular density or less complete vascular coverage.

Bootstrap resampling analysis

To quantify uncertainty in model performance, we performed participant-level stratified bootstrap resampling of the NYU external validation cohort. In each of 10,000 replicates, cases and controls were sampled separately with replacement, preserving the original cohort size and case–control distribution: 13 cases and 53 controls in the analytic NYU cohort and eight cases and 25 controls in the aspirin-recommended subgroup. All images from a participant were retained together within each replicate and image-level predictions were aggregated to the participant level before metric calculation. We calculated AUC, AP and threshold-based performance metrics for each replicate and derived nonparametric 95% confidence intervals from the 2.5th and 97.5th percentiles of the resulting bootstrap distributions. For threshold-based analyses, the model-specific operating thresholds were prespecified on the basis of clinical considerations, targeting a FPR of no more than 10%, and were fixed before evaluation of the NYU cohort. These thresholds were applied unchanged in every bootstrap replicate (Supplementary Note 7).

Clinical benchmarking analyses

Visionary AI was benchmarked against established clinical risk-stratification approaches for preeclampsia, including the FMF risk calculator and aspirin-eligibility criteria used in routine prenatal care. These analyses were designed to evaluate whether retinal vascular features provided predictive information beyond clinical risk factors and pregnancy-history variables.

FMF risk calculator

The official FMF²⁷ first-trimester risk calculator portal was used to estimate preeclampsia risk from maternal demographic, clinical and obstetric-history variables. The NYU cohort provided the most complete FMF comparison because all required maternal demographic, clinical and family-history variables were available from the electronic health record. These included maternal age, height, weight, racial origin, CHTN, type 1 diabetes, smoking status, systemic lupus erythematosus, method of conception, obstetric history and family history of preeclampsia.

FMF risk estimates were also generated for the CU cohort when available variables permitted. In the CU cohort, major maternal demographic and clinical-history variables were available but family history of preeclampsia, height, weight and a direct systemic lupus erythematosus variable were not available. When applicable, rheumatologic disease was used as a proxy for systemic lupus erythematosus. Because of these missing or proxy variables, CU FMF comparisons were treated as secondary benchmarking analyses and interpreted with caution.

FMF-derived probabilities were compared to Visionary AI predictions in the CU performance-optimized analysis, CU stability-optimized analysis and NYU external validation cohort. In addition, we evaluated a combined Visionary AI + FMF model in which Visionary AI predictions and FMF probabilities were integrated into a final prediction model. This analysis was used to assess whether FMF-derived clinical risk estimates added predictive information beyond the retinal vascular signal captured by Visionary AI.

NYU aspirin-eligibility benchmark

We also compared Visionary AI to the clinical risk-stratification approach used at NYU and CU to guide aspirin eligibility at approximately 12 weeks of gestation. Aspirin eligibility was abstracted from the electronic health record and reflected routine prenatal-care decision making on the basis of clinical risk factors, including IVF, maternal age over 35 years, BMI over 30 and other clinician-assessed risk factors used in standard obstetric care according to ACOG guidelines. Aspirin eligibility was treated as a binary clinical benchmark and compared to Visionary AI using threshold-based operating characteristics, including FPR, TPR, PPV and NPV.

Aspirin-recommended subgroup analysis

As an additional sensitivity analysis, we evaluated Visionary AI within the subgroup of NYU participants who were recommended aspirin during pregnancy. This subgroup represented a clinically enriched population expected to have elevated baseline risk on the basis of routine obstetric assessment. The stability-optimized CU model was applied to this subgroup without retraining or optimization using NYU outcome labels. Performance within the aspirin-recommended subgroup was compared to FMF-derived risk estimates to assess whether retinal vascular features retained predictive signal among individuals already identified as higher risk by clinical criteria.

Prevalence-adjusted PPV and NPV

Because PPV and NPV depend on outcome prevalence, threshold-based predictive values were interpreted in the context of each evaluation setting. In CU population-wide analyses, repeated sampled control sets were used and did not reflect the full underlying CU preeclampsia prevalence. Therefore, PPV and NPV estimated directly from sampled CU evaluation sets were not interpreted as real-world screening estimates. Prevalence-adjusted PPV and NPV were calculated using the observed CU preeclampsia prevalence of 4%.

In the NYU external validation cohort, cases and controls were evaluated together without repeated or balanced control sampling. However, the observed preeclampsia prevalence in the analytic NYU validation cohort was higher than the annual institutional preeclampsia prevalence at NYU Langone Health. Therefore, NYU threshold-based PPV and NPV were reported both at the observed validation cohort prevalence and after adjustment to the institutional annual prevalence estimate of 8%.

For prevalence adjustment, PPV and NPV were calculated from sensitivity, specificity, and the target prevalence using standard diagnostic-test formulas:

These prevalence-adjusted estimates were used to provide clinically contextualized screening performance estimates across baseline-risk settings.

RETFound and Inception-ResNet-v2 deep-learning benchmarks

To benchmark Visionary AI against direct image-level deep-learning approaches, we compared its performance to RETFound, a retinal foundation model, and Inception-ResNet-v2, a CNN model. Both models were trained and evaluated in the CU development cohort using participant-level LOOCV. In each LOOCV iteration, all retinal images from the held-out participant were excluded from model training, preprocessing decisions and model selection.

RETFound was initialized with publicly available pretrained weights and fine-tuned for preeclampsia prediction for up to 25 epochs using the developers’ publicly available scripts and default fine-tuning parameters. The model generated image-level probabilities, which were averaged across all available images from each participant to obtain participant-level predictions.

Inception-ResNet-v2 was evaluated as a transfer-learning CNN baseline and trained for up to 25 epochs. All layers were frozen except the final classification layer; thus, a new classification head was trained on fixed pretrained image representations without end-to-end fine-tuning of the convolutional backbone. Training used a weighted cross-entropy loss, with weights of 1 for controls and 4 for cases, reflecting the approximately 20% case prevalence in each combined population-wide training set. Image-level probabilities were averaged across all available images from each participant to obtain participant-level predictions.

For external evaluation, the CU-trained RETFound and Inception-ResNet-v2 models were applied to the independent NYU cohort without using NYU outcomes for model fitting, fine-tuning, preprocessing or model selection decisions, hyperparameter or epoch selection and calibration.

Feature-importance and interpretability analyses

Feature-importance analyses were performed to identify the retinal vascular representations contributing to Visionary AI predictions. Importance was evaluated at multiple levels of the stacked ensemble: individual features within base learners, base learners within the metalearner and aggregations of base learners into broader vascular feature categories.

Feature importance was extracted directly from RF and XGB models. For LR base learners, the absolute value of the model coefficients was used as the feature-importance measure. Feature-importance values were normalized within each model before aggregation. To aggregate importance across leave-one-out folds, normalized feature-importance values were averaged across folds; features or base learners not selected in a given fold were assigned an importance of zero for that fold.

For base learners, feature importance quantified the contribution of each vascular trait in a feature set. Similarly, metalearner importance was calculated from the trained metalearner to quantify the contribution of each base learner to the final participant-level prediction. For the importance heat maps, the base-learner importance was calculated by aggregating the importance of base learners trained on the same vascular feature set across model classes (LR, XGB and RF). Category-level importance was calculated by aggregating base-learner importance across broader retinal vascular categories, including graph topology, geometry, complexity and organization, nesting-tree features and mixed features.

For the stability-optimized model, feature-importance consistency was evaluated across leave-one-out folds and across nonoverlapping population-wide control settings. Recurrent base learners and feature categories were identified by examining which were repeatedly selected and assigned high importance across control settings. These analyses were used to assess whether model predictions were driven by stable retinal vascular representations rather than fold-specific or control-set-specific feature-selection artifacts.

Univariate vascular feature analyses

Univariate analyses were performed to characterize individual retinal vascular features associated with preeclampsia and related HDPs. These analyses were used for biological interpretation and hypothesis generation, not the primary basis for predictive performance claims, which were based on the multivariate Visionary AI models and independent validation analyses.

For each selected vascular feature, distributions were compared between cases and control groups using nonparametric statistical tests. Mann–Whitney U-tests were used to evaluate differences in feature distributions between groups, with emphasis on differences in central tendency. Kolmogorov–Smirnov tests were used to evaluate broader distributional differences between groups, including differences in shape, spread or cumulative distribution patterns.

For features represented as subject-level summaries, we also calculated effect-size measures where applicable, including odds ratios and fold-change risk ratios. LR models were used to estimate the association between individual vascular features and outcome status when appropriate. These analyses were performed for comparisons between preeclampsia cases and healthy controls, between preeclampsia cases and population-wide controls and, where sample size permitted, across preeclampsia subtypes and other HDPs.

Because many vascular features and disease comparisons were evaluated, Benjamini–Hochberg false discovery rate (FDR) correction was applied to the univariate vascular feature analyses. Both nominal P values and FDR-adjusted P values were reported where applicable. Features that did not remain significant after FDR correction were interpreted cautiously and treated as hypothesis generating, with emphasis placed on effect size, consistency of direction, convergence across feature categories and agreement with multivariate feature-importance analyses.

Univariate analyses were interpreted in conjunction with model-based feature importance. This approach allowed us to assess whether retinal vascular feature categories retained by Visionary AI corresponded to measurable differences in individual vascular features, while avoiding overinterpretation of any single univariate association.

Model performance evaluation

Model performance was evaluated using participant-level predictions. For cross-validated CU analyses, out-of-fold predictions from LOOCV were used to calculate performance metrics. For the NYU external validation cohort, predictions were generated by applying the CU-trained stability-optimized model to NYU participants without retraining, refitting, feature reselection or hyperparameter optimization using NYU outcome labels.

Primary model performance on the NYU cohort was summarized using the AUC and AP. AUC was used to evaluate discriminative performance across classification thresholds, whereas AP was used to summarize PR performance in the setting of class imbalance.

For analyses involving repeated CU population-wide control samplings, AUC and AP were calculated separately for each sampled evaluation set and summarized as the mean ± s.d. across samplings. The 95% confidence intervals were estimated from the s.e.m. across control settings. This evaluation strategy used multiple distinct population-wide control sets rather than upsampling preeclampsia cases, reducing the risk that performance estimates would be driven by repeated case duplication in the setting of low preeclampsia prevalence. By performing LOOCV across distinct population-wide control sets, these intervals summarize the variability of model performance across control definitions and provide an internal measure of robustness.

Threshold-based operating characteristics were also calculated, including TPR (sensitivity), FPR (1 − specificity), PPV, NPV and confusion-matrix counts. Threshold-based metrics were evaluated at prespecified classification thresholds and, where indicated, at thresholds corresponding to specified FPR operating points. Because PPV and NPV depend on outcome prevalence, predictive values from sampled CU evaluation sets were interpreted as sampled evaluation metrics and not as population-level screening estimates. Prevalence-adjusted PPV and NPV were calculated as described for clinical benchmarking.

Analyses involving subtypes, differential-diagnosis comparisons and demographic or clinical subgroups were interpreted descriptively when sample sizes were limited.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.