Added Value of Viscoelasticity for MRI-Based Prediction of Ki-67 Expression of Hepatocellular Carcinoma Using a Deep Learning Combined Radiomics (DLCR) Model

Simple Summary This study aimed to explore the added value of magnetic resonance elastography (MRE) in the prediction of Ki-67 expression in hepatocellular carcinoma (HCC) using a deep learning combined radiomics (DLCR) model. A total of 108 histopathology-proven HCC patients who underwent preoperative MRI and MR elastography were included. All the patients were divided into training and validation cohorts. An independent cohort including 43 patients was included for testing. A DLCR model was proposed to predict the expression of Ki-67 with conventional MRI (cMRI) as inputs. The images of shear wave speed (c-map) and phase angle (φ-map) derived from MRE were also fed into the DLCR model. Experimental results show that both c and φ values were ranked within the top six features for Ki-67 prediction with random forest selection, which revealed the value of MRE-based viscosity for the assessment of the tumor proliferation status in HCC. The model with all modalities (MRE, AFP, and cMRI) as inputs achieved the highest AUC of 0.90 ± 0.03 (CI: 0.89–0.91) in the validation cohort. The same finding was observed in the independent testing cohort with an AUC of 0.83 ± 0.03 (CI: 0.82–0.84). MRE-based c and φ-maps can serve as important parameters to assess the tumor proliferation status in HCC. Abstract This study aimed to explore the added value of viscoelasticity measured by magnetic resonance elastography (MRE) in the prediction of Ki-67 expression in hepatocellular carcinoma (HCC) using a deep learning combined radiomics (DLCR) model. This retrospective study included 108 histopathology-proven HCC patients (93 males; age, 59.6 ± 11.0 years) who underwent preoperative MRI and MR elastography. They were divided into training (n = 87; 61.0 ± 9.8 years) and testing (n = 21; 60.6 ± 10.1 years) cohorts. An independent validation cohort including 43 patients (60.1 ± 11.3 years) was included for testing. A DLCR model was proposed to predict the expression of Ki-67 with cMRI, including T2W, DW, and dynamic contrast enhancement (DCE) images as inputs. The images of the shear wave speed (c-map) and phase angle (φ-map) derived from MRE were also fed into the DLCR model. The Ki-67 expression was classified into low and high groups with a threshold of 20%. Both c and φ values were ranked within the top six features for Ki-67 prediction with random forest selection, which revealed the value of MRE-based viscosity for the assessment of tumor proliferation status in HCC. When comparing the six CNN models, Xception showed the best performance for classifying the Ki-67 expression, with an AUC of 0.80 ± 0.03 (CI: 0.79–0.81) and accuracy of 0.77 ± 0.04 (CI: 0.76–0.78) when cMRI were fed into the model. The model with all modalities (MRE, AFP, and cMRI) as inputs achieved the highest AUC of 0.90 ± 0.03 (CI: 0.89–0.91) in the validation cohort. The same finding was observed in the independent testing cohort, with an AUC of 0.83 ± 0.03 (CI: 0.82–0.84). The shear wave speed and phase angle improved the performance of the DLCR model significantly for Ki-67 prediction, suggesting that MRE-based c and φ-maps can serve as important parameters to assess the tumor proliferation status in HCC.


