Prediction of Microvascular Invasion in Hepatocellular Carcinoma With a Multi-Disciplinary Team-Like Radiomics Fusion Model on Dynamic Contrast-Enhanced Computed Tomography

Objective To investigate microvascular invasion (MVI) of HCC through a noninvasive multi-disciplinary team (MDT)-like radiomics fusion model on dynamic contrast enhanced (DCE) computed tomography (CT). Methods This retrospective study included 111 patients with pathologically proven hepatocellular carcinoma, which comprised 57 MVI-positive and 54 MVI-negative patients. Target volume of interest (VOI) was delineated on four DCE CT phases. The volume of tumor core (Vtc) and seven peripheral tumor regions (Vpt, with varying distances of 2, 4, 6, 8, 10, 12, and 14 mm to tumor margin) were obtained. Radiomics features extracted from different combinations of phase(s) and VOI(s) were cross-validated by 150 classification models. The best phase and VOI (or combinations) were determined. The top predictive models were ranked and screened by cross-validation on the training/validation set. The model fusion, a procedure analogous to multidisciplinary consultation, was performed on the top-3 models to generate a final model, which was validated on an independent testing set. Results Image features extracted from Vtc+Vpt(12mm) in the portal venous phase (PVP) showed dominant predictive performances. The top ranked features from Vtc+Vpt(12mm) in PVP included one gray level size zone matrix (GLSZM)-based feature and four first-order based features. Model fusion outperformed a single model in MVI prediction. The weighted fusion method achieved the best predictive performance with an AUC of 0.81, accuracy of 78.3%, sensitivity of 81.8%, and specificity of 75% on the independent testing set. Conclusion Image features extracted from the PVP with Vtc+Vpt(12mm) are the most reliable features indicative of MVI. The MDT-like radiomics fusion model is a promising tool to generate accurate and reproducible results in MVI status prediction in HCC.

Objective: To investigate microvascular invasion (MVI) of HCC through a noninvasive multi-disciplinary team (MDT)-like radiomics fusion model on dynamic contrast enhanced (DCE) computed tomography (CT).
Methods: This retrospective study included 111 patients with pathologically proven hepatocellular carcinoma, which comprised 57 MVI-positive and 54 MVI-negative patients. Target volume of interest (VOI) was delineated on four DCE CT phases. The volume of tumor core (V tc ) and seven peripheral tumor regions (V pt , with varying distances of 2, 4, 6, 8, 10, 12, and 14 mm to tumor margin) were obtained. Radiomics features extracted from different combinations of phase(s) and VOI(s) were cross-validated by 150 classification models. The best phase and VOI (or combinations) were determined. The top predictive models were ranked and screened by cross-validation on the training/ validation set. The model fusion, a procedure analogous to multidisciplinary consultation, was performed on the top-3 models to generate a final model, which was validated on an independent testing set.
Results: Image features extracted from V tc +V pt(12mm) in the portal venous phase (PVP) showed dominant predictive performances. The top ranked features from V tc +V pt(12mm) in PVP included one gray level size zone matrix (GLSZM)-based feature and four first-order based features. Model fusion outperformed a single model in MVI prediction. The weighted fusion method achieved the best predictive performance with an AUC of 0.81, accuracy of 78.3%, sensitivity of 81.8%, and specificity of 75% on the independent testing set.

