This retrospective study was approved by the institutional review board of our hospital, and the requirement of informed consent was approved for waiver. Patients were identified by searching through the picture archiving and communication system (PACS) database between September 2015 and September 2018. Preoperative gadobenate dimeglumine enhanced MRI were performed on 547 patients who were suspected of HCC. Two hundred and forty-two patients were initially excluded before SR or RFA due to tumor size (> 3 cm) ( Figure 1 ).

Patients were subsequently included according to the following criteria: 1) Patients with HCC confirmed by pathology or typical imaging findings (significant enhancement on arterial phase and wash-out in the portal venous or delayed-phase in multiphasic MRI); 2) Patients who underwent MRI examination within one month prior to SR or RFA; 3)Patients with complete clinical and laboratory data.

Exclusion criteria were as following: 1) Histopathologically diagnosed as other tumors rather than HCC (n=65); 2) Patients who underwent other treatments prior to SR or RFA (n=26); 3) Poor MR image quality, or lack of hepatobiliary imaging (n=32); 4) An insufficient follow-up time (<24 months) (n=87).

Patients were followed up regularly every 2-3 months for 2 years after SR or RFA and were monitored for recurrence by standard parameters including alpha-fetoprotein (AFP) level, liver function, CT and especially MRI. ER was defined as a new lesion with typical imaging characteristics in the remnant liver or other organs within 2 years after SR or RFA. Patients were followed until either ER or the end date of this study (September 30, 2020).

Clinical characteristics before SR or RFA were retrospectively obtained from an electronic medical record system. Demographic characteristics, history of hepatitis and cirrhosis, Child-Pugh score, AFP, alanine aminotransferase (ALT), aspartate aminotransferase (AST), and preoperative platelet count (PLT) were analyzed.

One month prior to SR or RFA, MRI was performed in all patients by using a 3.0 T MRI scanner (MAGNETOM Verio; Healthineers, Erlangen, Germany) with a dedicated phased-array body coil. The standard abdominal MRI protocol consisted of the following sequence: 1) axial T2-weighted fat-suppressed turbo-spin-echo (TSE): repetition time (TR)/echo time (TE), 4700/79 msec, slice thickness, 5 mm, slice gap, 1 mm, FOV, 21 mm×38 mm; 2) DWI (b=50, 800 sec/mm²) performed with a free-breathing single-shot echo-planar technique, TR/TE, 9965/73 msec, slice thickness, 5 mm, slice gap, 1 mm, FOV, 21 mm×38 mm. Corresponding ADC maps were automatically calculated by the MRI system; 3) multiphase contrast enhanced MRI, a 3D gradient echo sequence with volumetric interpolated breath-hold examination (VIBE), was performed before and after injection of Gadobenate Dimeglumine (MultiHance; Bracco, Shanghai, China) 0.2 ml/kg at a rate of 2 ml/sec followed by a 20 ml saline flush with the following parameters: TR/TE, 3.9/1.4 msec, slice thickness, 3 mm, slice gap, 1 mm, FOV, 25 mm×80 mm. Hepatic arterial phase (AP), portal venous phase (PVP), equilibrium phase images and hepatobiliary phase (HBP) were obtained at 20–30 sec, 70–80 sec, 180 sec and 90 min after contrast medium injection, respectively.

