The ethics committee of the First Affiliated Hospital of Chongqing Medical University approved this study (No:2022-63), and the demand for informed consent was exempted. Primary PDAC patients underwent surgery resection with standard regional lymph node dissection from January 2015 to September 2021 were collected. The enrolled criteria were as follows: (1) patients underwent arterial phase (AP) and venous phase (VP) scanning within 2 weeks prior to surgery; (2) pancreatic tumors were observed on CECT images; (3) the diagnosis of PDAC and LNM were confirmed pathologically; and (4) clinical data were completed. The exclusion criteria were as follows: (1) tumor diameter was at least 1.0 cm; (2) images with evident noise or severe motion artifacts; (3) treatment or biopsy before imaging scanning; and (4) other primary tumors existed. The patients were partitioned into a primary group (training group and validation group) and a test group at a proportion of 8:2 using random sampling. The patient selection flowchart is described in Figure 1 .

One 128-slice multidetector-row CT scanner (SOMATOM Definition Flash, Siemens Healthineers) was performed. The CT scanning parameters were set as follows: 120 kV; 300 mA; 0.7 pitch; collimation, 128 × 0.6 mm; beam collimation, 160×0.5 mm; matrix, 512×512; and gantry rotation time, 0.5 s. Non-ionic contrast agent (Ultravist 370, Bayer Schering Pharma) was injected into the antecubital vein using a pump injector (Medrad Mark V plus, Bayer). The injection dose was 1.2 ml/kg and the flow rate was 3.5 ml/s. Then, normal saline of 40 ml was injected to flush the tube. Unenhanced phase was scanned first. Approximately 15 s after the abdominal aorta reaching 100 HU, AP scanning was performed, and VP scanning was performed 30 s after the finish of the AP scanning.

Patient data were acquired from the electronic medical records. A total of 20 clinical, pathological, images and laboratory characteristics were evaluated referring to the World Health Organization and the American Joint Committee on Cancer (AJCC) TNM Staging System Manual, 8th Edition (24). Image characteristics were assessed by two radiologists with 8 and 10 years’ clinical experience, respectively. A consensus was reached when difference in opinion existed. The characteristics were classified as follows: (1) clinical characteristics: gender, age, abdominal pain, backache, pancreatitis, jaundice, operation method; (2) pathological characteristics: histological grade, duodenal invasion, surgical margin status, perineural invasion; (3) image characteristics: CT-reported tumor size, tumor location, clinical T stage, parenchymal atrophy, pancreatic duct dilatation, and common bile duct dilatation; (4) laboratory characteristics: carcino-embryonic antigen (CEA) level, carbohydrate antigen 19-9 (CA19-9) level, and total bilirubin (TBIL) level. Specific characteristics description can be referred in the Supplementary Material. Meanwhile, valuable characteristics were selected from above-mentioned characteristics using univariate and multivariable logistic regression analysis for clinical model building, except for that, we chose some feature, which is deemed meaningful in clinical experience for predicting tumor heterogeneity to further perfect this clinical model.

In this study, traditional manual segmentation layer by layer is needless. We designed one 3D stereochemical box as the VOI. First, the tumor was localized by two experienced radiologist without precisely segmentation. One radiologist localized the lesion and another radiologist checked the accuracy of location. Second, this 3D box, which contained the complete tumor and slight peritumoral tissue along height-axis, width-axis, and depth-axis in every slice, is determined according to the location mark. Specifically, the original images could be recognized as a 3D array, and the largest peripheral box including the tumor is also a 3D array, and the index corresponding to each pixel of the tumor 3D box in the original images could be determined via computer language. The corresponding original images were cropped out according to the index for retaining tumor 3D box, and the specific demand of box size only need to retain whole tumor as much as possible because deep learning algorithm is more robust and strong than traditional machine learning method, the requirement of precise segmentation in deep learning is not as strict as traditional radiomics. This method is consistent with An’s study about the localization and segmentation technology (26). Third, after subtracting the periphery regions out of tumor 3D box, original 3D images and 3D segmentation masks (ground truth) were resampled to a specified resolution of 5*224*224, and images were processed to (0,1) using min-max normalization. A complete workflow is displayed in Figure 2 .

