1. Introduction

The liver is a frequent target for distant metastases from several tumors. Liver metastases (LM) are associated with poor prognosis and may occur early in gastrointestinal malignancies because of hematogenous spread through the portal venous system^{[1][2][3][4][5]}. Selected patients with LM, mainly those with liver-only metastases, can be considered for aggressive systemic and loco-regional treatments to prolong survival expectancy and optimize quality of life. Several studies have focused on LM from colorectal cancer, for which significant progress has been achieved. Effective chemotherapy may lead to a relevant improvement in survival, exceeding 30 months in the most favorable reports^[6][7][8]. Liver resection in selected patients obtained 5-year survival rates as high as 50%^{[9] [10][11][12]}; percutaneous ablation gained consensus, as it can grant effectiveness approaching that of surgery in small LM^[13]. The treatment of non-colorectal LM is also evolving, but therapies other than chemotherapy are still less codified^[14][15][16]. 

Such an aggressive policy, including several therapeutic options, requires a precision medicine approach. The selection of the appropriate course of action should rely on an adequate understanding of tumor biology and robust clinical biomarkers. However, the availability of reliable prognostic indices is currently an unmet need. Pathologic details of LM can be identified only ex-post after resection. Response of LM to chemotherapy is strongly associated with prognosis^[17][18], but it is overestimated by standard imaging modalities^[19][20][21]. Genetic mutations are promising biomarkers, but they are still under evaluation^[22][23].  

In recent decades, we became aware that imaging contains a great amount of data, namely in the form of grey level patterns, which are invisible to the human eye^[24]. These texture characteristics can be correlated with pathology data and outcomes^[25], potentially allowing diagnostic and prognostic evaluation. The analysis of textural features in medical images, which rely on mathematical functions, such as histogram analysis and matrices, is termed radiomics^[24][26]. Recently, radiomic features have been standardized by the imaging biomarker research initiative^[27]. This technology is attractive because it could be used to extract biological data directly from the radiological images, without invasive procedures, thus sparing costs and time and avoiding any risk for the patients. It would ideally embody the concept of "virtual biopsy." For many tumors, radiomic analyses have already provided an accurate evaluation of biology, allowing the identification of indices correlated with clinical outcomes ^{[28][29][30][31]}. In LM, where multiple therapeutic options are often available, a radiomics-based approach could be used to attain the most appropriate treatment decision. Based on the available literature, the present systematic review aims to elucidate the contribution of radiomic analyses to the management of patients with LM.

2. General Characteristics of the Studies

Figure 1 depicts the selection process. After screening for duplicates and eligibility, 32 studies were included in the qualitative synthesis. More than half of the publications (n = 18, 56%) were published in the last eighteen months. Most papers (n = 28) described retrospective analyses, while four reported planned secondary analyses of prospectively acquired data^{[32][33][34][35]}. Nineteen authors analyzed computed tomography (CT)^{[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]}, eight magnetic resonance imaging (MRI)^{[51][52][53][54][55][56][57][58]}, three positron-emission tomography (PET)/CT^[59][60][61], and two multiple imaging modalities (CT and MRI; PET and MRI, respectively)^[62][63]. Various software applications were used for texture analysis, with these being custom-made in a large proportion of cases (n = 10).

Figure 1. PRISMA Flowchart of study selection.

For the qualitative synthesis, we distinguished four groups of studies according to their subject: 1) radiomics of colorectal LM; 2) radiomics of non-colorectal LM; 3) capability of radiomics to perform differential diagnosis of focal liver lesions, distinguishing LM from other tumors (benign and malignant); 4) technical aspects of radiomics of LM. In the first group (radiomics of colorectal LM), we further distinguished four subgroups of studies according to their endpoints: prediction of survival, prediction of response to chemotherapy, correlation with pathological data, and miscellaneous. For details, refer to Section 4.3. Figure 2 summarizes the organization of qualitative analysis. Most papers (n = 18) analyzed radiomics of colorectal metastases. Due to the heterogeneity of studies, some papers fitted into more categories.

Figure 2. Outline of the selected studies organization. CRC: colorectal cancer, GLSZM: gray level size zone matrices, GLCM: gray level co-occurrence matrices, GLRLM: gray-level run-length matrices.

3. Assessment of Study Quality  