Introduction
Liver cancer is the fourth most common cancer and the second leading cause of malignant death. As the most common histological type, hepatocellular carcinoma (HCC) accounts for the majority of incidence and mortality of liver cancer cases [1]. HCC is highly heterogeneous and exhibits various biological behaviors. Ki-67 is one of the most common used cell proliferation biomarkers representing tumor aggressiveness. Previous studies showed that expression of Ki-67 is valuable in selecting recipients for liver transplant [2], patients with high Ki-67 expression are often encouraged to receive preoperative adjuvant therapy to improve prognosis [3], and Ki-67-targeted therapy is an attractive and promising avenue in HCC treatment [4]. Thus, the precise evaluation of Ki-67 expression preoperatively may be beneficial to appropriate clinical decision. However, current clinical practice guidelines recommend imaging to establish diagnosis of HCC noninvasively rather than through confirmatory biopsy. Therefore, the noninvasive preoperative evaluation of Ki-67 expression in HCC is warranted to outline personalized treatment strategies in clinical practice.
With the rapid development of deep learning models and radiomics theories, they have become effective methods of predicting tumor aggressiveness [5]. Preliminary studies showed the feasibility of MRI-based radiomics or deep learning models in predicting the Ki-67 expression of HCC. Current results were mainly based on the assessment of vascularity derived from Gd-DTPA-enhanced images. Hu et al. [6] reported that histogrambased parameters from apparent diffusion coefficient (ADC) maps and arterial phase (AP) images could be used to determine the Ki-67 labeling index (LI) in HCC with an AUC up to 0.826. Similarly, Fan et al. [7] proposed a combined model including an arterial phasebased RAD-score and serum alpha-fetoprotein (AFP) level to predict Ki-67 expression preoperatively. The combined model (AUC, 0.863), that included the AP RAD-score and serum AFP level, demonstrated improved performance when compared with the single AP radiomics model (AUC, 0.813). Recently, emerging evidence confirmed that hepatobiliary agent-enhanced MR imaging might provide valuable information to characterize HCC tumor biology independent of vascularity. Li et al. [8] found that texture analysis on the hepatobiliary phase (HBP), arterial phase (AP), and portal vein phase (PVP) images were useful in predicting Ki67 LI. However, the application of the hepatobiliary agent was limited by severe hepatic dysfunction or cholestasis, the relatively weak arterial phase hyperenhancement, and the relatively high frequency of arterial phase artifacts. In addition, interpreting the complex associations between the deep learning or radiomics features and biologic processes of HCC still remains an enormous challenge.
From the view of biomechanics, MR elastography (MRE) can provide insight into the mechanisms governing liver biology. Elastography applies mechanical tension or stimuli to tissue, monitors the tissue's response to the induced stimulus, and uses it to reconstruct parameters that characterize the mechanical properties of the image-encoded response. A variety of factors, including cell types, extracellular matrix deposition, cellularity, and fluid transport, alter the mechanical properties of biological organization [9,10]. MRE can also be applied on the liver [11,12]. MRE can measure displacement due to propagating mechanical waves, from which biomechanical properties of tumors, such as stiffness, can be calculated [13]. However, stiffness quantified by MRE was always used for depicting diffuse liver disease, such as liver inflammation, fibrosis, and portal hypertension in chronic hepatitis patients [11]. To update, nearly all elastography techniques measure stiffness without taking into account the viscoelastic, anisotropic, and nonlinear properties of most living tissues. The complex shear modulus of biological tissues is a composite dimension consisting of a storage modulus and a loss modulus, representing the elasticity and viscosity, respectively.
Tomoelastography, as an advanced MRE technique [14], yields quantitative maps of shear wave speed (c in m/s) and loss angle (ϕ in rad) as proxies of tissue stiffness and viscosity, respectively. These two parameters, especially viscosity, have demonstrated great potential for tumor detection and aggressiveness prediction of prostate cancer, pancreatic cancer, and liver cancer [12,15,16]. For malignant tumors, increased cellularity and aggressiveness not only leads to the higher stiffness but also changes the tumor microenvironment, which have liquid properties [17,18]. Moreover, tomoelastography provides full field-ofview (FOV) quantitative maps with spatial resolution comparable to conventional MR images, which is suitable for analysis by deep learning or radiomics models. To update, the biomechanical properties of HCC have not been investigated for aggressiveness evaluation.
Our study aimed to investigate the value of viscoelasticity derived from tomoelastography for predictive capability of Ki-67 expression in HCC using an MRI-based deep learning combined radiomics (DLCR) model.

