This retrospective study was approved by the Ethics Committees of Changzhou First People’s Hospital and Suzhou Municipal Hospital. Clinical data, ultrasound, and CT images were collected from thyroid cancer patients treated at these hospitals between January 2022 and June 2024. Inclusion criteria: (1) pathologically confirmed primary classic PTC; (2) preoperative ultrasound and CT meeting analysis standards; (3) complete clinical data; (4) no prior thyroid surgery/ablation; (5) patients with PTC and/or concurrent benign thyroid conditions (nodular goiter, Hashimoto’s thyroiditis). Exclusion criteria: (1) non-classic PTC or other thyroid subtypes; (2) prior thyroid surgery/ablation; (3) history of head/neck cancer or familial cancer; (4) poor imaging quality; (5) incomplete clinical data; (6) non-curative surgery with persistent disease. A total of 799 PTC patients were enrolled, with data allocation as follows: Changzhou First People’s Hospital (training group, n=524; internal validation group, n=225) and Suzhou Municipal Hospital (external validation group, n=50).

Body mass index (BMI) was calculated as weight (kg)/height² (m²). Patients were classified as normal weight, overweight, or obese based on World Health Organization guidelines (8). Hashimoto’s thyroiditis was diagnosed based on elevated antibodies or ultrasound findings. Extrathyroidal extension (ETE) was defined as >25% tumor contact with the thyroid capsule (9). The largest lesion was used for pathological evaluation in multifocal cases. All diagnoses and lymph node statuses were confirmed pathologically. Preoperative laboratory tests included white blood cell count, platelet count, neutrophil count, lymphocyte count, monocyte count, and thyroid function, as well as inflammatory indices (lymphocyte-to-monocyte ratio (LMR), neutrophil-to-lymphocyte ratio, platelet-to-lymphocyte ratio, systemic immune-inflammation index).

All surgical procedures were performed by experienced thyroid surgeons following standardized protocols. Thyroidectomy procedures included: (1) Total thyroidectomy: complete removal of both thyroid lobes and isthmus; (2) Thyroid lobectomy: unilateral thyroid lobe removal with isthmus. Central lymph node dissection was performed in all patients, involving systematic removal of compartment VI lymph nodes. For patients with preoperative evidence of LLNM confirmed by fine-needle aspiration cytology, therapeutic lateral lymph node dissection was performed, involving systematic removal of levels II-V lymph nodes. A total of 102 patients underwent lateral neck dissection of levels II-V: 92 patients at Changzhou First People’s Hospital (61 training group + 26 internal validation group with confirmed LLNM, plus 5 patients with suspected but pathologically negative LLNM) and 10 patients at Suzhou Municipal Hospital with confirmed LLNM.

CT scans were performed with a Siemens Somatom Definition Flash dual-source CT scanner. Patients were positioned supine with slight neck hyperextension. The scanning range extended from the hyoid bone to the sternal manubrium, and if needed, to the aortic arch. Parameters included 120 kV, 200 mAs, 1.0 mm slice thickness, pitch 1.0, and a 200 mm × 200 mm field of view. Contrast-enhanced scans were performed using Iohexol (350 mg/ml iodine concentration), with dual-phase enhancement for arterial and venous phase imaging. Ultrasound was conducted using Philips iU22/EPIQ 5 or GE LOGIQ E9 systems. Experienced physicians obtained high-resolution tumor images and Doppler flow images, which were stored in DICOM format. Lymph nodes with suspicious features (e.g., round shape, absent echogenic hilum, microcalcifications) were classified as ultrasound-suspected LLNM. FNAC was then performed to confirm the histopathologic diagnosis of suspicious lateral lymph nodes. For patients with clinically suspicious LLNM confirmed by FNAC, thyroidectomy plus central neck dissection and therapeutic lateral neck dissection were performed.

To standardize image analysis and ensure cross-device reproducibility, all ultrasound and CT images were resampled to an isotropic voxel resolution of 1 mm³ using trilinear interpolation algorithms. For ROI segmentation, tumor regions of interest were manually delineated on contrast-enhanced arterial phase CT images by the two radiologists, as the arterial phase provides optimal tumor-to-background contrast for accurate boundary definition. Ultrasound images were min-max normalized to the (–1, 1) range, and tumor ROIs were manually delineated by two ultrasonographers using 3D-Slicer ( Supplementary Figure 1 ) following standardized protocols. A total of 874 quantitative features were extracted according to Image Biomarker Standardization Initiative (IBSI) guidelines, including morphological (25), first-order statistical (42), texture (729), and wavelet transform features (78) using 8 different wavelet decompositions (LLL, LLH, LHL, LHH, HLL, HLH, HHL, HHH). CT image analysis, independently performed by two blinded imaging physicians, involved spatial standardization using B-spline interpolation, Z-score normalization, and Gaussian filtering (σ=1.0 mm) for noise suppression. To mitigate potential inter-device variability across different scanner manufacturers (Philips, GE, Siemens), intensity harmonization using histogram matching was applied prior to feature extraction. LIFEx software extracted 1433 features following IBSI compliance standards, categorized into morphological (32), first-order statistical (75), texture (620), and filtered transform features (706). The filtered transform features included wavelet decompositions, Laplacian of Gaussian filters, and mathematical transformations (square, square root, logarithm, exponential). For reproducibility analysis of radiomics features, a random subset of 150 patients was selected for repeat feature extraction by the same two operators after a two-week interval to assess intra- and inter-observer reliability. Feature stability across different imaging platforms was assessed through intraclass correlation coefficient (ICC) analysis, with features demonstrating ICC >0.85 retained to ensure cross-device reproducibility.