The average radiomic quality score (RQS)^[64] was in the 10 ± 6.5 (range 1–22), roughly 25% of the maximum score (n = 36). Only four studies (16%)^{[34][35][39][40]} had a score higher than 18 (>50% of the maximum score). The main limitations in quality were the following: no cost-effectiveness analysis (32 studies, 100%); lack of open-data repositories (n = 31, 97%); no phantom calibration (n = 31, 97%); failure to include a calibration statistic (n = 30, 94%); lack of prospective design (n = 28, 87%); and missing validation cohort (n = 18, 56%). At the transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD) checklist^[65] (31 elements), studies had an average score of 18 ± 3 points (range 14–29), i.e., 58 ± 10% of the maximum possible score. According to the quality assessment of diagnostic accuracy studies (QUADAS-2)^[66], there was a high risk of a patient selection bias in 34% of papers because of selection/inclusion criteria in most cases. One-fourth of studies had a high risk of bias related to the index test or the reference standard, while only one study (3%) had a high risk of bias in flow and timing. The RQS and TRIPOD scores of studies are reported in Table 1. Details of QUADAS-2 assessment and summary of its findings are presented in Table 2 and Supplementary Table S1.

Table 1. Details of studies included in the review.

Table 1. Cont.

CRC: Colorectal Cancer, NET: Neuroendocrine Tumor, EC: esophageal cancer, BC: breast cancer, PVP: portal venous phase, FS: fat suppression, DWI: diffusion-weighted images, SPAIR: spectral-attenuated inversion recovery, ADC: apparent diffusion coefficient, GLCM: gray level co-occurrence matrix, GLRLM: gray-level run-length matrix, GTDM: gray-tone difference matrix, GLSZM: gray level size zone matrix, GLNU: gray-level non-uniformity, NS: not specified.

Table 2. QUADAS-2 evaluation of studies.

Before analyzing the studies' results in detail, it is helpful to elucidate terminology that is commonly used in radiomics. The definition of radiomic features investigated in the studies is detailed in Table 3. In addition, region of interest (ROI) is defined as the selected area or volume of any imaging modality to analyze for the extraction of radiomic features.

Table 3. Overview of the analyzed radiomic characteristics.

2.3. Radiomics of LM From Colorectal Cancer

2.3.1. Prediction of Survival

Ten studies analyzed radiomic features' ability to predict the outcome of patients with colorectal LM [32,33,38,40,43,46,47,59,61,62]. Four included training and validation datasets [33,37,40,46]. The endpoint was overall survival (OS) in nine studies [32,33,38,40,43,46,47,59,61], progression-free survival (PFS) in four [46,59,61,62], hepatic PFS in one [47], and event-free survival in one [59]. In six papers, patients underwent chemotherapy [32,33,40,43,46,61], in two, liver surgery [47,62], and in two, either chemotherapy or surgery [38,59]. Six studies analyzed imaging modalities before treatment [38,43,47,59,61,62], while four performed comparative analyses of imaging modalities before and after chemotherapy [32,33,40,46]. Two studies performed texture analysis not only of LM but also of non-tumoral liver [38,47]. Finally, six studies analyzed the prognostic role of radiomic features in comparison or combination with traditional biomarkers [33,38,46,47,59,62].

Four studies demonstrated an association between the homogeneity/heterogeneity of LM and survival. Ravanelli et al. reported lower OS and PFS in patients with a higher uniformity of LM at CT scan (cut-off ≥0.42; relative risk (RR) = 6.94; 95% confidence intervals (95%CI) = 1.79–26.79 for OS; RR = 5.05, 95%CI = 1.74–14.66 for PFS) [46]. Andersen et al. described an association between shorter OS and tumor homogeneity at CT (hazard ratio (HR) ranging from 1.5 × 10²⁰ to 1.3 × 10⁴⁹, according to the filter used) [32]. Comparing imaging before and after chemotherapy, Dercle et al. identified a radiomic signature associated with OS based on two measures of heterogeneity (spatial heterogeneity and Graytone Difference Matrix contrast, HR = 44.3, 95%CI = 6.4–307.7 for patients with high imaging quality; HR = 6.5, 95%CI = 1.8–23.6 for patients with standard imaging quality) [40]. In the validation setting, the radiomic signature predicted survival better than KRAS-mutational status and 8-week tumor shrinkage evaluated according to RECIST criteria (AUC = 0.80 vs. AUC = 0.67 for KRAS and AUC = 0.75 for RECIST, p < 0.001). Finally, in the study by Rahmim et al., at multivariable analysis, LM heterogeneity at ¹⁸F-FDG PET/CT was a predictor of shorter OS (included into a predictive model, HR = 4.29, 95%CI = 2.15–8.57) [59]. The authors also depicted a model including histogram uniformity, number of metastases, and metabolic tumor volume that was predictive of a shorter event-free survival (HR = 3.20, 95%CI 1.73–5.94, p < 0 .001) [59].