All MR images were independently reviewed by two experienced radiologists to assess the imaging features of the HCC in a single blind manner (unknown treatment for the patients and whether there were ER after SR or RFA). The two radiologists met to make final decisions by consensus when discordant cases occurred. The imaging features were selected according to the Liver Imaging Reporting and Data System (LI-RADS 2018) diagnostic algorithm including major features (nonrim arterial phase hyperenhancement (APHE), nonperipheral washout appearance (washout), enhancing capsule appearance (capsule), size and ancillary features (mild-moderate T2 hyperintensity, restricted diffusion, hepatobiliary phase hypointensity, etc.) (22). The maximum diameter of the tumor was measured at the level of the maximum cross-sectional area (if the lesion was clearly visible in other phases, do not measure in the arterial phase and DWI were avoided due to lesion size overestimation.

MaZda software (version 4.6, available at http://www.eletel.p.lodz.pl/mazda/) was used for texture analysis (TA) [24,25]. All MRI were transformed into Digital Imaging and Communications in Medicine (DICOM) format for compatibility.

All tumors were manually delineated by observer 1 (a radiologist with 7 years of experience in abdominal imaging), and the region of interest (ROI) was plotted on each cross section of the entire lesion. Texture features of the tumor were extracted from the T1 weighted, T2 weighted, AP, PVP and HBP image. HBP or T2 weighted imaging (in case of artifact) were used as the reference to delineate the tumor and were first segmented. Subsequently, the ROI was overlaid onto other phase images as required. If the tumor location had changed due to respiratory movement, fine adjustments were made to the ROI. The result of 101 features (derived from histogram, the absolute gradient, Gray Level Run-Length Matrix, Gray Level Co-occurrence Matrix, autoregressive model and wavelet transform) were generated from each three-dimensional segmentation, giving a total of 505 features for every lesion.

Twenty-six patients were randomly selected (12 ER and 14 NER) to explore the stability of each feature; observer 1 repeated tumor segmentation and observer 2 (with 9 years of experience in liver imaging) independently performed the segmentation in five image sequences (T2WI, T1WI, AP, PVP and HBP sequences). Intraclass and intergroup correlation coefficient (ICC) were used to evaluate the intra- and inter-observer repeatability of radiomics features, respectively. There is a good agreement of the feature extraction when the ICC is greater than 0.8.

To avoid overfitting and the curse of dimensionality, all radiomics features were loaded into the MAZDA feature selection package by image sequence (T2WI, T1WI and AP, PVP and HBP sequences) to select the most discriminative radiomics features between the ER group and the NER group. Feature selection algorithms included Mutual information [MI], fisher coefficient [Fisher], and classification error probability combined with average correlation coefficients [POE + ACC]. The three feature selection methods were supervised methods. Based on the MaZda’s automatic techniques, they were combined for the identification of 150 texture features in total, with the highest discriminative power for classification.

The least absolute shrinkage and selection operator (LASSO) method, which is suitable for the regression of high-dimensional data, was used to select the most useful predictive features among the final 150 texture features. LASSO is a regularization algorithm which can be used to eliminate irrelevant noises and do feature selection. It heavily relies on parameter λ, which is the controlling factor in shrinkage. The optimal value of λ is found by 10-fold cross-validation and 100 replicates to find the minimum mean squared error (minMSE) or minMSE + 1 standard error of minMSE backwards along the λ path (minMSE + 1SE), i.e., the largest λ-value such that the error is within 1SE of the minimum. The larger λ becomes, then the more coefficients are forced to be zero. Coefficients of some irrelevant variables are imposed to shrink towards zero, and every variable of which coefficient is non-zero is selected as the most significant predictor to be used in the model.

A combined model was created after analyzing potential factors in multivariate logistic regression analysis which included independent clinical risk factors, radiologic features and radiomics score. The nomogram based on the combined model was established to provide the clinician with a quantitative tool to predict individual probability of ER. To quantify the prognostic performance, area under the curve (AUC) for receiver operating characteristic (ROC) curve for the prediction model was then calculated with 95% confidence intervals (CIs). The optimal cutoff values from the maximum Youden’s index, as well as the corresponding sensitivity and specificity for discriminating ER and NER, were obtained from ROC curve analysis. The model was internally validated using 1000 bootstrap samples to avoid overfitting. Another validation of 5-fold cross-validation was also applied. The calibration curve provided a comparison between the expected and observed conversion probabilities.

Continuous variables were expressed as the mean ± standard deviation. The two-sample t-test or the Mann–Whitney U test for continuous variables, whereas the chi-square test or Fisher’s exact test was used as appropriate to compare the differences in categorical variables. Univariate and multivariate logistic regression analyses were performed to screen the independent risk factors of ER. Factors with a p value of 0.10 or less at univariate analyses were entered into the multivariate model. Odds ratios and 95% CIs were calculated. Model performance was assessed by model calibration. A p-value less than 0.05 (typically ≤ 0.05) is statistically significant. All statistical analyses were performed by using SPSS software (version 25.0; SPSS, Chicago, Ill) and R statistical software (version 3.6.3, https://www.r-project.org/).

A total of 111 patients including 93 (83.8%) male and 18 (16.2%) female were finally enrolled. These small HCC patients were treated with SR (n=45) or RFA (n=66) as initial therapy. A total of 53 (53/111, 47.7%) patients had confirmed tumor recurrence according to the two-year’ follow-up outcomes. Subsequently, patients were divided into the ER group (n=53) and the non-early recurrence (NER) group (n=58). Comparisons of baseline characteristics between small HCCs with or without ER are summarized in Table 1 . There was no difference in the rate of ER between SR and RFA (P=0.336).

There were significant differences in preoperative platelet count between the ER group and NER group (P=0.004). However, the other baseline characteristics did not differ between the two groups significantly. The optimal cut-off value was 157.5×10³/ml, corresponding to the maximum sensitivity and specificity of the ROC curve, where the AUC of the platelet count was 0.661 ( Figure 2 ). Patients with low preoperative platelet counts (<157.5×10³/ml) were significantly more likely to recur than those with high preoperative platelet counts.

The conventional MR imaging features between the ER and NER groups were described in Table 2 . We found that the enhancing capsule was statistically different between ER group and NER group (P=0.019). However, no significant differences were detected between the two groups across other imaging features, including tumor size, nonrim APHE, nonperipheral washout appearance, restricted diffusion, mild-moderate T2 hyperintensity, hepatobiliary phase hypointensity.

The interobserver ICC was >0.8, 0.5-0.79, <0.5 for 98%, 1% and 1% of the radiomic features extracted from five images sequences, respectively. The intraobserver ICC was >0.8, 0.5-0.79, <0.5 for 94%, 4% and 1% of the radiomic features extracted from six images sequences, respectively.

Nine potential features were selected by LASSO algorithm. These features were presented in the rad-score calculation formula. Features derived from non-enhanced T1 weighted images were reduced to 0, meaning a pretty poor robustness. The selected features in T2 weighted images included S(0, 0,1)Correlat. In AP images, the selected features included S(0,0,1)Correlat, Perc.01%3D, S(1,1,0)InvDfMom, GrNonZeros, S(0,1,0)InvDfMom, S(1,0,0)Correlat. And S(1,1,0)SumAverg in PVP images, Skewness3D in HBP images were also included. The rad-score for individual patients was calculated using the following formula: Rad-score = 5.670 + 0.487 × S(0, 0,1)Correlat (T2) - 5.272 × S(0,0,1)Correlat (AP) - 0.001 × Perc.01%3D(AP) - 4.071 × S(1,1,0)InvDfMom(AP) - 1.010 × GrNonZeros(AP) - 4.242 × S(1,0,0)Correlat(AP) - 0.901 × S(1,0,0)Correlat(AP) + 0.001 × S(1,1,0)SumAverg(PVP) + 5 × 10^-5 × Skewness3D(HBP). There were significant differences in rad-score between ER and NER patients (P=0.001), patients with ER generally presented higher rad-scores.

The area under the ROC curve of the rad-score was 0.979 and the optimal cutoff value was -6.15.

Table 3 summarized the risk factors found to be associated with ER in small HCC after univariate analysis, including the presence of liver cirrhosis (P=0.015), hepatitis B virus infection (P=0.083), presence of hepatobiliary phase hypointensity (P=0.089), the Child-Pugh score (P=0.083), preoperative platelet count (P=0.003), and the rad-score (P=0.001). No multicollinearity was found among all independent variables in the multivariate analysis. Eventually, multivariate logistic regression analysis further identified the rad-score and preoperative platelet count as the final independent predictors of ER of small HCC. The model that incorporated the aforementioned independent predictors was developed and further presented in the form of a nomogram ( Figure 3 ). The nomogram was used to provide clinicians with a quantitative tool for predicting the individual probability of ER.

The prediction ability of the final model was assessed using the AUC (estimated to be 0.981, 95% CI:0.957, 1.000, standard error 0.012) ( Figure 4A ) as well as the bias-corrected AUC, which was estimated using bootstrap with 1000 iterations and noted to be 0.980 (95% CI: 0.962, 1.000). The result of internal 5-fold cross-validation (AUC: 0.968, 95% CI: 0.916, 0.992, standard error 0.019) also showed favorable predictive efficacy. Figure 4B showed the calibration curve of the nomogram. The ideal curve fitted well with the calibration prediction curve, indicating the goodness-of-fit of the nomogram. Two case were provided to show nomogram’s ability to predict ER ( Figures 5 and 6 ).

The Nomo-score was calculated using the following formula: Nomo-score = 8.277 + rad−score × 1.336 – PLT(1, <157.5×10³; 0, >157.5×10³) × 1.767