All radiomics features were standardized using the zero-mean method, with highly correlated features (Spearman’s ρ>0.9) removed. Features with an intraclass correlation coefficient >0.85 were retained to ensure reproducibility.

T-tests or Mann-Whitney U tests were used to screen features with P<0.05 and |log2(Fold Change)|≥1. Subsequently, least absolute shrinkage and selection operator (LASSO) regression was applied to further reduce feature redundancy and optimize feature selection. The optimal λ value was determined using 10-fold cross-validation to identify key predictive features. Correlation heatmaps were generated to analyze relationships between radiomics features and clinical factors, ensuring they could provide complementary information for predicting LLNM.

To identify clinical features associated with LLNM, variance inflation factor (VIF) was calculated through collinearity analysis to exclude variables with multicollinearity. Logistic regression analysis with stepwise regression was then performed on the remaining variables to screen for independent risk factors.

Four machine learning models were constructed: Random Forest (RF), Gradient Boosting Machine (GBM), Support Vector Machine (SVM), and K-Nearest Neighbors (KNN). Models were trained using 10-fold cross-validation to avoid overfitting. While advanced deep learning approaches, including semi-supervised learning frameworks (10) and divide-and-conquer architectures (11), show promise in medical imaging, we selected traditional machine learning algorithms for better interpretability and performance with our dataset size. Hyperparameters were optimized through grid search. Model performance was evaluated by area under the curve (AUC), sensitivity, specificity, accuracy, and related metrics. DeLong tests compared model differences, while clinical utility was assessed using decision curves. SHapley Additive exPlanations (SHAP) and Local Interpretable Model-agnostic Explanations (LIME) were used to interpret model predictions, with SHAP providing global feature importance and LIME offering local explanations for individual cases.

Statistical analyses were performed using R (Version 3.5.3), SPSS (Version 25.0), and Python (Version 3.12.0). Categorical variables were compared using chi-square or Fisher’s exact tests. Continuous variables were compared using t-tests for normally distributed data and Mann-Whitney U tests for non-normally distributed data. Model evaluation included receiver operating characteristic curves, DeLong tests, Decision curve analysis (DCA), and SHAP and LIME analyses. P<0.05 was considered statistically significant.

In the training group, 146 patients (27.9%) were male and 378 (72.1%) were female, with a mean age of 44.0 ± 12.1 years. The internal validation group consisted of 51 males (22.7%) and 174 females (77.3%), with a mean age of 42.7 ± 11.6 years. The external validation group included 10 males (20.0%) and 40 females (80.0%), with a mean age of 44.2 ± 12.3 years. The incidence of LLNM was 11.6% in the training group, 11.6% in the internal validation group, and 20.0% in the external validation group ( Table 1 ). No statistically significant differences in clinical or pathological characteristics were observed between the training and internal validation groups (all P>0.05). Follow-up surveillance revealed that 7 patients developed contralateral residual thyroid recurrence. No patients who did not undergo initial lateral neck dissection developed subsequent lateral regional recurrence

To identify independent risk factors for LLNM, we first performed collinearity diagnostics. Variables demonstrating significant multicollinearity (VIF>10) included platelet count (VIF=99.79), neutrophil count (VIF=89.25), lymphocyte count (VIF=29.22), neutrophil-to-lymphocyte ratio (VIF=178.73), platelet-to-lymphocyte ratio (VIF=160.64), and systemic immune inflammation index (VIF=96.32). These variables were excluded from further analysis to enhance model stability.

Multivariate logistic regression analysis of the remaining variables identified five independent risk factors for predicting LLNM ( Table 2 ). Maximum tumor diameter was significantly associated with increased risk, with progressively higher odds for larger tumors (>1 to ≤2 cm: OR=2.494, 95% CI: 1.212–5.132, P=0.013; >2 to ≤4 cm: OR=7.851, 95% CI: 3.072–20.066, P<0.001; >4 cm: OR=13.032, 95% CI: 3.253–52.212, P<0.001). Multiple lesions (≥2 lesions: OR=2.846, 95% CI: 1.436–5.639, P=0.003) and tumor location in the upper pole (OR=5.181, 95% CI: 2.550–10.524, P<0.001) were also identified as independent risk factors. Additionally, decreased monocyte count (OR=0.004, 95% CI: 0.000–0.070, P<0.001) and a lower LMR (OR=0.524, 95% CI: 0.410–0.671, P<0.001) were significantly associated with LLNM.