Three studies showed an association between entropy and prognosis. Andersen et al. and Lubner et al. reported that the higher the entropy of LM, the better the OS (HR ranging from 0.16 to 0.63 in the Andersen et al. study, according to the filter used; HR = 0.65, 95%CI = 0.44–0.95 at coarse filter level in the Lubner et al. study) [32,43]. On the other hand, Beckers et al. reported some prognostic value of the ratio between entropy of LM and entropy of parenchyma (the higher the value, the shorter the OS, HR = 1.9, 95%CI = 0.95–3.78) [38].

Additional radiomic features have been reported. In the Simpson et al. study, LM correlation and contrast (combined into a single texture parameter) were associated with OS (HR = 2.35, 95%CI = 1.21–4.55) [47]. Dohan et al. analyzed imaging modalities before and after treatment and identified three predictors of OS: a decrease in the sum of the target liver lesions, high baseline density of dominant liver lesion, and drop in kurtosis [33]. Those three features (combined into a texture analysis score) evaluated after two months of chemotherapy had a strong association with OS (SPECTRA Score >0.02 vs. ≤0.02, HR = 2.82, 95%CI = 1.85–4.28 in the training dataset; HR = 2.07, 95%CI = 1.34–3.20 in the validation dataset). Radiomic score at two months had the same prognostic value of RECIST criteria after six months of chemotherapy. Shur et al. reported an association of minimal pixel value (negative prognostic factor, HR = 1.66, 95%CI = 1.28–2.16) and gray level size zone matrix (GLSZM) small area emphasis (positive prognostic factor, HR = 0.62, 95%CI = 0.47–0.83) with the PFS [62]. Finally, the following features have been associated with OS: standard deviation [32], LM density at CT scan [46], future liver remnant energy and entropy combined into a single linear predictor [47], ShapeSI4 (included in a radiomic signature) [40], and area under the curve of volume histograms at PET-CT [61].

The results of studies about radiomic features associated with the prediction of survival are summarized in Table 4 and Supplementary Table S2.

2.3.2. Prediction of Response to Systemic Chemotherapy

Ten studies analyzed the association between radiomic features and response to chemotherapy [32,33,37,38,40,45,46,54,56,61]. Evidence mainly regarded patients receiving anti-VEGF treatment. had training and validation datasets [33,37,40,46]. The reference standard was RECIST criteria in all but one study, which used tumor regression grade (TRG) [45]. Half of the studies analyzed the imaging modalities before and after treatment [32,33,40,45,46], while the other half considered only baseline imaging [37,38,54,56,61]. Four studies focused on targeted therapies associated with systemic chemotherapy [32,33,45,46].

In the Rao et al. paper, the entropy of LM after chemotherapy decreases in responders, while uniformity increases (entropy: −5.13 in responders vs. +1.27 in non-responders, OR = 1.34, 95%CI = 0.92–1.93; uniformity: +30.84 vs. −0.44, respectively, OR = 0.95, 95%CI = 0.89–1.01) [45]. Ravanelli et al. associated a good response with low baseline uniformity (cut-off ≥0.42; OR = 20, 95%CI = 1.85–217.4) [46]. In the study by Beckers et al., treatment success was slightly associated with higher entropy (6.65 ± 0.26 in responders vs. 6.51 ± 0.34 in non-responders, p = 0.08) [38]. The Zhang et al. analysis of T2 MRI images before chemotherapy showed that responding lesions had a higher variance and lower angular second moment (two measures of homogeneity) than non-responding ones (variance: 446.07 ± 329.60 in responders vs. 210.23 ± 183.39 in non-responders, p < 0.001; angular second moment: 0.96 ± 0.02 vs. 0.98 ± 0.01, respectively, p < 0.001) [56]. Dercle et al. built a signature, based on two measures of entropy, gray-tone difference matrix contrast and shape, which allowed to predict responsiveness to anti-angiogenic treatment (AUC = 0.80, CI95% = 0.69–0.94 for patients with high imaging quality; AUC = 0.72, CI95% = 0.59–0.83 for patients with standard imaging quality) [40]. Andersen et al. depicted LM modification after treatment with regorafenib. They observed data discordant with previous analyses (increase in entropy and decrease in uniformity), but none of the patients displayed a measurable response (85% had stable disease, while the remaining ones had progression) [32]. Considering skewness, in the study of Ahn et al., low baseline values (indicating a higher spread towards higher gray levels) were associated with response (0.02 ± 0.32 in responders vs. 0.33 ± 0.44 in non-responders, p = 0.001) [37]. One study demonstrated a skewness increase during treatment [32]. In opposition to CT and MRI, high entropy detected at ¹⁸F-FDG PET images before treatment predicted a worse response to therapy (AUC = 0.74, 95%CI = 0.52–0.97) [61].