Resnet as a celebrated CNN is proposed by Microsoft research. In 2015, Resnet won the image classification and object recognition competition in ImageNet Large Scale Visual Recognition Challenge (ILSVRC). The Resnet 18 represent 18 layers deep networks constructed by internal residual blocks, which were implemented by shortcut connection, thereby performing identity mapping and solving degradation problem caused by over-deeper layers. But Resnet 18, as one 2D CNN model, could only extract feature in single slice of 2D CT images without extracting the entirety in real 3D structure. Thus, Resnet 3D-18 was optimizing and upgrading on the basis of Resnet 18. Resnet 3D-18 could extract context features comprehensively and globally using automatic parameter learning, thereby avoiding losing stereoscopic information. Currently, Resnet 3D-18 is highly appreciated by its excellent learning performance and optimization ability, and increasingly applied in image segmentation, recognition, classification. Moreover, it has been reported that Resnet-3D model achieve better accuracy compared to 2D ones, and the deeper networks (34 layers) show little gain over 18 layer ones (33, 34), so the Resnet 3D-18 architecture was appropriate and selected for subsequent analysis.

Pretrained Resnet 3D-18 model on the ImageNet database with the method of self-supervised learning was implemented for feature extraction (35, 36). In previous studies of deep learning neural networks (25, 27), end-to-end learning methods were regularly used without feature selection procedure; however, these features are not identically significant to the issues. In this experiment, multiple attempts in feature selection method were conducted, including decision tree (DT) and least absolute shrinkage and selection operator (LASSO) model. Ultimately, the LASSO model was selected, which is consistent with prior deep learning radiomics studies (37). All extracted features are fed into the LASSO model. The loss function of LASSO is calculated as below:

Where N is the number of all samples, y represents the true label, and X represents feature, w represents weight of feature, and hyperparameter λ denotes penalty coefficient. Those features with non-zero weight were retained. The λ was searched in 0.005-0.02 using the traversal method. The larger the λ is, the stronger the effect of the regular term is, and more weights would be compressed to zero. The model error was calculated by putting the possible values into the model. The consequence of λ yielding the minimum error was selected as the optimal hyperparameter.

In this experiment, we have attempted six classifiers including logistic regression (LR), k-nearest neighbor (KNN), SVM, DT, random forest (RF), and fully connected neural network (FCN). Finally, the combination of LASSO+FCN achieving the best robustness and stability was selected, followed by the combination of DT + RF, which were easily influenced by randomness of the data set, thereby being excluded. Then, the features retained by LASSO were fed into FCN, which is a linear model to calculate the final probability score. The architecture of FCN was constructed by one hidden layer and one output layer with training epochs of 30. Limited Broyden-Fletcher-Goldfarb-Shanno (Lbfgs) optimizer was adopted to minimize the loss function—binary cross-entropy—which is calculated as follows:

Where y _i represents true label of sample i, the positive sample is 1, and the negative sample is 0. p _i   represents the probability of the sample predicted as a positive class. In hidden layer, neurons were set to 1,000, and the activation function used was rectified linear unit (Relu), which could make the output of some neurons be zero, thereby contributing to the sparsity of the network, decreasing the interdependence of parameters, and relieving the occurrence of over-fitting issues. In output layer, neurons were set to two, and the activation function was softmax function via mapping the values of the output layer to the 0-1 interval as a probability distribution where z represent values of the output layer. The larger the softmax value is, the better the model predicts. The formula of softmax is:

