CRC patients with lung metastases who underwent percutaneous RFA under CT guidance between August 2016 and January 2019 were enrolled in this study. Patients were recruited based on the following eligibility criteria: (1) histologically confirmed CRC; (2) ablated lung metastases with maximum diameter ≤3 cm; (3) chest enhanced CT examination within 4-6 weeks before RFA; (4) complete CT images during the procedure; (5) re-examination by chest enhanced CT at least 6 months after RFA; (6) technically successful ablation; (7) adequate normal organ function. Exclusion criteria, based on the European Society for Medical Oncology (ESMO) guidelines (38) were: (1) > 5 lung metastases; (2) maximum diameter > 3 cm; (3) other local or regional treatments such as radiotherapy before or after RFA; (4) incomplete clinical data; (5) second ablation (i.e., re-ablation). We allowed the inclusion of patients with multiple nodules and analyzed each nodule individually. A cohort of 198 patients with 479 lung metastases who received RFA was retrospectively selected ( Figure 1 ).

The following clinical data were retrieved: age at diagnosis, gender, serum tumor markers, including carcinoembryonic antigen (CEA) and cancer antigen 19-9 (CA19-9). The radiological data were recorded as follows: the location of pulmonary metastases, proximity to the heart, great blood vessels (diameter > 3 mm), pleura or diaphragm (within 1 cm) through the preoperative CT images; the IAH or pneumothorax were acquired.

All CT examinations (United Imaging uCT 760, Shanghai United Imaging Medical Technology Inc., China and Philips Brilliance 64 slice, Philips Medical Systems Inc., USA) were performed with a fixed tube current of 200 mA and a tube voltage of 120 kVp. The pixel spacing ranged from 0.684 to 0.748 mm, and the slice thickness was 1 mm or 1.5 mm. The intraoperative CT images were of fixed tube current of 200 mA, tube voltage of 120 kVp, and slice thickness of 1 mm or 3 mm. The image reconstruction method of both CT scanners is iterative reconstruction. Radiological follow-up consisted of chest-enhanced CT scans performed at 1, 3, 6, 12 months, and every 6 months after that. The shortest follow-up time was over 6 months.

This study was approved by the Institutional Review Board of the Ethics Committee of Fudan University Shanghai Cancer Center. Written informed consent was obtained from all patients.

RFA mainly utilizes 460 ~ 480 kHz high-frequency current to heat a tissue volume around a needle electrode and induce focal coagulative necrosis with minimal injury to surrounding tissues (39, 40). Here, RFA was performed using a radiofrequency applicator (MedSphere International), with the mode of temperature control or impedance control for choice. The power settings were adjusted according to the manufacturer’s protocols: 5 min for a 2.0 ~ 2.5 cm active tip at 30 W, 8 min for a 3.0 ~ 3.7 cm active tip at 50 W, and 10 min for a 4.0 ~ 4.7 cm active tip at 60 W, respectively.

All the operations were performed by three senior interventional radiologists (L.X., Y.W. and X.H. with over 10 years of experience in thoracic interventions under CT guidance). Depending on the location of the target nodules, patients were placed in a prone position, lateral position or supine position to ensure the best puncture site and entry route and avoid important structures, including ribs, interlobular fissures, and blood vessels. Lidocaine was administered at the puncture site to induce local anesthesia of the pleura. With CT monitoring, the radiofrequency electrode was punctured according to the predetermined direction and angle. The ablation was not performed until the CT scan confirmed that the electrode hooked the lesion. Considering tumor shape and size, one or two needle ablations with a constant antenna position were usually acceptable to achieve complete ablation. The operators strived to achieve ablation range greater than the lesions by at least 5 mm. If the intraoperative complications such as intra-alveolar hemorrhage (IAH) or atelectasis, made it impossible to determine the extent of ablation, at least 2 cycles of ablation would be performed to raise the impedance until ablation stopped. After completion of the RFA session, the ablation electrode was withdrawn, and a repeat CT (same parameters) scan was performed to evaluate whether the ablation zones covered the tumor and the occurrence of ablation-related complications, mainly including pneumothorax and hemorrhage.

Local efficacy was assessed by two radiologists who were blind to clinical data (H.C. and H.H. with over 5 years of experience) through chest enhanced CT examination at least 6 months after RFA according to mRECIST criteria (24, 41). If they had disagreements, it would be determined in consultation with the senior expert (W.L. with over 20 years of experience). The follow-up CT examination one month after ablation was taken as the baseline (42). Based on the mRECIST criteria, CR was defined if any of the following manifestations on CT were seen: the disappearance of the lesion, cavity, fibrosis or nodule without enhancement. If two consecutive CT examinations demonstrated the target lesions had irregular enlargement or enhanced solid components, they were classified as a non-complete response (non-CR).

