Preoperative MRI radiomic analysis for predicting local tumor progression in colorectal liver metastases before microwave ablation

Abstract

Purpose

Radiomics may aid in predicting prognosis in patients with colorectal liver metastases (CLM). Consistent data is available on CT, yet limited data is available on MRI. This study assesses the capability of MRI-derived radiomic features (RFs) to predict local tumor progression-free survival (LTPFS) in patients with CLMs treated with microwave ablation (MWA).

Methods

All CLM patients with pre-operative Gadoxetic acid-MRI treated with MWA in a single institution between September 2015 and February 2022 were evaluated. Pre-procedural information was retrieved retrospectively. Two observers manually segmented CLMs on T2 and T1-Hepatobiliary phase (T1-HBP) scans. After inter-observer variability testing, 148/182 RFs showed robustness on T1-HBP, and 141/182 on T2 (ICC > 0.7).

Cox multivariate analysis was run to establish clinical (CLIN-mod), radiomic (RAD-T1, RAD-T2), and combined (COMB-T1, COMB-T2) models for LTPFS prediction.

Results

Seventy-six CLMs (43 patients) were assessed. Median follow-up was 14 months. LTP occurred in 19 lesions (25%).

CLIN-mod was composed of minimal ablation margins (MAMs), intra-segment progression and primary tumor grade and exhibited moderately high discriminatory power in predicting LTPFS (AUC = 0.89, p = 0.0001). Both RAD-T1 and RAD-T2 were able to predict LTPFS: (RAD-T1: AUC = 0.83, p = 0.0003; RAD-T2: AUC = 0.79, p = 0.001). Combined models yielded the strongest performance (COMB-T1: AUC = 0.98, p = 0.0001; COMB-T2: AUC = 0.95, p = 0.0003). Both combined models included MAMs and tumor regression grade; COMB-T1 also featured 10^th percentile of signal intensity, while tumor flatness was present in COMB-T2.

Conclusion

MRI-based radiomic evaluation of CLMs is feasible and potentially useful for LTP prediction. Combined models outperformed clinical or radiomic models alone for LTPFS prediction.

Keywords:

• Liver
• metastasis
• magnetic resonance imaging
• ablation
• precision medicine

Introduction

Colorectal liver metastases (CLM) strongly affect prognosis in colorectal cancer and are frequently diagnosed along the course of disease history [1,2]. Thermal ablation by either radiofrequency (RFA) or microwaves (MWA) has been endorsed by international guidelines in the management of CLM, along with other loco-regional treatments (LRTs) [3]. In this scenario, thermal ablation is selected for curative treatment of CLMs up to 3 cm. Despite its widespread application and its uprising as a possible non-inferior alternative to surgical resection, oncological efficacy is still limited by the occurrence of local tumor progression, i.e., the disease recurrence abutting the ablation zone after evidence of a complete coverage of the tumor area [4]. Numerous attempts have been made for elaborating predictive models of local tumor progression-free survival, whose main consistent predictors are tumor size and minimal ablation margins [5,6]; other factors such as mutational status and the impact of pre-ablation systemic chemotherapy have shown non-univocal results [7,8]; the net result is that local tumor progression cannot be fully explained by pure clinical approaches.

Radiomics, i.e., the quantitative analysis of volumetric imaging data, offers a potential means by which the composition of tumors may be explored, including possible microscopic features of tumor aggressivity [9–11].

As precision medicine approaches in oncology progress, the application of radiomic-derived data on prognostic models have been implemented in several solid tumors, yielding interesting results when integrated with traditional clinical data [12–14].

In the context of CLM ablation, application of radiomic-derived models has offered interesting data regarding prediction of LTPFS when studied on CT images [15].

As MRI has higher sensitivity than CT for detection of CLM as well as a multiparametric characterization of lesions [16,17], radiomic analysis of MRI images may have a greater potential also for their prognostic characterization.

