This retrospective study was approved by the institutional review board of center 1 (Henan Province People’s Hospital, HPPH) and center 2 (Henan Cancer Hospital, HCH), and the requirement for informed consent was waived.

We searched hospital information system with cervical cancer between January 2013 and February 2022, and found 380 consecutive cervical cancer patients who accepted a complete NACT regimen and radical hysterectomy with definite pathological results in two hospitals. According to postoperative pathological reports, patients were labeled as LN positive [LN (+)] and negative [LN (-)] groups. The inclusion criteria were as follows: (i) patients with FIGO stage IB3-IVA; (ii) patients who underwent pelvic MRI for pretreatment evaluation before NACT; (iii) Patients with a b-value of 800 s/mm² on DWI; (iv) patients who received 1 to 3 cycles of NACT followed by hysterectomy and lymphadenectomy; (v) patients who had definite postoperative pathological results. The exclusion criteria were as follows: (i) patients with poor MRI quality; (ii)patients who were treated with NACT and brachytherapy or external beam radiotherapy simultaneously; (iii) patients who were judged to be inoperable after NACT and received CCRT; (iv) patients who underwent selective lymphadenectomy without hysterectomy after NACT.

Clinicopathological characteristics of patients, such as age, FIGO stage, NACT regimen, and tumor type were retrieved from their medical records. The LN status was used to adjust FIGO stage of the patients to conform with the 2018 staging system. Patients with LN (+) in MR images were classified as IIIC no matter what the original FIGO stage was (10). The LN was considered LN (+) if it met one of the following criteria: (a) a shortest diameter of more than 8mm; (b) ring enhancement showing necrosis inside the LN (23). The patient recruitment flowchart is shown in Figure 1 .

The NACT protocol included 1 to 3 cycles of platinum-based intravenous chemotherapy administered at three-week intervals. Gynecologists evaluated the feasibility of surgical resection and performed radical hysterectomy and lymphadenectomy three weeks after completion of NACT. If pelvic node disease was detected during intraoperative examination or if bulky aortic nodes were identified on preoperative imaging or during surgery, aortic lymphadenectomy was performed.

All the resected surgical samples including the uterus and LNs were formalin fixed, paraffin embedded, and then histopathologically diagnosed by the pathologists of the two hospitals. LN status was extracted from the pathological reports.

Baseline MRI data (within 1 week before the NACT regimen started) were collected from two 3.0-T MR Platforms from Siemens Healthineers (MAGNETOM TrioTim 3.0 T, Skyra 3.0 T) and GE healthcare (discovery MR 750 3.0 T, signa HDxt 1.5 T). The MRI protocols of each hospital included sagittal T2-weighted imaging (Sag-T2WI) sequence, axial diffusion-weighted imaging (Ax-DWI) sequence, and the delayed phase of dynamic-contrasted-enhanced sagittal T1-weighted imaging (Sag-T1C) sequence. The comprehensive parameters of MR image acquisition of each hospital are detailed in Supplementary Table S1 .

All images were analyzed independently by two radiologists, who were blinded to LNM results.

The pretreatment MRI sequences were collated for tumor segmentation and feature generation. The regions of interest (ROI) were manually delineated on the cervical tumor slice by slice using ITK-SNAP (version 3.8.0 software: www.itksnap.org) containing the tumor mass on Sag-T2WI, Ax-DWI and Sag-T1C. To ensure accurate tumor boundaries, ROIs on all slices were carefully delineated manually by a radiologist with 13 years of experience in abdomen and pelvis imaging, who was blinded to the pathological findings.

Then, the ROI-based features of images from 30 randomly selected patients were re-extracted to test the features’ reproducibility for extraction in a blind fashion by the same radiologist. Only features with ICCs >0.75 were considered stable and selected for subsequent analysis.

The process of feature extraction and selection was carried out using the in-house software (PyRadiomics package version 2.3.0). First, all images were normalized using the max–abs scaling algorithm. Second, resampling and interpolation were also applied. The Sag-T1C images were reconstructed to 1mm×1mm×1mm voxel size, the Ax-DWI and Sag-T2WI images were reconstructed into 1mm×1mm×4mm image voxel sizes using a neural network-based B-spline surface construction approach to convert 2D ROI to 3D volume of interest (VOI). A total of 2264 features were extracted from each labeled ROI, including first order statistics features, shape-based (2D and 3D) features, gray level co-occurrence matrix (GLCM) features, gray level dependence matrix (GLDM) features, neighborhood gray tone difference matrix (NGTDM) features, gray level run-length matrix (GLRLM) features, and gray level size zone matrix (GLSZM) features.

After 3 months, 40 patients [30 LN (-) and 10 LN (+)] were randomly selected and resegmented by the same radiologist (QW) to assess intra-observer reliability. Only features with intra-observer ICCs > 0.75 were considered stable and selected for subsequent analysis.