Other features have been associated with response: high mean attenuation [37]; narrow standard deviation [37]; high baseline density of dominant liver lesion [33]; and mean values of histogram parameters for apparent diffusion coefficient maps [54].

The results of studies about radiomic features associated with the prediction of response to chemotherapy are summarized in Table 4 and Supplementary Table S3.

2.3.3. Prediction of Pathology Data

Three studies evaluated the association between radiomic features and pathology data [39,43,45]. Lubner et al. demonstrated that entropy, mean of positive pixels, and standard deviation are inversely associated with tumor grading (p = 0.007 for entropy, p = 0.002 for mean positive pixels, and p = 0.004 for standard deviation), while skewness and kurtosis showed a trend for an inverse association with KRAS mutation (p = 0.04 for skewness, and p = 0.058 for kurtosis) [43]. Cheng et al. reported that growth patterns of LM (desmoplastic, replacing, and pushing) can be successfully discerned on CT images by using second-order radiomic features, in particular gray level size zone matrix and gray level non-uniformity (AUC = 0.926, 95%CI = 0.875–0.978 in the training dataset; AUC = 0.939, 95%CI = 0.859–1.000 in the external validation dataset) [39]. In Rao et al.’s paper, the delta in entropy and uniformity values between pre- and post-chemotherapy imaging modalities were predictors of TRG values in patients receiving oxaliplatin-based chemotherapy with or without bevacizumab (entropy variation: −5.13 in TRG 1–2 vs. +1.27 in TRG 3–5, OR = 1.34, CI95% = 0.92–1.93; uniformity variation: +30.84 vs. −0.44, respectively, OR = 0.95, CI95% = 0.89–1.01), while RECIST criteria were not [45].

2.3.4. Other Papers

Three additional papers studied the radiomic features of colorectal LM. Reimer et al. analyzed the evaluation of response of LM undergoing trans-arterial radio-embolization [55]. In post-treatment MRI, higher kurtosis in arterial and portal phases and higher skewness in portal phase identified patients with a progressive disease earlier than standard RECIST criteria. Li Y et al. reported a model based on radiomic features of the primary tumor and LM before resection (heterogeneity, entropy, energy of vertical wavelet, emphasis) that was able to predict the future appearance of further LM [42]. Wagner et al. analyzed CT and PET-CT imaging of primary tumor and LM [60]. They demonstrated that colon cancer and LM have different skewness and kurtosis at both imaging modalities (CT and PET-CT), while colon cancers with or without LM have similar features.

Table 4. summarizes the data of studies dealing with colorectal LM.

Table 4. Studies on radiomics in patients with liver metastases from colorectal cancer.

Table 4. Cont.

Table 4. Cont.

R: Retrospective, P: Prospective, CRC: Colorectal Cancer, NET: Neuroendocrine Tumor, LM: liver metastases, RF: radiofrequency, XELOX: Capecitabine and Oxaliplatin, FOLFIRI: Fluorouracil and Irinotecan, FOLFOX: Fluorouracil and Oxaliplatin, TARE: trans-arterial radioembolization, OS: overall survival, HDFS: hepatic disease-free survival, EFS: event-free survival, PFS: progression-free survival, RECIST: response evaluation criteria in solid tumors, HGP: histological growth patterns, GLSZM: gray level size zone matrix, GTDM: gray-tone difference matrix, ADC: apparent diffusion coefficient, AUC-ISH: area-under-the-curve of cumulative SUV/Volume histograms, NS: not specified.

4. Radiomics of Non-Colorectal LM