The objective of this research is to assess the capability of MRI-derived radiomic features (RFs) to predict local tumor progression-free survival (LTPFS) in patients with CLMs treated with microwave ablation (MWA).

Materials and methods

Patient population and study design

An institutional prospectively collected database of all patients who received microwave ablation of CLMs (September 2015–February 2022) was retrospectively reviewed. Inclusion criteria for analysis were the following: oligometastatic colorectal cancer, ≤ 10 CLMs, ablation performed alone or in combination with surgical resection to achieve liver disease eradication, availability of preoperative (<30 days) MRI with gadoxetic acid. Patient with a follow-up shorter than 6 months or with evidence of residual disease at first follow-up imaging (30–40 days) were excluded.

Since the objective of the study was strictly related to LTPFS, which is a CLM-related endpoint, each CLM was analyzed independently in patients with multifocal disease.

Clinical variables

Preoperative clinical and radiological data were recovered for all patients in the cohort. Clinical data included primary tumor TNM staging, RAS mutational status, localization along the lower gastrointestinal tract, clinical risk score [18], history of previous CLM-directed therapies, history of systemic therapy. In patients undergoing concomitant resection and ablation, tumor regression grade (TRG) on histological specimens of resected metastases was also collected.

For each CLM, collected data included size, localization, synchronicity vs metachronicity relative to the diagnosis of primary CRC, radiomic-derived data on T1 and T2 images. Ablation time and power were recorded for each procedure and expressed as energy over tumor size (J/mm)

Image acquisition

MRI scans were acquired using a 1.5-T MRI (Achieva DStream, Philips Medical Systems, Best, The Netherlands) equipped with a multichannel coil.