INTRODUCTION
Hepatocellular carcinoma (HCC), is the fifth most common cancer and third leading cause of cancer-related death worldwide (1). Despite advancements in medical technology leading to great achievements in the treatment of HCC, the prognosis of HCC remains poor, with a 5-year recurrence rate of 70% with hepatectomy and 35% following liver transplantation (2)(3)(4). Microvascular invasion (MVI) is defined as the presence of tumor cells in portal veins, large capsule vessels, or in the vascular space lined by endothelial cells (5,6). Evidence has shown that MVI is an independent predictor of recurrence and poor outcome following surgical hepatic resection (4,(7)(8)(9).
Currently, pathological examination is the gold standard for identifying MVI of HCC following operation or biopsy collection. This approach, however, is unreliable in case of sample contamination or for ineffectively done preoperative needle biopsies due to intratumoral heterogeneity (10). Furthermore, needle biopsy may increase the risks of unintended tumor bleeding or implantation metastasis (11). Therefore, there is a pressing demand for an accurate and non-invasive method for early prediction of MVI.
To date, the identified clinical indicators of MVI include: Alpha-fetoprotein (AFP), prothrombin induced by vitamin K absence or antagonist II (PIVKA-II), and other serological markers (12). However, the prediction performances of these markers were unsatisfactory due to poor sensitivity and specificity. Morphological and imaging characteristics such as large tumor size, unsmooth tumor margins, multinodular tumor morphology, rim enhancement on the arterial phase, and peritumoral hypointensity on the hepatobiliary phase of Gadoliniumexthoxybenzyl-diethylenetriamine-pentaacetic (Gd-EOB-DTPA)enhanced MRI, have also been associated with the presence of MVI (13)(14)(15)(16)(17)(18). These characteristics, on the other hand, are easily prone to inter-observer variations (19), as evidenced by the inconsistent interpretations of the conventional CT/MRI images for MVI prediction in previous studies (20,21).
Radiomics is a new method for disease diagnosis and prognosis prediction. Radiomics have exhibited great potential in predictive/discriminative models by integrating diseaserelated imaging features with clinical, pathological, and genetic data (22)(23)(24). Progress has been made in applying CT or MRI based radiomics for investigating MVI in HCC (25)(26)(27)(28). However, most studies only included a single phase of dynamic contrastenhanced (DCE)-CT or MRI. Exclusion of the other phases is likely to omit useful information which may impair MVI prediction. Also, research has focused mainly on the intratumoral region. Therefore, key image information beyond the tumor core might be lost since MVI often occurs in regions neighboring the tumor/nontumor interface (6,29). Image texture from the peripheral liver parenchyma, such as a settled 5mm or 10mm distance from the tumor edge, have shown encouraging predictive ability in MVI (26,28,30). To our knowledge, no study has comprehensively reported on MVI predictive performances using different combinations of DCE-CT/MRI phase(s) with varying distances from the tumor margin.
Predictive performance is also closely related to the prediction model (or classifier) used. Different classifiers are built on different mathematical models and thus generate inconsistent performances with the same classification task (31). An ensemble of classifiers, a process analogous to disease diagnosis by a multidisciplinary team (MDT), produces more reliable and accurate predictions compared with a single classifier (32)(33)(34). This study hypothesizes that this principle also applies with the fusion of various predictive models to yield enhanced and reproducible performance in MVI prediction.
The main study objective was to investigate the performance of radiomics analysis for MVI prediction in HCC. This study also aimed to identify the predominant phase(s) and the most relevant tumor periphery range for MVI prediction using noninvasive DCE-CT. Miscellaneous predictive models are established by considering different combinations of phases, tumor peripheral margins, feature selection methods and classifiers. The predictive performances of each model, as well as the final fusion model obtained through a multi-disciplinary team (MDT)-like fusion method, were to be explored.

Patient Cohort
This study was approved by the Institutional Review Board of Guangzhou First People's Hospital and the requirement for informed consent was waived based on the nature of a retrospective study. A total of 212 patients who underwent preoperative DCE-CT for newly diagnosed HCC from January 2016 to April 2020 at Guangzhou First People's Hospital were considered for inclusion in the study. The inclusion criteria were: 1) Pathologically confirmed HCC; 2) preoperative quadriphasic DCE-CT performed, and 3) complete preoperative lab tests. The exclusion criteria were: 1) Patients who had received anticancer therapy including chemoembolization, radiofrequency ablation, or transcatheter arterial chemoembolization (n=98); and 2) time interval between DCE-CT scan and surgery of more than two weeks (n=3). Finally, a total of 111 HCC patients (MVI positive: n=57 and MVI negative: n=54) were enrolled in this study.

Imaging and Histopathology
Preoperative DCE-CT were performed on multiple scanners with four phases following intravenous injection of the contrast agent, including phase 1-early arterial phase (EAP), 18-25 s; phase 2-late arterial phase (LAP), 35-40 s; phase 3portal venous phase (PVP), 50-60 s; and phase 4-equilibrium phase (EP) 120-250 s. The detailed imaging parameters are shown in the Supplementary Materials.
All surgical specimens were examined by one pathologist (W.S. Ding, with 14 years of experience in pathological diagnosis of hepatocellular carcinoma) to confirm the MVI status of the resected tumor.