The CNN training process included two aspects (backward propagation and forward propagation). First, five times 5-fold cross-validation was set for training model. The primary group was stochastically shuffled and divided into 5-fold averagely. In every time among five times, 1-fold was regarded as validation group successively, and the other 4 folds were regarded as training group in this order. Then, the images were fed into the network, and predictive result through forward propagation displayed in the network’s output layer. Simultaneously, the model parameters were updated and decided via backward propagation until achieving the minimum difference between the true label and predictive result. In addition, we used early stopping, which is a technique use to stop training before overfitting occurs. The use of early stopping can obtain the best generalization performance and prevent overfitting by intercepting the model with the best results in the whole process of model training. In total, three 3D CNN models including AP DLR model, VP DLR model, and AP+VP DLR model were separately built. The predictive ability of three CNN models was exhibited among every time in 5-fold cross-validation. The general prediction result was fused among these 5-fold groups by averaging the scores. The independent test group was employed to evaluate the model performance. It should be noted that no patients in the test group directly or indirectly participate in the model training process, thereby avoiding data leakage.

In our study, the following five quantitative indicators were calculated: the area under the curve (AUC), accuracy, precision, recall, and F1 score. Due to the imbalances between the LNM group and non-LNM group, we used the AUC as our principle evaluation indicator, followed by accuracy. Furthermore, the Radiomics score (Rad score) were calculated through a linear combination of features’ weighted coefficients. To make a comparison of diagnostic performance between the 3D CNN models with physician-level accuracy, CECT images of all 139 patients were respectively reviewed by two abdominal radiologists (senior radiologist and junior radiologist) following double-blind principle. The site of lymph node resection is consistent with that of images measurement. The accuracy, sensitivity, specificity, positive predictive value (PPV) and negative predictive value (NPV) were evaluated in performance of radiologists.

Continuous variables were analyzed using Student’s t test or the Mann–Whitney U test; Categorical variables were analyzed by chi-square test or Fisher’s exact test. Wilcoxon rank sum test was used to compare Rad scores in negative and positive groups. Delong test was applied to evaluate the discrimination ability among AUCs, and P<0.05 was regarded as statistically significant. Multivariable logistic regression analyses use the likelihood ratio test with Akaike’s information criterion as the stopping rule. All statistical tests were executed with SPSS software (version 25.0), R software (version 4.0.5), and Python software (version 3.8.0). In this study, the pre-processed and feature extraction approaches were conducted using SimpleITK, numpy, and scikit-learn package; the PyTorch 1.0 configuration was arranged to build the neural network. The training process was conducted on Ubuntu OS with an Intel Xeon E5 2687W V3, NVIDIA GeForce 1080ti GPU, and 16 × 8GB of RAM.

A total of 139 patients (99 men, 40 women) were recruited for our study, and were partitioned into a primary group (n = 111) and a test group (n = 28). The LNM rates in every group were 35.1% (39/111) and 32.1% (9/28), respectively, and P>0.05 ensuring grouping consistency. In the primary group, histological grade, duodenum invasion, and CT-reported tumor size were notably different between non-LNM and LNM cohorts ( Table 1 ). There were no significant differences in the test group. Some studies suggested that CA19-9 level is an independent predictive factor for LNM, and the higher CA19-9 levels indicating a worse patient’s condition (38). Thus, CA19-9 is enrolled for subsequent analysis. Univariate and multivariable logistic regression analysis demonstrated that histological grade, CT-reported tumor size, and CA19-9 were independent predictors of LNM in PDAC ( Table 2 ). PDAC patients with LNM were more likely to have higher histological grade (OR,0.175; 95% CI, 0.061 to 0.499), larger CT-reported tumor size (OR, 1.182; 95% CI, 1.069 to 1.307), abnormal CA 199 level (OR, 4.139; 95% CI, 1.300 to 13.175).