The protocol included but was not limited to axial T2 weighted Tubo Spin Echo (TSE) Multivane XD images (2000 ms repetition time (TR), 100 ms echo time (TE), 400 × 400 mm field of view (FOV), 560 × 560 matrix, 5 mm section thickness, 0.5 mm gap, turbo factor = 39 and hepatobiliary-phase T1-weighted Fast Field Echo (FFE) fat sat DIXON (5,5 ms repetition time (TR), 3,8 ms echo time (TE),375 × 310 mm field of view (FOV), 400 × 400 matrix, 3 mm section thickness, obtained 20 min after injection of Gadoxetic acid disodium (Primovist, Bayer Schering Pharma, Berlin, Germany, 0.1 ml/kg) and 0.9% saline chaser (35 ml) administered as a bolus at a flow rate of 0.8 ml/s.

Ablation technique

All ablations were performed by experienced interventional radiologists either by the percutaneous approach or intraoperatively, as described previously by De Cobelli et al. [8]. A 2450 MHz/100W Microwave Emprint^TM ablation system with Thermosphere^TM Technology (Medtronic, Minneapolis, MN, USA) was used. The MWA protocol was tailored to tumor size to achieve complete tumor coverage seeking at least 5 mm minimal ablation margins.

Tumor recurrence and follow-up

As per institutional protocol, a contrast-enhanced CT was acquired 1 month after ablation to define technique efficacy; subsequently, follow-up consisted in cross-sectional imaging (CT/MRI) every 3 months for the first year, then every 6 months.

Minimal ablation margins were measured according to the methodology provided by Wang et al. [5]. First, the distance between the tumor and the nearest reliable landmark (including vessel bifurcations, benign lesions, capsula, surgical staples) was measured on preoperative imaging on all directions; then, the distance between the same landmark and the ablation zone was measured on the portal phase images of the CT obtained one month post-ablation.

The difference between the two distances was termed ‘Ablation Margin’; the smallest value among all the ablation margins calculated on all directions was termed ‘minimal ablation margin’ (MAM). MAM cutoffs of 5 and 10 mm were used in the analysis.

The definition of Local Tumor Progression (LTP) followed the standard terminology by Ahmed et al., i.e., appearance of disease foci abutting the ablation zone at follow-up [4].

Texture analysis

Segmentation – The segmentation of CLMs was performed by two independent observers (A.D.C, diagnostic and interventional radiologist, and F.C., senior resident in Diagnostic Radiology; the two readers had respectively 7 and 5 years of experience in interpretation of abdominal MRI images) on T2 and T1-Hepatobiliary phase (T1-HBP) sequences. Segmentations were performed manually with a dedicated software (MIM Software Inc., Cleveland, OH, USA). The dual manual segmentation was performed to allow for subsequent quantification of inter-observer variability. Images and segmentations were exported as DICOM files to proceed with subsequent analysis.

Radiomic features extraction – Import of DICOM files exported to Matlab platform (MathWorks, Natick, MA, USA)was performed through the Computational Environment for Radiological Research (CERR) [19,20]. Radiomic features were extracted with SPAARC Pipeline for Automated Analysis and Radiomics Computing (SPAARC) [21,22], following the International Biomarker Standardization Initiative (IBSI) guidelines [23]. Resampling of all images at cubic voxels of 1.5 × 1.5 × 1.5 mm³ with bilinear interpolation was performed directly from SPAARC settings pre-processing tools to minimize directional bias, according to IBSI recommendations. Image rebinning was implemented to optimize the RFs extraction process, as well as to minimize noise. Sixty-four bins were selected, in accordance with previous literature [24,25]. SPAARC allowed extraction of 182 first and higher order RFs: morphology, Statistical, Intensity Histogram (IG), Grey Level Run Length 3D_average (GLRL3D_avg), Grey Level Run Length 3D_combined (GLRL3D_comb), Grey Level Co-occurrence Matrix 3D_average (GLCM3D_avg), Grey Level Co-occurrence Matrix 3D_combined (GLCM3D_comb), Grey Level Size Zone Matrix 3D (GLSZM3D), Neighbors Grey Tone Difference Matrix 3D (NGTDM3D), Grey Level Distance Zone Matrix 3D (GLDZM3D).

Inter-observer variability – To test the robustness and reliability among features extracted from different contours the Intra-Class Correlation Coefficient was used (ICC) in accordance with other similar studies [26,27]. ICC is a metric based on the analysis of variance (ANOVA), ranging between 0 and 1 (which indicates, null and perfect reproducibility between observers, respectively). The ICC of a single feature was calculated from the feature values obtained from all contours and all patients, using a one-way random single-measure model. The threshold of ICC > 0.70 was considered as an indicator of ‘very good’ inter-observer agreement between, so that each radiomic feature with ICC <0.70 was discarded for subsequent analysis.

RFs redundancy limitation – A correlation-based filter was implemented in order to minimize the risk of redundancy: a Spearman coefficient threshold was set at 0.80 to discriminate redundant (≥0.80) and independent RFs (<0.80). The latter were included in further analysis; among the former, only the one with the lowest Univariate Logistic Regression p-value associated to LTPFS was selected.

Signal analysis– In this study all the scanner and the acquisition/reconstruction protocols used are the same for all patients. Consequently, it is likely that the signal normalization, which is a well-known issue in MRI images, can be overlooked. However, to exactly quantify the possible impact of a bias due to the signal variations on the radiomic features, two different ROIs were selected of both images’ series T1 and T2 representing a normal tissue as reference and a target. The first ROI was placed out of the tumor in the abdominal aorta, at the level of the L2 vertebra, where the signal is expected to be the same from each patient supposing there are not effects due to MRI scanner. Then another ROI was segmented in the lesion to investigate the impact of scanner also in the lesion. In addition, these structures have different range of gray-level intensities, covering different regions of the whole image histogram. Radiomic features were extracted from both regions and their variations around the mean value were quantified looking at the relative standard deviation. Subsequently, the changes of the values of individual features were analyzed in relation to fluctuations of the median intensity signal of the aorta. These trends were depicted in dispersion graphs because in cases of signal instability or other effects, a similar variation in the median signal from the aorta should also be reflected in the radiomic features extracted from the lesion itself, thus indicating an observable trend.

Statistical analysis

Standard descriptive statistics were utilized to describe patients and CLMs characteristics. Categorical variables were expressed as absolute number and presented as percentages; continuous variables were presented as mean ± standard deviation.

Machine learning-based bootstrap ranking was employed for variable selection, as previously described [9,13], to identify the best possible models among all possible variables combinations excluding variables that were found to be redundant and radiomic features not robust to segmentation, i.e., with ICC < 0.70. The Bootstrap was run 1000 times on the original dataset to obtain different populations with the same information content as the original. A maximum of three covariates were included in the Cox regression, in accordance with the number of events [28–30]. All the possible combinations of three variables were tested with Cox Logistic Regression on each bootstrapped sample, and for each triplet, the frequency with which the overall model P-value was found < 0.05 out of 1000 was calculated. A ranking of those appearing most frequently with overall model Pvalue <0.05 was done, and the combination with higher frequency was chosen as the most ‘robust’ model. The first combination in the ranking including only clinical variables was retained to develop the most reliable clinical model (CLIN-mod). The same was done selecting two models including only radiomic variables for the T2 and T1-HBP MRI sequences (RAD-T1; RAD-T2) and two combined models including radiomic and clinical variables together (COMB-T1; COMB-T2). In this way, the variables for the 7 models were selected showing the best robustness and performance in terms of p-value and AUC among all available combinations. Therefore, for each model the Cox Regression was finally applied on the variables selected on the original dataset. A prognostic index was calculated according to the formula of risk of Multivariate Cox Regression: $$\text{Pindex}=-\hspace{0.17em}\mathit{\text{ln}}\left(\frac{S}{\mathit{\text{So}}}{e}^{-{\displaystyle {\sum }_1^N\mathit{\text{Bi}}*\mathit{\text{Xi}}}}\right)={\displaystyle \sum _1^n\mathit{\text{Bi}}*\mathit{\text{Xi}}}$$

Where Bi are the coefficients estimated by Cox Regression and Xi are the variables included in the model, which range from 1 to N. It is represented by the linear combination of each variable included, weighted with its Cox Regression coefficient. To evaluate the capability of Pindex to discriminate CLMs according to LTPFS, a threshold value was calculated as the Youden criterion value of ROC analysis, then all metastases were dichotomized in two risk classes. Finally, the Kaplan-Meier test was run to compare the two risk categories. Models’ performances were quantified in terms of: C-index (that is the AUC value of the overall fit model), Sensitivity, Specificity, Kaplan-Meier p-value and hazard ratio. Analyses were performed using homemade Matlab codes.

To test the models given the small sample size, the ‘internal validation’ criterion suggested by the TRIPOD guidelines was followed [31]. The original sample was bootstrapped 1000 times again and subsequently, the P-indexes of the 7 models were calculated directly on the new populations and were tested with COX and ROC analysis on the 1000 samples to determine the median value and interquartile ranges of the p-value, C-index, and HR to assess the robustness of the results obtained in the original sample. It is worth noting that in this second bootstrap made for testing, the coefficients of the 7 models were fixed at those values obtained by the COX regression of the training. Therefore, they were simply tested on the 1000 samples to obtain the confidence intervals for the performance parameters obtained on the original dataset. However, in the first bootstrap used for selecting the best models, the coefficients of the Cox regressions were free.

Results

The study cohort included 43 patients with 76 CLM, consisting of 25 males (58%) and 18 females (42%). Mean age was 69 years (range 46 – 88); N = 6 patients were treated percutaneously (14%), whereas 37 (86%) were treated by a combined surgical-ablative approach during laparoscopy or laparotomy. Regarding notable CLM-related characteristics, mean tumor size was 12 mm, and minimal ablation margins (MAMs) were <5 mm in 23 nodules (30.3%). Population characteristics are summarized in Table 1.

Table 1. Population characteristics.

– Baseline characteristics
Patients (N = 43)
GENDER		AllRas mutation	10 (23.3%)	LOCATION OF PRIMARY TUMOR
Male	25 (58%)	AllRas wt	22 (51.1 %)	Cecum	3 (6.9%)
Female	18 (42%)	Missing	11 (25.6%)	Ascending Colon	5 (11.6%)
Age (years)	69 (46-88)	KRAS mut	9 (20.9%)	Transverse Colon	2 (4.7%)
>5 Lesions	22 (51.1%)	KRAS wt	23 (53.5%)	Descending colon	10 (23,3%)
EHD at time of ablation	14 (32.5%)	Missing	11 (25.6%)	Sigmoid Colon	10(23.3%)
INDICATION FOR ABLATION		NRAS mut	2 (4,6%)	Rectum	11 (25.6%)
Non-Surgical Candidate	3 (6.9%)	NRAS wt	30 (69.8%)	Missing	2 (4.6%)
Surgery Refusal	1 (2%)	Missing	11 (25.6%)
	Combined resection/ablation 35 (81.4%)			0	11 (25.6%)
TECHNICAL APPROACH				1	19 (44.2%)
Percutaneous	6 (14%)	Clinical Risk score		2	7 (16.3%)
Laparoscopic	14 (32.5%)	0	1 (2.3%)	Missing	6 (13.9%)
Laparotomic	23 (53.5%)	1	0	GRADE OF PRIMARY TUMOR
IMAGING GUIDANCE		2	6 (13.9%)	G1	0
Ultrasound	43 (100%)	3	19 (44.2%)	G2	22 (51.2%)
CT	0 (0%)	4	1 (2.3%)	G3	9 (20.9%)
LESIONS ABLATED PER SESSION		5	0	Missing	12 (27.9%)
Single	22 (51.2%)	Missing	16 (37.2%)	Previous Hepatic Resection	7 (16.3%)
Multiple	21 (48.8%)		NODAL STATUS OF PRIMARY	Post-ablation Chemotherapy	23 (53.5 %)
CLM (N = 76)
Size at ablation (mm)^†	11.8 ± 4.2 (4-30)	Subcapsular Location	10 (13.2%)	Synchronous	46 (60.5%)
Size at CSI	11.8 ± 2.8 (4-28)	Proximity to Vessels >3mm	25 (32.9%)	Metachronous	30 (39.5%)
Largest lesion size at any time	13.8 ± 7.1 (5-30)	Chemotherapy after lesion discovery and before ablation	54 (71%)	Time from lesion discovery to ablation (days)	155.1 ± 107.5 (0-779)
Previous LRT	41 (19.2%)	Time from last hepatic resection to ablation	406.7 ± 382.5 (5 – 1464)	Response to pre-ablation chemotherapy (RECIST1.1)
TARE	3 (1.4%)	Time from last chemotherapy cycle to ablation (days)	68,17 ± 60,55 (0 – 381)	PD	28 (29.5%)
RFA	28 (13.1%)			SD	11 (11.6%)
Resection	10 (4.7%)			PR	56 (58.9%)

After a median follow-up of 14 months, local tumor progression was observed in 19/76 CLM (25%).

After Inter-observer variability testing, 148/182 radiomic features (81.3%) proved to be robust on T1-HBP and 141/182 radiomic features (77.5%) on T2 (ICC > 0.7).

Concerining signal analysis, the relative standard deviation of the radiomic features included in the model extracted from the aorta was around 2%, indicating signal stability in blood. As expected, the relative standard deviation of the radiomic features included in the model extracted from the lesion was larger, up to 14%, as suggested by lesion heterogeneity by itself.

The changes of the values of individual features were analyzed in relation to fluctuations of the median intensity signal of the aorta. In cases of signal instability, a similar variation in the median signal from the aorta and the lesions is expected. However, no correlation trends were detected between the variation of radiomic features and those of the signal in the aorta, thus excluding signal instability. Dispersion graphs depicting these trends were reported in the Supplementary Materials.

Multivariate analysis

Following ranking procedure, the best clinical model (CLIN-mod) was identified and was composed of MAMs, intra-segment progression and primary tumor grade (Table 2). CLIN-mod exhibited moderately high discriminatory power in predicting LTPFS (AUC = 0.89, p = 0.0001), Figure 1. Upon Cox multivariate regression and Kaplan Meyer curve analysis, CLIN-mod categorized the risk of local tumor progression with a good performance [HR = 7.56 (95% CI = 2.42 − 23.61; p = 0.0005] (Figure 2).

Figure 1. ROC curves for prediction of local tumor progression according to clinical model (green), pure radiomic model (red), and combined model (blue) on T1-HBP images (figure a) and T2 images (figure b), respectively.

Figure 2. Kaplan Meyer curve for low-risk (green) and high risk group (red) as identified by the clinical model on local tumor progression.

Looking at the ranking procedure applied to imaging-derived data, two pure radiomic models (RAD-T1, RAD-T2) and two mixed models (COMB-T1 and COMB-T2) for T1-HPB and T2 sequence images were obtained, respectively.

RAD-T1 included IH₁₀, NGTDM3D-busyness and Morphology-COMshift (Table 3) and predicted LTPFS with quite good discriminatory power (AUC = 0.83, p = 0.0003), Figure 1, as well as a good performance in categorizing the risk of LTP on Cox multivariate regression and Kaplan-Meyer curve analysis [HR = 4.29 (95%CL = 1.54–11.96) p = 0.005, Figure 3].

COMB-T1 included minimal ablation margins, TRG and IH-percentile10 (Table 3) with a very strong discriminatory power for predicting local progression (AUC = 0.98, p = 0.0001, Figure 1). The Kaplan Meier analysis confirmed significant separation of survival curves in the groups identified by COMB-T1 (Figure 4, p = 0.0152)

Figure 3. Kaplan Meyer curve for low-risk (green) and high risk group (red) as identified by the T1-HBP-based (figure a) and T2-based (figure b) radiomic model on local tumor progression.

Figure 4. Kaplan Meyer curve for low-risk (green) and high risk group (red) as identified by the T1-HBP-based (figure a) and T2-based (figure b) combined model on local tumor progression.

RAD-T2 included GLSZM3D-smallZoneEmphasis, GLCM3D-sumEntropy, IH-entropy (Table 4) and showed a good discriminatory power for LTPFS prediction (AUC = 0.79, p = 0.001, Figure 1). Cox multivariate regression and Kaplan Meier analysis confirmed its ability in stratifying the risk of LTP [HR = 4.14 (95%CL = 1.52–11.23) p = 0.005, Figure 3]. COMB-T2 comprised MAMs, TRG and flatness (Table 4); it demonstrated a very high performance in discrimination of local progression (AUC = 0.95, p = 0.0003). The Kaplan Meier analysis confirmed significant separation of survival curves in the groups identified by COMB-T2. (Figure 4, p = 0.0078)

Regarding the internal validation, the bootstrap procedure confirmed the robustness of results as reported by Table S1, S2 and S3 of Supplementary Materials.

Table 2. Characteristics of the clinical model.

Clinical model
Variables	p-value (model)	p-value	HR	AUC	PPV	NPV	Accuracy	Sensitivity	Specificity	b	p-value (KM)	HR	95%CL
Minimal ablation margin >5 mm	0.0001	0.0009	0.16	0.89	0.5	0.94	0.76	0.85	0.73	−1.84	0.0005	7.55	2.41 − 23.61
Intrasegment recurrence		0.0062	5.27							1.66
Tumor grade (G)		0.085	2.96							1.08

Table 3. Characteristics of T1-based models.

Radiomic T1-HBP-based model
Variables	p-value (model)	p value	HR	AUC	PPV	NPV	Accuracy	Sensitivity	Specificity	b	p-value (KM)	HR	95%CL
IH10	0.0003	0.0333	0.87	0.83	0.34	0.95	0.64	0.87	0.59	−0.14	0.0054	4.29	1.54 − 11.96
NGTDM3D busyness		0.0037	1052.4							6.96
Morphology_COMshift		0.0699	0.01							−4.33
Combined T1-HBP-based model
Minimal ablation margin >5 mm	0.0001	0.0093	0.04	0.98	0.26	1	0.64	1	0.59	−3.27	0.0152	*_	*_
TRG		0.0017	0.17							−1.75
IH10		0.0980	0.85							−0.16

*The HR of the Pindex with the separation given by the Youden criterion of ROC analysis is not computable because of the sparsity of events in the lower risk class.

IH₁₀: 10^th percentile of histogram intensity; NGTDM3D: neighbors grey tone difference Matrix 3D; TRG: tumor regression grade.

Table 4. Characteristics of T2-based models.

Radiomic T2-based model
Variables	p-value (model)	p-value	HR	AUC	PPV	NPV	Accuracy	Sensitivity	Specificity	b	p-value (KM)	HR	95%CL
GLSZM3D-smallZoneEmphasis	0.0014	0.0026	0.000002	0.79	0.35	0.92	0.64	0.81	0.59	−13.24	0.0053	4.14	1.52 − 11.23
GLCM3D-avg-sumEntropy		0.0059	0.0040							−5.51
IHentropy		0.0264	317.48							5.76
Combined T2-based model
Minimal ablation margin >5 mm	0.0003	0.0118	0.02	0.95	0.26	1.00	0.63	1.00	0.57	−3.90	0.0078	*_	*_
TRG		0.0021	0.18							−1.71
Morphology-flatness		0.0790	0.001							−6.63

* The HR of the Pindex with the separation given by the Youden criterion of ROC analysis is not computable because of the sparsity of events in the lower risk class.

GLSZM3D: Grey level size zone matrix 3D; GLCM3D: grey level co-occurrence matrix 3D; IH: histogram intensity; TRG: tumor regression grade.

Discussion

In the present study, MRI-derived radiomic data showed a prognostic value when integrated to traditional predictive factors for local recurrence of colorectal liver metastases after microwave ablation.

After semi-manual contouring by two observers, more than 80% of extracted radiomic features exhibited a high concordance rate (>0.7) on both T2 and HBP-T1 sequences. This finding highlights the technical feasibility of MRI-radiomic analysis for colorectal liver metastases, since a large amount of radiomic features were retained as potential robust predictors. Most radiomic studies on CLM mainly focus on CT scan, mainly due to its higher availability. However, given the increasing role of locoregional treatments in management of CLM, clinical application of MRI in this setting has become predominant, therefore assessing prognosis from MRI-derived radiomics might have a stronger clinical impact. In particular, the relatively high concordance among observers may have been positively affected by the higher sensitivity of MRI in detection and characterization of CLMs, even when dealing with smaller lesions.

Furthermore, the possibility to assess different imaging sequences allows to achieve multiparametric tumor characterization also in the domain of texture features.

The potential superiority of MRI over CT has been partially demonstrated by Shur et al. when applied to tumor recurrence after surgical resection [32].

Mixed models were the best at predicting local tumor progression-free survival: in particular, COMB-T1 comprised MAMs, TRG and IH₁₀; in COMB-T2, the clinical covariates were retained, while IH₁₀ was replaced by flatness.

Minimal Ablation Margins, in particular below 5 mm, are a known to predict local tumor progression [6,33]. In our study, we followed the methodology described by Wang et al. to calculate MAMs [5]. Recent studies have emphasized the critical role of 3D softwares in the accurate assessment of ablation zones and margins post-procedure. These techniques could be implemented in future studies since they have been demonstrated to enhance the accuracy of these fundamental assessments, accounting for organ motion and deformation [34–37].

TRG is an interesting addition to the local prognostic model. Tumor regression grade reflects the degree of pathological response to systemic chemotherapy and was calculated on resected CLMs in all those patients who underwent a concomitant surgical and ablative intervention. Lower TRG indicates higher degree of response and has demonstrated to be associated with better clinical outcome in patients resected for colorectal liver metastases [38]. In our study cohort, the role of TRG in local tumor progression of concomitantly ablated CLMs may suggest a generally higher degree of local aggressivity of the oligometastatic disease.

In the combined model derived from T1-HBP images, IH₁₀ was negatively related to the event; in other words, a tendency toward lower signal intensity is associated to a metastasis more susceptible to local progression. Heterogeneous enhancement on HBP images has been described in colorectal liver metastases [39] and its occurrence may be at least partially explained by expression of organic anionic transporting polypeptide 1B3 (OATP1B3) negatively impacting prognosis [40]. In our study cohort, metastases did not qualitatively exhibit enhancement given their relatively small size (mean 1.2 cm), but this result may unveil a tumor microstructure more prone to local aggressiveness.

In the T2 scenario, flatness is the most predictive radiomic feature. This radiomic feature, ranging from 0 to 1, expresses the ratio between maximal and minimal CLM diameter. A perfectly spherical volume has flatness equal to 1, whereas smaller values of flatness express increasingly flatter objects. As a result, since this study showed the flatness negatively correlated to the events, it suggests that the flatter the volume the higher the risk of local progression .This result is sound considering the utilization of thermosphere technology which produces a spherical ablation zone and may result in better covering of more spherical tumors, regardless of their size [41].

Pure radiomic models achieved a good prognostic performance, even though they did not perform as well as the pure clinical one. This could be explained by the presence in the clinical model of intra-segment progression, i.e., appearance of tumor foci in the same liver segment but not abutting the ablation zone. Since intra-segment progression is observed during follow-up, it cannot function as an actual predictor, but its relationship with LTP, as proposed by a previous study [8], may be a sign of a locally aggressive behavior of the tumor.

This study has several limitations. Despite the internal validation reported, the monocentric nature of the study, with limited sample size, requires replication in other data sets for external validity. Furthermore, the retrospective nature of the study implied high degree of variability in potentially confounding variables (e.g., previous chemotherapeutic regimens), as well as inconsistent availability of potentially prognostic variables (e.g., RAS mutational status, carcino-embronic antigen (CEA) at the moment of ablation).

Conclusion

The application of T1-HBP and T2-derived radiomic data for local prognosis of CLM after microwave ablation is a feasible and has a potential to provide an added value to a clinical prognostic model. An external validation in a larger cohort is warranted to confirm consistency of the results. If validated, a perspective application of the mixed model could be an interesting add-on to address therapeutic options when considering patients for a local therapy of oligometastatic disease. The large amount of data coming from CLM texture characterization may potentially be the object of correlation studies with other data coming from novel technologies, including characterization of tumor genome, immune landscape and the application of metabolic imaging biomarkers, particularly through the use of 18 F-FDG PET/CT, which has already demonstrated added value for evaluation of ablation efficacy and treatment guidance [42–45]. Multi-level integration of such diverse data may eventually allow the extraction of even more complex and informative prognostic models.

Ethical approval

This study was performed in line with the principles of the Declaration of Helsinki. Approval was granted by the Ethics Committee of IRCCS Ospedale San Raffaele (Date: 12/09/2019; No: 137/INT/2019).

Consent to participate

Informed consent was obtained from all individual participants included in the study.

Author’s contributions

All authors contributed to study conceptualization. Methodology, data curation and formal analysis were performed by Angelo Della Corte, Martina Mori, Francesca Calabrese, Domenico Santangelo and Alessandro Pellegrini. The first draft of the manuscript was written by Angelo Della Corte and Martina Mori and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Acknowledgments

Stephanie Steidler contributed to language editing.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author, A.D.C, upon reasonable request.

Additional information

Funding

The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.