Four papers focused on non-colorectal LM. A single study assessed CT-based radiomic indices in LM from esophageal cancer [41]. The study found that the characteristics of pre-treatment CT related to heterogeneity and gray-level intensity, such as wavelet gray level co-occurrence matrix correlation and gray level distance zone matrix with large dependence emphasis, were predictors of response to chemotherapy. Two studies explored radiomic analyses in LM from neuroendocrine tumors (NET)^[44][63]. Martini et al. analyzed a small series of patients (n = 49) but observed a number of associations: pancreatic NET had lower skewness and higher mean HU than non-pancreatic ones; entropy in the arterial phase was negatively associated with PFS in pancreatic NET and with OS in non-pancreatic NET; kurtosis was associated with lower OS in pancreatic NET, while skewness with higher one^[44]. Weber et al. investigated the correlation between parameters derived from the somatostatin receptor agonist (⁶⁸Ga-DOTATOC) PET and MRI with the proliferation index Ki67^[63]. Entropy and dissimilarity (from both PET and MRI) had a direct correlation with Ki67, while homogeneity had an inverse one. Moreover, it was possible to distinguish G1 and G2 LM on the basis of entropy, homogeneity, and dissimilarity (on PET data only). Finally, Trebeschi et al. reported heterogeneity-related radiomics parameters as predictors of response to immunotherapy in LM of melanoma and non-small-cell lung carcinoma^[49].

2.5. Differentiation of LM From Other Hepatic Lesions

Four studies investigated whether radiomic features could discriminate LM from other hepatic lesions. Jansen et al. analyzed metastases, primary hepatic tumors, and benign lesions (adenomas, cysts, and hemangiomas) on MRI images^[52]. A model using, among other features, the time to peak histograms and the sum of squared variance could distinguish different liver lesions. In the paper by Gatos et al., selected texture characteristics (inverse different moment, sum variance, and long-run emphasis) could differentiate metastases, hepatocellular carcinoma (HCC), and benign lesions^[51]. Li et al. tested a model to distinguish hemangiomas, LM, and HCC, using second-order features (gray level co-occurrence matrix, gray-level run-length matrix, and intensity-size zone matrix)^[53]. In their model, no feature combination could differentiate the three types of lesion at the same time. Differential diagnosis of the two malignant entities (LM and HCC) required a more complex model, with a higher number of features, than differential diagnosis between benign and malignant lesions (LM vs. hemangiomas or HCC vs. hemangiomas). Finally, a study by Song et al. identified kurtosis, variance, and inverse difference moment as distinguishing criteria between benign and malignant hypervascular lesions^[48]. Only the latter study used pathology data as the reference standard for all the analyzed patients.

Figure 3 provides an overview of the potential contribution of radiomics to the management of patients with LM.

Figure 3. Present and potential future contribution of radiomics to clinical practice.

5. Influence of Technical Features on Radiomic Analyses

Some studies set off to investigate whether acquisition or reconstruction parameters could influence the values of texture analysis indices. Ahn et al. tested three different reconstruction modes of CT images, i.e., filtered back-projection, iterative reconstruction model, and hybrid iterative reconstruction. The reconstruction method affected numerous parameters, including entropy, homogeneity, skewness, kurtosis, and gray-level co-occurrence matrices ^[36]. Lubner et al. compared the effect of 2D and 3D reconstruction on radiomic parameters by performing a Bland–Altman analysis on a subset of 20 patients^[43]. The results were similar for the two methods. Those results were confirmed by a further investigation by Ahn et al.^[37]. The latter study also compared the influence of different CT scanners, ranging from 8 to 64 rows, on the radiomic parameters, without finding any significant difference. Similar results were reported in the MRI setting: Peerlings et al. used the concordance correlation coefficient to test the reproducibility of an array of first- and second-level radiomics parameters over time (multiple MRI) and different MRI systems on apparent diffusion coefficient maps, finding good stability with most parameters^[58]. Conversely, two studies reported that radiomic parameters derived from CT scans are affected by slice thickness setting and ROI size^[34][35]. Dercle et al. demonstrated that ROI area size, metastatic site, and the individual characteristics of image acquisitions should be considered as confounding factors in the evaluation of tumor entropy^[35].

Inter-observer agreement was assessed by four studies [33,52,54,60]. Although they used different indices, such as K-statistic, intra-class correlation, and correlation index (r-value), all studies reported a substantial or excellent agreement among different readers. Finally, one study by Chatterjee et al. devised a method to reduce the rate of false discovery when analyzing radiomic parameters in small datasets^[57].