This prospective study was conducted in accordance with the declaration of Helsinki and received approval from the Ethics Committee at The First Affiliated Hospital of Shandong First Medical University & Shandong Provincial Qianfoshan Hospital (YXLL-KY-2022[043]).

Patients who were diagnosed with TC and had undergone thyroidectomy with newly discovered suspicious lymph nodes between January 2022 and July 2023 were considered for inclusion in this study. The inclusion criteria were as follows: (1) Postoperative patients with PTC who develop new abnormal cervical lymph nodes. (2) The lymph nodes had successfully undergone fine needle aspiration (FNA) and obtained clear cytological diagnosis results. The exclusion criterion was lymph nodes without clear cytological results because of inadequate cytological specimen or inconclusive result. All participants in this study provided informed consent prior to their participation.

Conventional grayscale ultrasound, SWE, and SMI images were obtained using an Aplio i900 or i800 system (Canon Medical Systems Corporation, Tochigi, Japan) equipped with a linear array transducer with a bandwidth of 5 to 18 MHz. Two radiologists (YH and FL) with experience of 3 and 8 years in ultrasound performed lymph node examinations. Two or more grayscale images were acquired for each lymph node. We evaluated the size, location, echogenicity of the cortex and visualization of the hilar fatty tissue.

After that, SMI and SWE were performed with minimum pressure from the transducer. The mean value of the vascular index (VI) and SWE for each lymph node was calculated from three measurements of each node obtained from longitudinal views. The parameters of SMI were shown as follows: velocity scale was 1.5 cm/s. For qualitative analysis, vascular pattern of lymph node was observed and reported. For quantitative analysis, VI was measured at the image with the strongest blood signal. VI is the ratio of number of pixels of vascular signal to the whole lymph node (tracing the outline of lymph node) and was automatically calculated and displayed on the screen (Fig. 1A). Margin of each lymph node was manually traced and the outline was on behalf of region of interest (ROI). We then switched to SWE mode and stiffness of lymph node was measured. Shear waves generated from the probe and displayed in a real-time color map. The average stiffness of lymph node (SWEmean) was assessed by placing ROI (tracing the outline of lymph node) and was expressed in kPa (Fig. 1B).

According to 2021 Korean Thyroid Imaging Reporting and Data System [17], we selected suspicious lymph nodes. Lymph nodes with any of the following features are defined as suspicious: cystic areas, hyperechoic foci (calcifications), cortical hyperechogenicity (focal/diffuse), or abnormal vascularization (peripheral/diffuse). We carefully evaluated the ultrasound features of the lymph nodes, including the presence or absence of a lymphatic hilum, the presence or absence of cortical hyperechogenicity, the presence or absence of calcification, the presence or absence of cystic components (Fig. 2) and the vascular pattern. There are four acknowledged vascular patterns, including hilar, peripheral, mixed and avascular or spot vascularity (Fig. 3). Mixed or peripheral vascularity was considered a malignant feature.

Due to the subjective nature of ultrasound feature evaluation, two ultrasound physicians with more than five years of work experience independently reviewed the images and evaluated the features. Consistency analysis was performed on the results. The final ultrasound features of the lymph nodes were determined based on the results of these two evaluations. In case of disagreement, a third senior physician was consulted for a final determination.

Under the guidance of color ultrasound, puncture site was selected. After routine disinfection, a sterile drape was placed and 2% lidocaine was used for local anesthesia. A 23G biopsy needle (BN-MAR-1/GA1; GMT Medical, Beijing, China) was inserted into the lymph node avoiding blood vessels. The aspiration was performed for three times and smear was made. After that, smears were sent for cytological examination. To prevent complications, all patients were monitored for 30 minutes after the procedure and underwent an ultrasound follow-up before discharge.

The gold standard of cytological diagnosis for distinguishing benign from metastatic lymph nodes were that thyroid follicular epithelial cells with papillary nuclear features were found in the cytological smears.

The normality of the distribution of the continuous data was verified using the Shapiro–Wilk tests. Patients were divided into two groups: benign and malignant. Continuous data (such as age, size, VI, and SWEmean) were expressed as mean ± standard deviation or median with range (interquartile range). The independent samples t-test or Mann–Whitney U test was used to compare the two groups. Percentages of categorical data were calculated and Fisher’s exact test or the χ² test were used to calculate the significance of the difference of categorical variables such as sex, laterality and location of lymph nodes, grayscale ultrasound features between the metastatic group and benign group. The reproducibility of the ultrasound feature evaluation was tested using Kappa consistency analysis. The diagnostic capacity of parameters was evaluated using receiver operating characteristic (ROC) curve analysis. The diagnostic accuracy of parameters was showed by calculating the area under the ROC curve (AUC). We calculated sensitivity, specificity, positive and negative predictive values and diagnostic accuracy. Using ultrasound features as the independent variable and the pathologic diagnosis as the dependent variable, linear regression analysis was performed to diagnose multi-collinearity of ultrasound indicators with statistical significance (p < 0.05) between the benign group and the metastatic group. Variables that could be included in the model were selected, and further Fisher stepwise discriminant analysis was used to eliminate variables that cannot be included. Some ultrasound features were then selected through Fisher’s stepwise discriminant analysis. The stepwise process involves iteratively adding and removing variables based on their Lambda values. In each step of the stepwise discriminant analysis, the variable with the smallest Lambda proceeds to the next step. A comprehensive diagnostic model for metastatic lymph nodes based on multimodal ultrasound was established. And self-testing method and leave one cross validation were performed for testing and ROC curve was constructed to evaluate the effectiveness of the diagnostic model.