The Volume of Interest Delineation
All images in each phase were stored in DICOM format and anonymized. Delineation of the target volume of interest (VOI) was performed by the ITK-SNAP software (http://www.itksnap. org) on the CT images slice-by-slice on phases 1-4 ( Figure 1). Visible tumor margins were first manually delineated to obtain the volume of tumor core (V tc ). This procedure was conducted by two investigators (W.L Zhang and R.M Yang, with 4 and 15 years of experience in radiological diagnosis, respectively) who lacked prior knowledge of the patients' MVI status. The conformity of delineated VOIs was measured by the Dice similarity coefficient. The two delineated VOIs with Dice index greater than 0.9 were averaged to yield the final VOI. Discrepancies on the lesion boundary (Dice < 0.9) were resolved by further discussions until mutual consensus were reached. The V tc was then extended to different distances (2; 4; 6; 8; 10; 12; 14mm) from the tumor margin, to obtain seven VOIs of the tumor periphery (V pt ), which were automatically generated with a morphological dilation algorithm. This process was not entirely isotropic as the expansion would stop on encountering large vessels (vessel caliber ≥2mm), bile ducts, or liver margin. All the manual steps allowed slight adjustment to acquire tailored VOIs for each phase. The VOI delineation was performed on the largest lesion for patients with multiple lesions.

Radiomics Feature Extraction and Analysis
Radiomic features were extracted from each VOI using an open-source python package Pyradiomics (https://pyradiomics. readthedocs.io/en/latest/index.html). There were 94 features in total extracted from the candidate features set including: 1) First order features (n=19); 2) gray level co-occurrence matrix (GLCM) features (n=24); 3) gray level size zone matrix (GLSZM) features (n=16); 4) gray level run length matrix (GLRLM) features (n=16); 5) neighboring gray-tone difference matrix (NGTDM) features (n=5); and 6) gray level dependence matrix (GLDM) features (n=14). Please refer to the Pyradiomics documentation (32)(33)(34) for their detailed definitions. Feature extractions were performed on each of the four phases, and the corresponding obtained features were combined and categorized into four groups: This resulted in a total of 15 types (4 + 6 + 4 + 1) of different combinations of features, which then served as the input for a specific predictive model.

Multi-Disciplinary Team-Like Prediction Modeling
The patient samples were divided chronologically into a training/ validation set (n = 88) and an independent testing set (n = 23). In this research, a typical prediction model was developed on a feature selection strategy and a classifier, and was then crossvalidated by a tenfold cross-validation (CV) using the training/ validation set (90% training, 10% validation). In each step of the ten-fold CV, a specific feature selection method screened out an optimal subset of the features to train a particular classifier. Fifteen feature selection methods and ten classifiers were investigated and their possible combinations resulted in 150 different prediction models.
The 15 types of features were extracted from different combinations of VOIs (V tc , V pt or V tc +V pt , note here V pt is obtained with seven different tumor periphery distances), and then were fed into each of the above mentioned 150 models. We yield totally 33750 (15×150×(1 + 7+7)) models to be compared on the training/validation set. The predictive powers of all the evaluated models were quantified by the area under the receiver operating characteristic (ROC) curve (AUC), accuracy (ACC), sensitivity (SEN), and specificity (SPE).
The best phase and V pt for MVI status prediction were determined by comparing the AUC values. The top-3 ranked models were identified by ten-fold CV on the training/validation set from the best phase and V pt . These top-3 ranked models were then fused by two fusion methods, i.e., the plurality voting (PV) and the weighted fusion (WF), to generate the final model that was then verified on the independent testing set.
Each of the top-3 predictive models was regarded as a clinical specialist providing a prediction of MVI. The PV gives a consensus prediction based on the highest number of votes.

Statistical Analysis
All statistical analyses were conducted using the SPSS 25.0 software (IBM SPSS Corporation, USA) and python 3.6.2 (Python Software Foundation (USA, https://www.python.org/ downloads/). Baseline patient characteristics were analyzed via univariate analysis. The categorical variables were presented as numbers and proportions and analyzed via the Chi-square test. Two-sided p values less than 0.05 were considered statistically significant. Comparisons between the 15 feature types were conducted using independent samples Kruskal-Wallis test with Bonferroni correction for adjusting for significant level in pairwise comparison.

Demographics
The study cohort comprised of 57 MVI-positive and 54 MVInegative patients who met the inclusion criteria. No statistically significant differences were seen in age, sex, presence of cirrhosis, hepatitis B and C virus infection, serum AFP, WBC, RBC, neutrophil count, Hb, ALB, PLT, PT, INR, AST, ALT, CB, TB, Scr, ALP, and Child-Pugh class in the training/validation and independent testing set ( Table 1).

Optimal Setting Determination
All the established predictive models were comprehensively compared to determine the optimal phases, tumor periphery, and VOIs combinations. Figure 2 shows the prediction comparison results on 15 phase combinations (F 1 pha , F 2 pha , …, F 1;2;3;4 pha ) and three VOIs combinations (V tc , V pt and V tc +V pt ). The AUCs were the mean values averaged for all the seven tumor peripheral distances and all the 150 classifiers. The best performance was seen in F 3 pha (portal venous phase, PVP) and this conclusion was consistent for both V tc (AUC=0.78) and V tc +V pt (AUC=0.82). Models inclusive of the PVP phase had the second (F 3;4 pha with V tc + V pt , AUC = 0.78) and third-best (F 2;3;4 pha with V tc +V pt , AUC=0.76) performances. In terms of VOIs combinations, the V pt generally achieved better predictive performance than the V tc in almost all phasic combinations (except for in F 3 pha andF 4 pha ). The V tc +V pt had improved prediction as compared with V tc or V pt alone in most of the 15 phase. combinations (except in F 2 pha , F 4 pha , F 1;3 pha , F 2 pha , F 2;3 pha , F 1;2;3 pha ). Using the optimal phase (F 3 pha ) and VOIs combination (V tc + V pt ), the optimal tumor peripheral distance was determined by comparing the predictions on F 3 pha and V tc +V pt with the aforementioned seven different peripheral distances. As shown in Figure 3, the predictive accuracy (in terms of the mean AUC averaged over all the 150 classifiers) increased gradually (maximal at 12mm) as larger tumor peripheral distance was involved.

Microvascular Invasion Prediction Performance
The training/validation set was cross-validated with features from F 3 pha and V tc +V pt (with 12mm tumor peripheral distance) on the 150 prediction models (classifiers + feature selection). The performances of all models were ranked, and the top-3 models were respectively combinations of "Random forest" & "t_score"; "Random forest" & "f_score" and "k-Nearest Neighbor" & "f_score". The predictive performances in terms of AUC of the top-3 models were 0.788, 0.776, and 0.775 on the training/validation set (with ten-fold CV), and 0.792, 0.78, and 0.803 respectively on the independent testing set ( Table 2). Fusion of the top-3 models further improved the predictive accuracy to AUC = 0.795 with the PV method and AUC=0.811 with the WF method. Comparison between the top-3 models and the two fusion methods were quantified by the NRI metric, with a positive NRI value indicating superiority. It was observed that PV or WF fusion outperformed prediction using any of the top-3 models (lower-left corner in Figure 4).

Top-Ranked Features
We also counted the number of times for each feature (with F 3 pha , V tc +V pt(12mm) ) being selected as the top features in the 150 prediction models using ten-fold CV ( Figure 5). The five most frequently selected features included: four first-order features (10 th percentile (20%), mean (11.31%), median (11.31%) and root mean square values (11.03%)) as well as one texture feature GLSZM-based Gray Level Non-uniformity Normalized (GLNN) (16.28%).

DISCUSSION
Image features extracted from the portal venous phase (PVP) with V tc +V pt(12mm) were shown to be the most reliable features for MVI prediction in HCC. An AUC of 0.81 on the independent testing set was achieved by the WF fusion method by integrating the top-3 models ranked from the training/validation set. Also, one GLSZM-based texture feature GLNN and four first-order features were found to be most associated with MVI.
Previous studies using conventional CT/MR imaging reported the common manifestations indicative of positive MVI, as pseudo-capsule, unsmooth tumor margins, rim or peritumoral enhancement in the arterial phase, and peritumoral hypointensity in the hepatobiliary phase. Chou et al. reported that unsmooth tumor margin which has been confirmed to be correlated with pathologically extra nodules, multinodular fusion, or infiltrative margin, is indicative of positive MVI in HCC lesions. A correlation between focal extra nodules on CT images and MVI in the pathologic specimens was also found (20). Research work by Matsui et al. suggested that peritumoral regions had a concentration of tumor drainage vessels which presented as corona enhancement in CT hepatic arteriography (CTHA) and CT arterio-portography (CTAP) (37). Nishie et al. reported that MVI positive HCCs, especially those with lesion diameter < 3 cm, tended to have larger area of peritumoral enhancement (due to peritumoral hemodynamic change) than MVI negative HCCs (38). Although conventional radiological manifestations are known to provide hints regarding MVI status in HCC, inconsistent results have been reported in previous studies (perhaps owing to interobserver variations) and thus standard diagnostic consensus has not been reached (16,20). Indeed, traditional radiological imaging analysis can be integrated with radiomics analysis for prediction modeling. However, radiographic diagnosis confirmed by naked-eye observation is usually limited by one's visual perception, which is insensitive to subtle image differences. Furthermore, diagnostic performance of traditional radiological imaging analysis is closely related to the radiologist's clinical experiences and can be easily biased by subjectivity. While radiomics provides an auxiliary alternative for radiologists to explore more hidden image patterns in characterizing diseases. Several studies have attempted to predict MVI status via radiomics analysis on DCE-CT or MRI (26-28, 30, 39, 40). For instance, Ma et al. extracted textural/non-textual features from the DCE-CT arterial phase (AP), PVP, and delay phase (DP) of the tumor core for MVI prediction. The PVP-based radiomics model was reported to achieve an AUC of 0.783 and 0.793 in the training and validation datasets, respectively (27). Similarly, Xu et al. compared radiomic features from the entire-volumetric interest (VOI entire ), 50% of the entire tumor volume (VOI 50% ) and a 5mm annular region neighboring the tumor surface (VOI penumbra ). It was shown that VOI entire and VOI penumbra (with AUC 0.841 and 0.829 in the training and validation sets respectively) outperformed VOI 50% in MVI prediction, which was consistent with this study (28). However, the combinational effect of VOI entire +VOI penumbra was not addressed. Nebbia et al. used multi-parametric MRI radiomics to predict MVI from the tumor core, at a fixed (10mm) peritumoral region, as well as the tumor core + peritumoral region (39). However, the fixed 10mm peritumoral margin was heuristically proposed which may omit MVI occurring in regions beyond this distance. The HCC pathological specimen collection standard may additionally support this as it includes liver tissue within a 2cm range from the tumor margin (6).
Previous investigations were conducted on a single (or two) DCE CT/MRI AP, PVP, or DP phase(s) and have not addressed  the combination of phases (26)(27)(28). This study is the first attempt, to our knowledge, to provide a comprehensive analysis aimed at substantiating the predominant role of the PVP phase (27). More diagnostic attention should be emphasized on PVP. Also, EAP had limited contributions to MVI prediction. Exemption of an EAP scan is therefore recommended in routine enhanced CT to reduce the radiation exposure. One GLSZM-based feature GLNN and four first-order features were shown to exhibit strong predictive capabilities of MVI. The GLNN measures the variability of gray-level intensity values in the image, with a lower value indicating greater similarity in intensity values (41). Our results implied that the MVI positive HCCs are associated with more heterogeneity within the V tc +V pt(12mm) than the MVI negative HCCs. This can be explained by the underlying HCC hemodynamic mechanism, that is, tumor cells may implant to the surrounding normal liver parenchyma from the tumor-draining vessels (37,38). Furthermore, the MVI positive groups demonstrated higher 10 th percentile, mean, median and root mean square values. Perhaps, this may be attributed to the hyper-attenuation resulting from excretion of CT contrast agent via the tumoral drainage vessels to the surrounding normal liver parenchyma in PVP.
This study also employed a novel method of using classifier fusion for MVI prediction modeling, such as a process analogous to disease diagnosis by a MDT, produce more reliable and accurate predictions compared with a single classifier. This was based on two reasons: 1) First, there is large performance differences between different classifiers, which has been confirmed in our previous study (31). It is not practical to select a suitable classifier for a given task from a large pool of classifiers, since different classifiers are built on different mathematical grounds; 2) Second, fusion of classifiers has been proved to generate more stable and reproducible classification performance than an individual classifier, and is effective in improving classification/prediction accuracy in decision-making (32,32,(42)(43)(44).

LIMITATIONS
Our study has several limitations. First, this was a retrospective and single-center study with a relatively small sample size due to an inclusion requirement of the EAP phase. However, EAP is seldom performed in the local institutions. Future investigations  should thus include more participants across different centers to confirm this study's findings. Secondly, manual VOI delineation is time-consuming and has uncertainties on subsequent radiomics analysis and prediction modeling (45). Semi-or automatic segmentation methods are expected to generate more consistent and reproducible results.

CONCLUSION
In conclusion, this study demonstrated the feasibility of noninvasive MVI prediction via CT radiomics analysis and a MDT-like fusion-based radiomics prediction modeling. Image features extracted from the portal venous phase on the tumor  core and within 12 mm of the tumor peripheral region may be considered as potential quantitative imaging biomarkers.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article AUTHOR CONTRIBUTIONS WZ, RY, and FL conducted the literature search. WZ, RY, XW, XZ, and XJ designed the study. WZ, GL, WD, and HW collected the data. WZ, RY, XZ, FL, SL, and AC analyzed the data. All