In this study, we attempted experiments comparing DLR model with feature selection procedure or without feature selection procedure. Without feature selection, the learning curves were overfitting in the training group, and were not converged in validation group and test group in all three models. However, learning curves improved significantly and reached perfect fitting with feature selection. The training curves of three models could be found in Supplementary Figure S1 . This finding suggested that feature selection may be an important method in deep learning researches, especially in small data sets. Adding feature selection procedure might achieve better performance.

The Resnet 3D-18 model separately extracted 512 deep learning features from AP and VP. Then, the most representative and significative characteristics were retained by LASSO model with λ setting as 0.01 ( Figure 3 ). A total of 45, 49, and 59 features were selected from AP model, VP model, and fusion of AP+VP model, respectively. Indeed, more features were retained due to a large base of total deep learning features, which is greatly increased by the hidden layer structure of deep neural networks. These features with non-zero weights were ultimately assigned to construct DLR model. The feature heatmap and Rad score were plotted in Figure 4 , and specific Rad score computational formula, including corresponding coefficients, was in Supplementary Material. No matter in any of the three 3D CNN models, significant difference was observed in the Rad score between patients in LNM and non-LNM groups (all P< 0.01).

The performance of different 3D CNN models are showed in Table 3 and Figure 5 . Overall speaking, the AP + VP DLR model reached the optimal ability for identifying LNM in PDAC with an AUC of 0.995 (95% CI: 0.989-1.00) and an accuracy of 0.969 in the training group; an AUC of 0.940 (95% CI:0.910-0.971) and an accuracy of 0.883 in the validation group; an AUC of 0.949 (95% CI:0.914-0.984) and an accuracy of 0.836 in the test group, followed by AP DLR model with an AUC of 0.962 (95% CI: 0.951-0.972) and an accuracy of 0.926 in the training group; an AUC of 0.884 (95% CI: 0.800-0.968) and an accuracy of 0.821 in the validation group; an AUC of 0.872 (95% CI: 0.823-0.921) and an accuracy of 0.736 in the test group; and the VP DLR model reached an AUC of 0.967 (95% CI:0.955-0.979) and an accuracy of 0.903 in the training group; an AUC of 0.884 (95% CI:0.829-0.938) and an accuracy of 0.784 in the validation group; an AUC of 0.844 (95% CI:0.820-0.867) and an accuracy of 0.764 in the test group. The AP model and VP model achieved similar performance, while both were not as good as AP+VP model.

Moreover, we could observed that all models reached decent outcome with most AUCs ranging from 0.92-1.00 in every fold among 5-fold cross-validation. The ROC curves including every fold curve and average curve in training group and test group in all models kept stable and slight wobble. The performance comparison of the different 3D CNN models in training, validation, and test group were displayed in Figure 6 . No matter in training, validation, and test groups, all statistical indicators in AP+VP model were the highest than that of other models, except in test group, precision is lower in AP+VP model than AP model. Delong test demonstrated no significant difference was observed in test group with AUCs ranging from 0.844 to 0.949 (all P>0.05), denoting the robustness and consistency of AP, VP, and AP+VP models.

The performance of radiologists and clinical model were showed in Table 4 and Figure 7 . In primary group, the accuracy of senior radiologist was 0.838, which was slightly lower than that of AP DLR model and VP DLR model with an accuracy of 0.874 and 0.844, respectively. In test group, the accuracy of senior radiologist was 0.821, which is also higher than that of AP DLR model and VP DLR model with an accuracy of 0.736 and 0.764, respectively. However, the predictive ability of senior radiologist was lower than that of AP+VP DLR model in both primary group and test group. The predictive ability of junior radiologist was lower than that of AP DLR model, VP DLR model, and AP+VP DLR model in primary group and test group. A clinical model was constructed by concatenating histological grade, CT-reported tumor size, and CA19-9 level. The clinical model achieved an AUC of 0.747 (95% CI:0.657-0.837) in primary group and an AUC of 0.737 (95% CI:0.549-0.925) in test group ( Figure 6 ), which is lower than all 3D CNN models and radiologists.