Patients
This retrospective study was approved by our institutional review board and local ethics committee in Ruijin Hospital, Shanghai, China. Written informed consent was obtained from each participant. Between June 2020 and June 2021, 184 patients with suspected HCCs were initially enrolled to establish the predictive model. For patients with multiple lesions, the largest one was selected to avoid clustering effect. The inclusion criteria were as follows: (1) age ≥ 18 years; (2) patients who underwent conventional MRI and preoperative tomoelastography examinations. The exclusion criteria were: (1) poor image qualities (n = 3); (2) prior chemotherapy or radiotherapy (n = 24); (3) no histopathological results (n = 6) ( Figure 1). Finally, 108 patients (number of males, 93; number of females, 13; age, 59.6 ± 11.0 years) were included. All the demographics and clinical characteristics of patients were recorded (Table 1). They were randomly divided into a training cohort (76 males; 11 females; age, 59.4 ± 11.1 years) and a testing cohort (17 males; 4 females; age, 60.38 ± 10.20 years) with a ratio of 8:2. Whereafter, between July 2021 and November 2021, 43 HCC patients (36 males; 7 females; age, 60.01 ± 11.10 years) who met the above criteria were included as an independent validation cohort. The median time of interval between liver resection and MRI examination was 2 days, ranging from 1 to 3 days.

Tomoelastography
Tomoelastography was performed using a 1.5 T MR scanner (Magnetom Aera, Siemens, Germany). Continuous harmonic vibrations at frequencies of 30, 40, 50, and 60 Hz were induced by an external pressurized-air driving system. Two rear pneumatic actuators (0.6 bar) and one front pneumatic actuator (0.4 bar) were positioned near the liver. To ensure a stable mechanical state, apply vibration 5 s before imaging. Three-dimensional wavefields were acquired using a single-shot spin-echo echo-planar imaging (EPI) sequence with flowcompensated motion encoded gradients (MEG). Fifteen consecutive transverse slices of images were acquired with an FOV of 384 × 312 mm 2 and a resolution of 3 × 3 × 5 mm 3 (corresponding to a matrix size of 128 × 104 pixels) during the free breathing period. The other imaging parameters were a slice thickness of 5 mm, an echo time (TE) of 59 ms, a repetition time (TR) of 2050 ms, a GRAPPA parallel imaging factor of 2, an MEG frequency of 43.48 Hz for 30, 40, and 50 Hz vibrations and 44.88 Hz for 60 Hz vibration, and an MEG amplitude of 30 mT/m. The total time for the MRE measurements was approximately 3.5 min. Technical success of tomoelastography was evaluated by visual assessment of the shear wave images by an experienced operator.
Multifrequency wave field data were processed using the processing pipeline [19]. High spatial resolution maps of shear wave velocity (c) and phase angle (ϕ) over the entire FOV were obtained. Since c is proportional to the square root of the storage modulus, while ϕ continuously changes between 0 (solid properties) and π/2 (liquid properties), it can be considered to represent stiffness and fluidity, respectively.

Histopathological Analysis
Histopathological analysis was performed by a pathologist with 10 years of experience in liver pathology, who was blinded to the radiological and clinical findings. Immunohisto-chemistry staining for Ki-67 protein was performed. The Ki-67 expression was assessed by noting the percentage of positively stained cells with brownish yellow nuclei. Then, HCC was classified as low Ki-67 expression (Ki-67 ≤ 20%) and high Ki-67 expression (Ki-67 > 20%) groups.

Image Preprocessing
DCE and DWI images were registered to T2-weighted image using ANTs toolbox [20]. The B-spline interpolation was applied to the conventional MR images to normalize the image resolution to 0.9 × 0.9 × 4.0 mm 3 and FOV to 360 × 360 mm 2 . The tumor region of interest (ROI) was manually annotated by an abdominal radiologist (with five years of experience) on T2-weighted images, c-map, and ϕ-map using ITK-SNAP software [21].
Clinical data were also fed into the predictive model, including alpha-fetoprotein (AFP), platelet count, prealbumin level, alanine aminotransferase (ALT) level, aminotransferase aspartate (AST) level, total bilirubin, direct bilirubin, albumin level, etc. To avoid overfitting, LASSO regression (λ = 0.001) with 5-fold cross-validation was applied to reduce the dimensionality of features. Finally, the top six features were filtered out for further modeling. After features reduction, an SVM model (with penalty parameter = 1.0, kernel = 'rbf', order of polynomial kernel function = 3) was employed for the final prediction of high Ki-67 expression. The different network architectures (i.e., Inception-Resnet, Inception, Resnet, VGG16, VGG19, and Xception) were compared in the DLCR model based on conventional MRI (including T2, DWI, and DCE). After comparison, the best model was determined (Figure 2).