The z-score standardization method was used for feature normalization. The K-Best method was first used to remove features with low correlation with classification labels in order to reduce computational complexity and prevent overfitting. And then LASSO was used to remove redundant features. Finally, the radiomics features significantly correlated with LNM of LACC were selected.

The Sag-T2WI model, DWI model, and Sag-T1C model were constructed based on the selected radiomics features from Sag-T2WI, Ax-DWI, and Sag-T1C sequence, respectively. In order to build a more individualized predictive model, the correlation analysis was performed among all features of the three sequences, and the highly correlated features (r >0.75) were excluded, then the selected features were fused with clinical characteristics to construct a combined model. All of the models were developed using random forest classifier based on the training dataset, and the performance of which were evaluated and compared in the validation dataset.

The development and validation of models were performed with InferScholar platform version 3.5 (InferVision, Beijing, China). The study framework is shown in Figure 2 .

Statistical analyses were performed using IBM SPSS software (version 28.0, Chicago, IL, USA). The differences between continuous clinical variables were evaluated by the Mann Whitney U test. The spearman correlation analysis was used to measure the correlation between arbitrary features. The visualization of the features used to build the combined model is illustrated using heatmap drawn by the R language with “ComplexHeatmap” package. The ROC and DCA were plotted using the R language with the “rms” package (version 4.2.1) and the “rmda” package (version 1.6), respectively. Delong’s test was used to compare the differences among AUCs of different models and different MRI instruments. The Hosmer-Lemeshow test was used to evaluate the consistency between the predicted and actual positive rate of LN in the validation set. A two-sided P <0.05 was considered statistically significant.

The study included 282 patients diagnosed with LACC. The training set comprised 171 patients from HPPH, consisting of 43 LN (+) and 128 LN (-) cases. Simultaneously, the validation set involved 111 patients from HCH, with 28 LN (+) and 83 LN (-) cases. Table 1 presents the clinical characteristics of the patients. Importantly, only the FIGO stage showed statistical significance between the LN (+) and LN (-) groups in both datasets (P < 0.001).

A total of 5 Ax-DWI features (5 textural features), 4 Sag-T1C features (4 textural features) and 4 Sag-T2WI features (1 first order feature and 3 textural features) were selected ( Table 2 ). And these features were utilized to construct three single-sequence models, respectively. Spearman correlation analysis was used to eliminate redundant features in the following way: when a pair of features showed a coefficient >0.75, the feature with the lower P value was remained. As a result, the Feature 13, wavelet_glcm_wavelet-LHL-Autocorrelation, was excluded. Subsequently, 4 Ax-DWI features, 4 Sag-T1C textural features, 1 Sag-T2WI first order feature, and 3 Sag-T2WI textural features were employed to construct the combined model ( Figure 3 ). According to Mann Whitney U test, only FIGO stage (P <0.001) used for the construction of the combined model.

Leveraging the remaining features, a combined model was developed, incorporating the FIGO stage. The ROC curves illustrating the performance are shown in Figure 4 . The AUCs value, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) of each model are summarized in Table 3 .

Notably, among the single-sequence models, the Ax-DWI model demonstrated the best predictive performance with AUCs of 0.855 (95%CI: 0.791-0.919) in the training set and 0.753 (95%CI: 0.638-0.867) in the validation set. There were no significant differences among AUCs of different single-sequence models in both sets (training set, Ax-DWI vs Sag-T2WI P =0.238, Ax-DWI vs Sag-T1C P=0.243, Sag-T2WI vs Sag-T1C P=0.921; validation set, Ax-DWI vs Sag-T2WI P =0.236, Ax-DWI vs Sag-T1C P=0.200, Sag-T2WI vs Sag-T1C P=0.838).

Furthermore, the diagnostic efficacy of the combined model, integrating radiomics features and the FIGO stage for LNM in LACC patients, surpassed that of among other models with AUCs of 0.889 (95%CI: 0.833-0.945) in the training set and 0.859 (95%CI: 0.781-0.936) in the validation set. However, according to the Delong test, only the difference between the combined model and the DWI model was not statistically significant in both the training (P =0.268) and validation sets (P =0.059). As shown in Figure 4 , the model based on FIGO stage exhibited superior performance in the validation set with AUC of 0.806 (95%CI: 0.742-0.870) and displayed significant differences from the Sag-T2WI and Sag-T1C models (P <0.05).

Detailed comparisons between different models are presented in Supplementary Table S2 . Additionally, no significant differences were observed among different MR image instruments (all P >0.05) ( Supplementary Table S3 ).

The DCA of the prediction model in the validation set is depicted in Figure 5 . If the threshold probability is below 80%, the combination model exhibited a higher net benefit compare to the FIGO stage and other single-sequence models in distinguishing LN (+) from LN (-).

The Hosmer-Lemeshow test indicated that in the validation set, all models showed good agreement between the predicted LNM probability and the observed incidence (Sag-T1C model, P =0.904; Ax-DWI model, P =0.078, Sag-T2WI model, P =0.127; Combined model, P =0.720).