Differential feature analysis using independent sample t-tests (for normally distributed data) or Mann-Whitney U tests (for skewed distributions) was performed to compare radiomics features between groups with and without LLNM. Features with P<0.05 and |log2(Fold Change)|≥1 were selected. This initial screening reduced the 874 ultrasound radiomics features to 100 and the 1433 CT radiomics features to 97 ( Figures 1a ).

LASSO regression was subsequently applied to further eliminate redundant features while balancing model simplicity and predictive performance. Using 10-fold cross-validation to determine the optimal λ value (corresponding to minimum binomial deviance), we identified 4 key ultrasound features (1 first-order feature and 3 texture features) and 11 key CT features (1 morphological feature and 10 texture features) ( Figures 1b ).

Correlation analysis using heatmaps evaluated the relationships between radiomics features and clinical factors (tumor diameter, number of lesions, tumor location, monocyte count, LMR and LLNM). Results demonstrated weak correlations between clinical factors, except for LLNM, and radiomics features ( Figures 2a ), indicating that these features could provide complementary information for constructing more robust prediction models.

Based on the 5 clinical features, 4 ultrasound radiomics features, and 11 CT radiomics features, we constructed four machine learning models: RF, GBM, SVM, and KNN. Comprehensive evaluation revealed that the GBM model demonstrated superior overall performance, achieving AUCs of 0.973, 0.803, and 0.975 in the training, internal validation, and external validation sets, respectively ( Table 3 ).

In the training set, the GBM model achieved an accuracy of 0.914, with 0.894 specificity and 0.957 sensitivity. In the internal validation set, these metrics were 0.725, 0.889, and 0.794, respectively. In the external validation set, the GBM model demonstrated excellent performance with an accuracy of 0.900, specificity of 0.950, and sensitivity of 0.800 ( Table 3 ). Although the RF model showed a comparable AUC in the external validation set (0.955 vs. 0.975), the GBM model exhibited better overall diagnostic metrics and F1 score (0.842 vs. 0.870). DeLong tests revealed statistically significant differences between the GBM model’s AUC and other models (P<0.05) ( Figures 3a ). Combined with the higher AUC values, this indicates that the GBM model achieved superior discriminative performance compared to the RF, SVM, and KNN models. DCA further validated these findings: across clinically relevant threshold ranges, the GBM model (yellow line) maintained the highest net benefit, followed by KNN, SVM, and RF models ( Figure 3 ).

The observation-prediction probability scatter plots revealed the predictive characteristics of each model ( Figures 4a ). The RF model ( Figure 4 ) exhibits a characteristic “striped” pattern due to its ensemble voting mechanism, with good overall separation but notable uncertainty in the mid-probability range. The GBM model ( Figure 4 ) demonstrates a more continuous probability distribution with clearer class separation, reflecting its superior calibration and generalization capability. The SVM model ( Figure 4 ) displays a polarized prediction pattern, clustering probabilities at extremes, which corresponds to its lower sensitivity (0.404) in internal validation due to misclassifications, particularly in the external validation set. The KNN model ( Figure 4 ) produces a “stepped” probability pattern, indicating limited discrimination ability and aligning with its lower accuracy (0.567) and F1 score (0.519) in external validation. These probability distributions visually corroborate the performance metrics in Table 3 , further supporting the GBM model as the most reliable classifier for LLNM prediction.

SHAP analysis was employed to interpret the GBM model’s prediction process. SHAP analysis ( Figure 5 ) identifies LMR as the most influential feature, followed by wavelet-based texture features (wavelet-LHH_glszm_GrayLevelNonUniformity, wavelet-LLL_glszm_ZoneVariance), morphological characteristics (wavelet-HHH_glszm_SizeZoneNonUniformity), and exponential transformations of gray-level matrices. The color gradient highlights how feature values impact prediction outcomes. Figure 5 visualizes the model’s decision path, mapping cumulative feature contributions from the base value (0.4) to final probabilities (0–1), revealing complex feature interactions with stronger positive contributions toward LLNM prediction.

LIME analysis ( Figure 5 ) demonstrates clear probability stratification, with distinct low- and high-risk regions and a minimal overlap between positive and negative cases, supporting robust risk discrimination. These findings indicate that while LMR is the dominant predictor, the GBM model’s superior performance arises from the integration of clinical and radiomic features, enhancing LLNM prediction beyond conventional visual assessment.

To facilitate clinical translation of our multimodal prediction model, we developed an interactive web-based calculator interface ( Figure 6 ). This platform integrates all identified predictive features into a streamlined clinical workflow, including: (a) dropdown menus for clinical features (tumor diameter, lesion number, tumor location); (b) input fields for laboratory parameters (monocyte count, LMR); (c) DICOM image upload functionality for ultrasound and CT images; and (d) real-time LLNM probability calculation. This implementation provides healthcare providers with a practical tool for preoperative risk stratification in PTC patients.