In order to avoid the data bias due to the difference in scanning spacing and slice thickness between preoperative and immediately postoperative CT images, the following preprocessing steps were adopted: the CT images were uniformed to a common resolution of 1 mm × 1 mm × 1 mm by B-spline interpolation algorithm, and then the window width was adjusted within the range of - 1200 Hu to 600 Hu and the intensity was scaled within the range of 0 ~ 255. After normalization of all CT images, the samples containing pulmonary metastases were trimmed to 3D cubes with the size of 40 mm × 40 mm × 40 mm. Finally, the gray values of sample cubes were normalized between 0 and 1 ( Figure 2 ).

To objectively and accurately delineate the target lesions in the preoperative CT images and the boundary of the ablation area immediately after RFA, A 3D U-Net model (43) was used to segment the lesions and ablation region automatically, and two junior radiologists (H.C.and H.H.) verified the segmentations and made the necessary adjustments to guarantee the accurancy and repeatability ( Figure 3 ). If they had disagreements, it would be determined in consultation with the senior expert (W.L.).

For each segmented preoperative lesions and ablation region, 1252 radiomics features were extracted through the open-source feature toolboxes PyRadiomics (44) and PREDICT. The radiomics features were comprised of 13 intensity features, 35 shape features, 9 orientation features and 507 texture features which contained 144 Gray Level Co-occurence Matrix (GLCM) features, 16 Gray Level Size Zone Matrix (GLSZM) features, 16 Gray Level Run Length Matrix (GLRLM) features, 14 Gray Level Dependence Matrix (GLDM) features, 5 Neighborhood Grey Tone Difference Matrix (NGTDM) features, 156 Gabor filters features, 39 Laplacian of Gaussian (LoG) filters features, 39 Local Binary Patterns (LBP) features (32, 44–47), and 688 wavelet features.

In order to reduce unnecessary, redundant information and complexity in the process of calculation and modeling, the maximum correlation and minimum redundancy algorithm (MRMRA) (48) was used for features selection. There are five common variants under the MRMRA framework (49): mutual information difference (MID), mutual information quotient (MIQ), F-test correlation difference (FCD), F-test correlation quotient (FCQ), and random-forest correlation quotient (RFCQ). The formulas were as follows:

Assuming that there were m features in total, for a given feature Xi, i∈ (1, 2,…, m), the importance of the feature could be determined by MRMRA, commonly in the following five forms:

where, Y is the category label corresponding to the variable, S is the selected feature set, |S| is the size of the feature set, Xs∈S is a feature outside the feature set S, and Xi represents a feature that was not currently selected; the functionI(·,·) represents mutual information,p(·,·) is the Pearson correlation coefficient, F(·,·) is the Fstatistics, and I _RF (·,·) is the random forest feature importance score. Since inconsistent results of various methods under different super parameter conditions, we utilized the above 5 methods to filter features. The frequency of the top 5, top 10, and top 15 features was counted in the importance ranking, and the experiments were conducted from 5 to 15 features with the highest frequency to obtain the best performance, and eventually to confirm the 9 selected features.

As the Gaussian mixture model (GMM) had good performance in the evaluation of sample distribution and similarity in high-dimensional space (50–52), we used distance metric learning based on the Gaussian mixed model (DML-GMM) rather than random splitting to divide data according to our previous research results (53). We demonstrated that when the sample size was large, there was little difference between random splitting and the DML-GMM model. As for a smaller sample size, however, the DML-GMM model could obtain more stable results than random splitting. Therefore, the log-likelihood of the extracted radiomic features was calculated by DML-GMM model to describe the distribution, the data was split into multiple clusters and then was divided into 5 groups by stratified sampling. One group was selected as the testing set and the remaining 4 groups were used as the training set. Actually, we did five-fold cross-validation and chose a single split data including 386 lesions for the training set and 93 for the testing set.

Due to the unbalanced distribution of case counts in CR and non-CR, we adopted the oversampling method (synthetic minority over-sampling technique, SMOTE) (54) to mitigate the biased impact of data imbalance on the models during training.

Clinical model: all the clinical and radiological features were included in the univariate Logistic regression analysis, after which the variables with P< 0.1 were included in the multivariate analysis. Finally, the independent factors with P< 0.05 were selected for modeling. The random forest technique was a regression tree technique which utilized bootstrap aggregation and randomization of predictors to achieve a high degree of predictive accuracy (55). Since the random forest algorithm has been proven to be effective and superior in building clinical models (56, 57), our clinical model was also built on it.