cMRI-Based DLCR Model with Tomoelastography
The cMRI-based DLCR model incorporated traditional MRI modalities. On this basis, both the c-map and ϕ-map were integrated to establish the model to compare prediction effects. To visualize the contributions of each feature, a random forest-based feature ranking was implemented (Figure 3).

Statistical Analysis
Continuous data were expressed as mean ± standard deviation. Categorical data were expressed as number and percentages. Receiver operating characteristic (ROC) analysis was conducted to evaluate the performance of two models in predicting high Ki-67 expression. Comparisons between the AUCs were conducted using the Delong test. A two-tailed Student's t-test with a p-value of <0.05 was considered statistically significant.

Demographics and Clinical Characteristics
In the training cohort, the high Ki-67 expression group included 40 patients while the low Ki-67 expression group included 47 patients. In the testing cohort, the high Ki-67 expression group included 9 patients while the low Ki-67 expression group included 12 patients. In the independent testing cohort, the high Ki-67 expression group included 17 patients while low Ki-67 expression group included 26 patients. Clinical characteristics were not significantly different between the training and validation cohorts. The detailed information is summarized in Table 2. represented significant difference between the two groups. The other demographic and laboratory parameters showed no significant difference between the two groups.

Optimization of cMRI-Based DLCR Models
When comparing the CNN models based on the six network architectures, the Xception architecture showed the best performance for classifying high Ki-67 expression, with an AUC of 0.80 ± 0.03 (95% confidence interval ( (Table 3). Thus, the Xception architecture was selected as the final deep learning feature extraction. By adding clinical data as input, the AUC achieved to 0.84 ± 0.03 (95%CI: 0.83-0.85) in the validation cohort and 0.74 ± 0.02 (95%CI: 0.73-0.75) in the independent testing cohort.

Comparison of DLCR Models with/without Tomoelastography
Based on the training cohort, both the cMRI-based DLCR model and the DLCR model with tomoelastography were established. When combined tomoelastography with cMRI and clinical data, the model achieved a higher AUC of 0.90 ± 0.03 (95%CI: 0.89-0.91) than the AUC of 0.84 ± 0.03 (95%CI: 0.83-0.85) for the DLCR model without tomoelastography in the validation cohort.

Contribution of Predictive Efficacy
Both c and ϕ values were ranked within the top six features for Ki-67 prediction with random forest selection (Figure 3b). The top two features were related to the viscosity parameter ϕ, which revealed the value of MRE-based viscosity for the assessment of the tumor proliferation status in HCC. Two representative cases with different Ki-67 levels are shown in Figure 4. The lesions with higher Ki-67 expression tended to have higher ϕ.

Discussion
Our study demonstrated that the MRE-based viscoelasticity map improved the performance of the DLCR model for Ki-67 expression in HCC. The results suggested that biomechanical properties, elasticity, and viscosity could provide additional information for the assessment of the HCC proliferation status. To the best of our knowledge, this is the first study to investigate contribution of biomechanics for HCC aggressiveness using machine learning methods.
According to the AUC, our model achieved better performance than the previously mentioned studies, and the reason can be explained because more imaging biomarkers were introduced in the predictive model. Our model used both deep learning and radiomics, and innovatively incorporated MRE. Compared to the models that only included conventional MRI and routine clinical data for the prediction of Ki-67 expression, the proposed DLCR model combined cMRI as well as MRE provided more information of the tumor characteristics. The tomoelastography that we used in this study provided high-resolution parameter maps, which can reflect the mechanical properties of the tissue accurately instead of only the total water content in the conventional MRI [28].
It is worth noting that several studies have found the value of elasticity in predicting the Ki-67 expression of breast cancer. Cha et al. [29] had found that higher elasticity value was associated with Ki-67 and the invasive size of the tumor in breast cancer. Choi et al. [30] discovered that Ki-67 positivity, high nuclear grade, large tumor (invasive) size, and so on, were associated with a significantly high shear wave elasticity ratio. However, its application in the prediction of HCC Ki-67 expression has not yet been investigated. In addition, those studies of breast cancer only measured the shear wave velocity by ultrasound, while we used two-dimensional tomoelastography in the DLCR model, which can better capture the textures for modeling and provide another vital biomechanical parameter ϕ, representing the viscosity of tissue.
MRE quantified the stiffness of tissues with structural and elastic characteristics, which was widely reported for tumor detection and classification [31,32]. However, the previous studies only reported the simple relationship between tumor aggressiveness and MREbased measurements, while we used c and ϕ maps as inputs of the DLCR model in order to get more detailed and spatial information about the tumor biomechanical properties. Our results showed that tomoelastography images, especially the ϕ-map, improved the performance of the DLCR model significantly. The features extracted by CNN and the zone variance value extracted by radiomics from the ϕ-map ranked the top two contributions in the DLCR model for the prediction task. In other words, the viscosity of HCC, which represents the fluidity properties, showed a good predictive ability of tumor aggressiveness.
The elevated fluidity was supposed to be associated with higher vessel density and the increased intracellular fluid mobility of tumor cells [30]. Higher Ki-67 expression is associated with faster progression and poorer prognosis for HCC patients. During tumor progression, increased vascularity and cellularity, a denser extracellular matrix, and higher interstitial pressure influenced the tumor stiffness and fluidity [33,34]. Moreover, constant elevating blood supply for maintaining tumor growth and development also had an impact on the tumor fluidity. Similarly, previous studies had found that the fluidity of tumors, which was related to angiogenesis and intercellular friction, have a significant contribution to the development and proliferation of malignant tumors [12,30]. In addition, solid tumors become aggressive by metastatic spread, which requires partial fluidification so that cancer cells can move. These findings show the close relationship between tumor fluidity and aggressiveness and can help to explain why the features from the ϕ-map ranked the top two contributions in the DLCR model.
Nevertheless, the present study had some limitations. First, the sample size was not large enough. Although we tested our model on an independent cohort, more validations are still needed, especially multicenter validations. Second, there is currently no standardized Ki-67 expression level threshold in HCC, which varies from 10% to 30% in studies. We treated the Ki-67 prediction task as a binary classification problem (with a threshold of 20%) according to previous studies [35,36] in order to make sure that the two groups had roughly the same number of people. A better solution might be dividing the Ki-67 expression into several subgroups. Third, the use of deep learning models was more targeted, which may not have widespread application. The complexity of the Inception architecture made it harder to make changes to the network. Xception performed better on the small dataset. As the number of layers increased, the VGG network showed a degradation problem. We only compared six widely used network architectures for the prediction task. More recently proposed architectures such as lightBGM can be used to further improve the prediction performance. Fourth, the morphological features from conventional MRIs, such as the LI-RADS classification [37] and other pathological features like MVI, are very important in evaluating the HCC aggressiveness and the patients' prognosis. Further work about the other prognostic factors for the HCC prognosis from the point of view of the biomechanical can be done.

Conclusions
In conclusion, our study proposed a noninvasive DLCR model and explored the added value of multifrequency MRE in the prediction of Ki-67 expression in HCC patients. The experimental results proved that the shear wave speed and phase angle improved the performance of the DLCR model significantly, suggesting that MRE-based c and ϕ-maps can serve as important parameters to assess tumor proliferation status in HCC.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cancers14112575/s1, Table S1: Scan parameter of the multiparametric MRI; Table S2: Parameter of the CNN networks. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patient(s) to publish this paper.
Data Availability Statement: Data presented in this study are available, by request, from the corresponding authors, due to matters of privacy.