All statistical analysis was performed using SPSS (version 25.0; IBM Corp., Armonk, NY, USA) and MedCalc version 15.8 statistical software (MedCalc Software bvba, Ostend, Belgium). A p value < 0.05 was considered to indicate statistical significance.

A total of 84 lymph nodes from 69 patients were eligible for inclusion in this study, including 51 females and 18 males with a mean age of 43.6 years (range, 20–74 yr). All patients included in our study were diagnosed with PTC. Table 1 displays the characteristics lymph nodes. There were no significant differences in the location, size of lymph nodes between the two groups. Among the 48 benign lymph nodes, 16 (33.33%) presented with a lymphatic hilum, while only 1 of the 36 malignant lymph nodes (2.78%) did. Cortical hyperechogenicity was obsevred in 8 of 48 benign lymph nodes (16.67%) and 24 of 36 (66.67%) metastatic lymph nodes. Calcification was detected in 3 of 48 benign lymph nodes (6.25%) and 8 of 36 (22.22%) metastatic lymph nodes. Cystic area was found in 4 of 36 (11.11%) metastatic lymph nodes. The hilar vascularity was present in 32 of the 48 benign lymph nodes (66.67%), while 28 of the 36 metastatic lymph nodes (77.78%) showed peripheral or mixed vascular signal. All these differences were statistically significant (p < 0.05). The VI and SWEmean values of metastatic lymph nodes were significantly higher than those of benign lymph nodes (p = 0.002, p < 0.001). The cutoff values for VI and SWEmean were 12.9 and 21.4 kPa, respectively. The kappa coefficients for sonographic features (the presence or absence of a lymphatic hilum, cortical hyperechogenicity, calcification, cystic area, and vascular pattern) between the two observers were approximately 0.688, 0.845, 0.950, 0.851, and 0.950, respectively. This analysis demonstrates that the ultrasound features were evaluated with good reproducibility and consistency.

Figure 4 illustrates the AUCs for various sonographic features. Among the gray-scale sonographic features, the vascular pattern exhibited the highest AUC. The sensitivity, specificity, positive predictive value, negative predictive value, and accuracy of vascular pattern in SMI were 77.78%, 85.42%, 80.00%, 83.67%, and 82.14% respectively. Table 2 lists the diagnostic performance of various sonographic features in distinguishing between benign and metastatic lymph nodes.

The pathological type of lymph nodes was used as the dependent variable, and the features of multimodal ultrasound with statistically significant differences (p < 0.05) were assigned as independent variables (Table 3). To avoid strong linear relationships between variables, which may lead to inaccurate parameter estimation in the Fisher discriminant equation, collinearity diagnosis was performed before constructing the diagnostic model. The values of tolerance of all the variables were > 0.1, and the values of variance inflation factor were all < 10 (Table 4). Therefore, it was believed that there was no multi-collinearity between the above variables, and Fisher discriminant analysis can be used for the study. Furthermore, Fisher’s stepwise discriminant analysis was performed for analysis. The stepwise process involved iteratively adding and removing variables based on their Lambda values. In each step of stepwise discriminant analysis, the variable with the smallest Lambda would proceed to the next step (Supplementary Table 1). The standardized canonical discriminant function coefficients for the four selected variables, namely lymphatic hilum (X1), cortical hyperechogenicity (X2), vascular pattern (X4), and SWEmean (X7), have been calculated to determine the magnitude of each variable’s influence on the dependent variable. The order of the variables by their influence, from the most to the least, was as follows: SWEmean (X7), cortical hyperechogenicity (X2), vascular pattern (X4), and lymphatic hilum (X1). Discrimination equations based on their classification function coefficients were established. Benign group: Y1 = −3.461 + 2.423X1 + 0.321X2 + 1.620X4 + 0.109X7. Malignant group: Y2 = −8.053 + 0.414X1 + 2.600X2 + 2.504X4 + 0.192X7. If Y1 < Y2, the lymph node would be diagnosed as metastatic lymph node.

The diagnostic model was validated using both self-validation and leave one cross-validation methods, and the results showed high accuracy, with accuracy rates of 83.3% and 79.8% (Supplementary Table 2, 3). The diagnostic performance of the model was evaluated by constructing an ROC curve, and the results showed that the area under the curve was 0.833, with a 95% CI of 0.736–0.906. The sensitivity (83.33%) and specificity (83.33%) were both good (Fig. 5).