Radiomics model: A random forest, which was the most common classifier used for radiomics features classification, contained multiple decision trees, and the total output result was determined by the subcategories of each decision tree. When processing high-dimensional data, it had a strong ability of anti-interference and anti-overfitting, especially for unbalanced medical data. Several studies have confirmed that the random forest model could be used to predict the survival rate, recurrence risk, and efficacy evaluation of lung cancer patients (58–60).Multimodal data fusion models: the random forest model was integrated based on radiomics and clinical models ( Figure 4 ). The weighted fusion strategy (61) adopted in our study was decision level fusion (late fusion) (62). This level of fusion allowed features from different representations to be combined in the same format of representation, which had more and better scalability and flexibility (63). The exact formula was confidence=ω ₁·confidence _image +ω ₂·confidence _clinical · The fusion prediction score was calculated to obtain the final prediction result.

In order to evaluate the performance of various models, we validated them on an independent test dataset, drew receiver operating characteristic curves (ROC) with the obtained prediction scores and calculated the corresponding area under curve (AUC). The difference in the predictive performance of models was compared by the Delong test (64). Meanwhile, the accuracy ( ACC = TP + TN TP + FP + TN + FN ) , sensitivity ( Sensitivity =   TP TP + FN ) , and specificity ( Specificity =   TN TN + FP ) were also calculated. In the formula, TP was true positive, FP was false positive, TN was true negative and FN was false negative.

Statistical analyses were performed using IBM SPSS (version 26.0, Chicago, USA). Man-Whitney U test was used for continuous variables which were presented as mean ± standard deviation (SD). Chi-square or Fisher test was used for categorical variables. All statistical tests were conducted at a two-sided significance level of P<0.05. All the medical image processing procedures and evaluation processes were performed on Python 3.6. In order to build the models and calculate the evaluation scores, we used publicly available packages such as SimpleITK, PyTorch, scikit-learn, numpy, and scipy.

A total of 198 patients with 479 lung metastases from CRC were enrolled; the detailed demographic characteristics are listed in Table 1 . After RFA treatment, there were 379 CR lesions and 100 non-CR lesions. Due to the small sample size, we analyzed each lesion individually in the same patient with multiple metastases as the recent literatures (20, 65–67). Through the GMM method, 386 lesions (305 CR and 81 non-CR) were selected to constitute the training dataset, and 93 lesions (74 CR and 19 non-CR) were chosen to constitute the independent testing dataset ( Figure 1 ). There were 227 lesions (47.4%)< 10 mm, and most lesions (399, 83.3%) were not close to the mediastinum or great vessels (diameter greater than 3 mm), but close to the pleura or diaphragm (287, 59.9%). The incidences of IAH and pneumothorax were 25.9% (124/479) and 24.0% (115/479), respectively.

Univariate logistic regression analysis in Table 2 showed that CEA, CA19-9, lesion location (including upper lobe of the right lung, right lower lobe, and left lower lobe), and intra-alveolar hemorrhage (P< 0.1) could completely identify ablated lesions. Furthermore, multivariate regression analysis demonstrated that CA19-9 (odds ratio [OR] = 1.007, P< 0.001) and lesion location (including right upper lobe [OR = 1, P = 0.005], right lower lobe [OR = 2.997, P = 0.003], and left lower lobe [OR = 2.498, P = 0.011]) were independent risk factors for incomplete ablation. These two clinical variables were used to construct a clinical model.

In order to prevent the model from overfitting because of the small sample size, 5 to 15 vital features with the highest scores were selected by MRMRA, and the five forms of MRMRA features importance scores were calculated separately, and compared with the default important feature s of the random forest model as the benchmark. The results demonstrated that the important features selected by MRMRA in the form of MID, MIQ, FCQ, and RFCQ had better performance than the features automatically selected by random forest, and the experimental model with 9 selected features had achieved better stability and smaller deviation. The selected feature results are shown in Table 3 .

The AUC values of each model were calculated in an independent testing dataset, and the DeLong test compared the corresponding P values ( Tables 4 , 5 and Figure 5 ). When radiomics features were integrated with clinical variables, and the coefficient of the radiomics model was 0.7 and the coefficient of clinical model was 0.3, the resulting AUC value achieved the highest (0.921) with the statistically significant difference (P values of 0.043) compared with the clinical model alone (0.830). In addition, the accuracy, sensitivity, and specificity of this multimodal data fusion model were also the best (82.6%, 80.3%, and 81.4%, respectively).

Figure 6 presents one example of a patient with post-lung RFA CR, with a nodule in contact with a vessel, complicated by IAH. In contrast, Figure 7 illustrates another example with post-lung RFA non-CR.