Texture Analysis Using Semiquantitative Kinetic Parameter Maps from DCE-MRI: Preoperative Prediction of HER2 Status in Breast Cancer

Objective To evaluate whether texture features derived from semiquantitative kinetic parameter maps based on breast dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) can determine human epidermal growth factor receptor 2 (HER2) status of patients with breast cancer. Materials and Methods This study included 102 patients with histologically confirmed breast cancer, all of whom underwent preoperative breast DCE-MRI and were enrolled retrospectively. This cohort included 48 HER2-positive cases and 54 HER2-negative cases. Seven semiquantitative kinetic parameter maps were calculated on the lesion area. A total of 55 texture features were extracted from each kinetic parameter map. Patients were randomly divided into training (n = 72) and test (n = 30) sets. The least absolute shrinkage and selection operator (LASSO) was used to select features in the training set, and then, multivariate logistic regression analysis was conducted to establish the prediction models. The classification performance was evaluated by receiver operating characteristic (ROC) analysis. Results Among the seven prediction models, the model with features extracted from the early signal enhancement ratio (ESER) map yielded an area under the ROC curve (AUC) of 0.83 in the training set (sensitivity of 70.59%, specificity of 92.11%, and accuracy of 81.94%), and the highest AUC of 0.83 in the test set (sensitivity of 57.14%, specificity of 100.00%, and accuracy of 80.00%). The model with features extracted from the slope of signal intensity (SIslope) map yielded the highest AUC of 0.92 in the training set (sensitivity of 82.35%, specificity of 97.37%, and accuracy of 90.28%), and an AUC of 0.79 in the test set (sensitivity of 92.86%, specificity of 68.75%, and accuracy of 80.00%). Conclusions Texture features derived from kinetic parameter maps, calculated based on breast DCE-MRI, have the potential to be used as imaging biomarkers to distinguish HER2-positive and HER2-negative breast cancer.


INTRODUCTION
Breast cancer is one of the most common cancers in women, and breast cancer alone accounts for 30% of new cancer cases in females (1). The status of human epidermal growth factor receptor 2 (HER2) is a biological factor that influences breast cancer survival. The 5-year relative survival rate has increased to 91% largely due to improvements in treatment, such as aromatase inhibitors for hormone receptor-positive tumors and trastuzumab for HER2-positive tumors. Of patients with hormone receptor-positive tumors, 81% receive hormonal therapy (2). With the development of new therapeutic drugs, better responses are seen with more specific pharmaceutical treatment options based on different molecular markers (3,4). Therefore, it is very important to accurately identify HER2 status to individualize treatment. At present, HER2 amplification status is determined by immunohistochemistry (IHC); tumors are considered to be HER2-positive if the IHC analysis is scored as 3, whereas tumors are considered to be HER2-negative if scored as 0 or 1. For tumors with IHC scores of 2, further analysis by fluorescence in situ hybridization (FISH) is needed to detect the amplification status of the HER2 gene. However, these methods require invasive biopsies, and are also subject to sampling errors due to intratumoral heterogeneity (5). Moreover, FISH examination is costly and time-consuming. Thus, it would be clinically beneficial to develop a cost-and time-effective, accurate, noninvasive method to detect HER2 status.
Breast dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is the most widely used and clinically proven imaging technique in breast cancer, and it has high sensitivity for detecting breast lesions. It provides anatomical information as well as hemodynamic information of the tumor with a high spatial resolution (6). Previous studies on breast DCE-MRI have indicated that morphological and kinetic characteristics were associated with benign and malignant tumors, response to neoadjuvant chemotherapy, and histopathological factors of breast cancer (7)(8)(9)(10)(11).
Texture analysis has been widely applied to characterize the spatial distribution of gray level intensities in images, capturing image patterns which are usually unrecognizable or unresolved by the human eye (12)(13)(14). This approach aims to extract highthroughput information to characterize image heterogeneity in specific target regions (15,16). The most commonly used texture features can be layered by the statistical order of the voxel information encoded within the target regions, including firstorder (also called histogram), second-order (gray level cooccurrence matrix and run-length matrix), and high-order (structural and transformed) texture features, proving to be helpful in assessing tumors. Earlier studies on rectal cancer revealed that texture features were useful for prediction of pathological complete response after neoadjuvant chemotherapy (17)(18)(19). Moreover, histogram features have been shown to be useful in evaluating tumor heterogeneity in glioma and cervical cancers (20,21).
In previous studies, texture features derived from mammography and multidetector computed tomography images have been applied and shown to potentially identify HER2 status in patients with breast cancer (22,23). However, DCE-MRI is recognized as the most common and effective method in breast cancer imaging. Montemurro et al. (10) showed that Fischer's score, which included three morphological, two functional, and five DCE-MRI features, was inversely associated with HER2-overexpression. Another study demonstrated that texture features from DCE-MRI were predictive of HER2 status (24). Semiquantitative kinetic parameter maps provide a technique for leveraging the preand post-contrast acquisitions, and can reflect kinetic information for breast cancer. A recent study demonstrated that the model based on texture features from semiquantitative kinetic parameter maps was able to discriminate sentinel lymph node status (25). To the best of our knowledge, no previous study has investigated the association between HER2 status in breast cancer and texture features extracted from semiquantitative kinetic parameter maps calculated from breast DCE-MRI.
Thus, the aim of this study was to evaluate whether features derived from semiquantitative kinetic parameter maps could be used to identify HER2 status in patients with breast cancer.

Study Population
This retrospective study was approved by our institutional review board (NO.2019PS175K) and the requirement for informed consent was waived. From January 2019 to January 2020, female patients with histologically confirmed breast cancer who underwent breast DCE-MRI were reviewed with our picture archiving and communication system (PACS). The inclusion criteria were as follows: (1) visible breast lesion on DCE-MRI; (2) histologically confirmed breast cancer; and (3) exact HER2 amplification status determined by IHC/FISH examination. The exclusion criteria were: (1) patients who underwent a biopsy before MRI examination; (2) patients who received neoadjuvant chemotherapy before MRI; and (3) insufficient MRI quality due to obvious motion artifacts.
Finally, a total of 102 patients were enrolled in this study retrospectively. Of these patients, 48 were HER2-positive and 54 were HER2-negative. The clinical characteristics collected using the PACS included age, maximum tumor diameter, estrogen receptor status, progesterone receptor status, Ki-67 status, histological grade, and histological type. Patients were randomly divided into a training set (n = 72, 34 HER2-positive and 38 HER2-negative) and a test set (n = 30, 14 HER2-positive and 16 HER2-negative) at a proportion of 70% and 30%, respectively. Figure 1 shows the workflow of this study.

MRI Acquisition
All patients received a pretreatment breast DCE-MRI at our institution using a 3.0 Tesla MR scanner (Ingenia, Philips Medical System, Best, Netherlands) equipped with a dedicated 7-channel bilateral breast coil with patient in a prone position. First, an axial fat-saturated T1-weighted precontrast scan based on the VIBRANT-VX technique was acquired. Then, eight axial contrast-enhanced fat-saturated T1-weighted scans were acquired after the intravenous bolus injection of a contrast agent (Magnevist, Bayer Healthcare Pharmaceuticals, Berlin, Germany) with a dose of 0.15 mmol per kg body weight. The imaging parameters were as follow: repetition time, 4.14 ms; echo time, 2.10 ms; flip angle, 12°; slice thickness, 2.00 mm; spacing between slices, 1.00 mm; field of view, 340 × 340 mm 2 ; matrix, 380 × 380. Eight subtraction sequences were obtained by subtracting the precontrast scan from each of the eight postcontrast scans.

Image Processing and Semiquantitative Kinetic Parameter Calculation
Two breast radiologists, each with over 6 years of experience, were blinded to HER2 status of patients and invited to help review the images. Slices with the maximum tumor diameter were chosen in consensus. The third-phase subtraction image, the eight phases of postcontrast images, and the precontrast image of the slice were downloaded and used for subsequent processing.
The lesion area was first delineated automatically on the third-phase subtraction image using an in-house software programmed with MATLAB 2018a (Mathworks, Natick, MA, USA). Seven semiquantitative kinetic parameter maps were calculated on the lesion area, respectively. The seven kinetic parameters included the initial percentage of enhancement (E initial ), the percentage of peak enhancement (E peak ), the early signal enhancement ratio (ESER), the maximum slope of increase (MSI), the second enhancement percentage (SEP), the signal enhancement ratio (SER), and the slope of signal intensity (SI slope ). The calculation formulas of the parameters are as follows: where SI 1 and SI 0 represent the signal intensities of the first postcontrast image and the precontrast image, respectively.
where SI peak represents the peak signal intensity value of the contrast enhancement.
where SI 2 represents the signal intensity at the second postcontrast time point.
where SI i and SI i+1 stand for the signal intensity of a certain phase and the following phase respectively, with i ranges from 0 to 7.
where SI 8 is the signal intensity at the eighth postcontrast time point.
where SI mean is the mean value of the signal intensity at the first two postcontrast time points.

Texture Feature Extraction
All texture feature extraction was performed using an in-house software developed in MATLAB 2018a. Fifty-five texture features were derived from each kinetic parameter map, including histogram features, gray level co-occurrence matrix (GLCM) features, gray level run-length matrix (GRLM) features, and discrete wavelet transformation (DWT) features. The details of these features are provided in Table 1. GLCM parameters were calculated from four GLCMs corresponding to a distance of one pixel and four angles (0°, 45°, 90°, 135°), and the mean value of each feature over the four GLCMs was utilized. GRLM parameters were calculated from four GRLMs corresponding to four angles (0°, 45°, 90°, 135°), and the mean value of each feature over the four GRLMs was utilized. DWT parameters were calculated for two layers and three directions (horizontal, vertical, diagonal) to produce low and high frequency components. For example, harr_L represented the low frequency component using harr wavelet, and harr_DH2 represented the diagonal high frequency component of the second layer using harr wavelet.

Model Construction and Statistical Analysis
The clinical characteristics and kinetic parameters of the patients were statistically analyzed using SPSS 22.0 (IBM, Corp). Categorical variables included estrogen receptor status, progesterone receptor status, Ki-67 status, histological grade, and histological type, and these variables were compared between HER2-positive and -negative groups using the chisquare test or Fisher's exact test. For quantitative data including age, maximum tumor diameter, kinetic parameters, and texture features, the independent sample t-test was utilized when the data was normally distributed with homogeneous variance, and the Mann-Whitney U test was used when the data was not normally distributed. A two-sided P value less than 0.05 was considered statistically significant. The data from 72 patients in the training set were used for feature selection and model construction. Feature selection was performed using MATLAB 2018a. Separately for each kinetic map, Pearson's correlation analysis was first performed among features in the training set. Highly correlated features with coefficients greater than 0.95 were marked, and the ones with higher correlations with other features were removed. Then, the least absolute shrinkage and selection operator (LASSO) was used to select features with nonzero coefficients among the remaining features by 10-fold cross-validation. After removal, the features were randomly divided into 10 groups. At each feature selection loop, one group of features was chosen as the validation set and the remaining groups were used as the training set. The optimal subset of features for prediction was generated after each loop, and this process was repeated for all ten folds. All selected features were recorded for further analysis.
The multivariate logistic regression analysis using forward stepwise selection was applied with entry of the selected features to establish the prediction model. Spearman's correlation analysis was performed to evaluate the correlation between texture features contained in the model and HER2 status. The performance of the trained model was assessed through the area under the receiver operating characteristic (ROC) curve (AUC). The sensitivity, specificity, and accuracy were calculated correspondingly. The optimal threshold was chosen according to the maximum Youden index. The established prediction model was further tested on the test set using the same threshold determined on the training set. The corresponding AUC, sensitivity, specificity, and accuracy were also calculated.

Characteristics of the Study Population
A total of 102 patients (51.60 ± 10.10 years) were included in this study. The detailed clinical and histopathological characteristics between HER2-positive and -negative groups are listed in Table 2. There was no statistical difference between the two groups with respect to age (P = 0.57), maximum tumor diameter (P = 0.26), histological grade (P = 0.17), or histological type (P = 0.91). The two groups showed significant differences in terms of estrogen receptor status (P < 0.01), progesterone receptor status (P = 0.02), and Ki-67 status (P = 0.04). Figure 2 shows two randomly selected cases used to display the results of lesion segmentation along with seven semiquantitative DCE maps and corresponding pathological results.

Performance of the Prediction Model
The comparison results of the average value of seven kinetic parameters in the lesion area between HER2-positive and -negative groups is provided in Table 3. There were no significant differences in the average value of seven kinetic parameters between the two groups. Table 4 presents the logistic regression models obtained from the training set. Table 5 shows comparison results of texture features included in models in the training set between HER2-positive and -negative groups. Short Run Emphasis derived from E initial , ESER, MSI, and SER maps were significantly different between HER2-positive and -negative patients (P < 0.01, < 0.01, < 0.01, and < 0.01, respectively). Contrast (P < 0.01) and harr_HH2 (P < 0.01) from E peak maps, autocorrelation (P < 0.01) from SEP maps, and gray level nonuniformity (P < 0.01) from SI slope were also significantly different between the two groups.
The performance of the prediction models is summarized in Table 6. Among the seven prediction models, models with features extracted from the ESER map yielded an AUC of 0.83 in the training set [95% confidence interval (CI), 0.72-0.91; sensitivity of 70.59%, specificity of 92.11%, and accuracy of 81.94%], and the highest AUC of 0.83 in the test set (95% CI, 0.64-0.94; sensitivity of 57.14%, specificity of 100.00%, and accuracy of 80.00%). The model with features extracted from the SI slope map yielded the highest AUC of 0.92 in the training set (95% CI, 0.84-0.97; sensitivity of 82.35%, specificity of 97.37%, and accuracy of 90.28%), and an AUC of 0.79 in the test set (95% CI, 0.59-0.91; sensitivity of 92.86%, specificity of 68.75%, and accuracy of 80.00%). The corresponding ROC curves of the models with features extracted from the seven kinetic parameter maps are shown in Figure 3.

DISCUSSION
In this study, the correlation between texture features and HER2 status in breast cancer was investigated using a texture analysis of seven semiquantitative kinetic parameter maps based on breast DCE-MRI. The results demonstrated that texture analysis based on DCE-MRI images has the potential to discriminate HER2   status in breast cancer. HER2 is a cell surface receptor expressed in normal breast cells that controls growth, division, and repair of breast cells (26). HER2-positive breast cancer is considered an aggressive disease, because the amplification of the HER2 gene results in an abnormally high amount of HER2 gene expression and HER2 proteins per cancer cell. Therefore, HER2-positive cancers promote the rapid growth and division of cancer cells, and the prognosis is generally poor (27). Trastuzumab treatment can be beneficial for breast cancer with HER2 amplification and overexpression, and therefore the HER2 status serves as a guide for treatment and is a crucial indicator of prognosis (28).
In recent years, there have been some studies on the association between radiomic features and HER2 status (22-24, 26, 29-31). Several studies investigated the relationship between HER2 status and radiomic features in gastric cancer and colorectal cancer (29-31). Li et al. (29) built and validated a CT-based radiomics nomogram for HER2 status prediction, which showed good performance. Zhou et al. (22) reported that mammography radiomics features can effectively diagnose HER2 status of patients with breast cancer, most notably with a model built using a combination of features from cranial caudal and mediolateral oblique views. One study indicated that    radiomics features from multidetector computed tomography images were associated with HER2 status in patients with breast cancer (23). Another study developed support vector machine models based on radiomic features from fat-suppressed T2weighted images and DCE-MRI and, using a combination of these features, noninvasively evaluated the HER2 status of patients with breast cancer. The model based on the combination of fat-suppressed T2-weighted images and DCE-MRI exhibited the best performance for predicting the HER2 status of patients with an AUC of 0.86 and an accuracy of 79.5% in the primary cohort, and an AUC of 0.81 and an accuracy of 78.3% in the validation cohort (24). However, no previous studies have explored the relationship between HER2 status in breast cancer and radiomic features derived from semiquantitative kinetic parameter maps based on breast DCE-MRI.
In the present study, the association between texture features from semiquantitative kinetic parameter maps and HER2 status in breast cancer was investigated. Fifty-five texture features were extracted from each of the seven semiquantitative kinetic parameter maps for each patient. Logistic regression models using forward stepwise selection were developed and validated to predict different HER2 status in breast cancer. Among the seven prediction models based on texture features from E initial , E peak , ESER, MSI, SEP, SER, and SI slope maps, two of the prediction models showed relatively good performance. The model built using features from the ESER map yielded an AUC of 0.83 in the training set, and the highest AUC of 0.83 in the test set. The model with features extracted from the SI slope map yielded the highest AUC of 0.92 in the training set, and an AUC of 0.79 in the test set. The texture features selected in the seven models included mean, kurtosis, contrast, autocorrelation, gray level nonuniformity, short run emphasis, high gray level run emphasis, and three DWT features. Contrast represents local variations presented in maps. Autocorrelation detects repetitive patterns of texture elements. Gray level nonuniformity measures the similarity of the gray level throughout the lesion area. Short run emphasis reflects the distribution of short runs. Compared with HER2-negative breast cancer, semiquantitative kinetic parameter maps of HER2-positive breast cancer showed higher contrast, autocorrelation, and gray level nonuniformity, as well as lower short run emphasis in the training set. The manifestation of these features indicated that semiquantitative kinetic parameter maps of HER2-positive breast cancer may show more heterogeneity and higher texture complexity than HER2-negative breast cancer.
Fusco et al. (32) calculated 10 semiquantitative kinetic parameters, maximum signal difference (MSD), time to peak between wash-in and wash-out segments, wash-in slope (WIS), wash-out slope (WOS), wash-in intercept, wash-out intercept, area under the curve of wash-in, area under the curve of washout, area under the curve of wash-in and wash-out, and standardized index of shape (SIS) as well as 50 textural features to predict breast cancer therapy response. The results demonstrated that SIS achieved the highest AUC value (0.93), suggesting that the joint feature from semiquantitative parametric maps may obtain the best diagnostic performance. In our study, we evaluated the texture features based on seven independent semiquantitative kinetic parameter maps.
A recent study showed that deforming autoencoder convolutional neural networks based on 3TP (three-channel images representing a given slice at three different time points, uniquely identified by means of the three time points) slices of DCE-MRI could be developed to discriminate malignant from benign lesions, and good diagnostic performance was achieved (33). In comparison, the performance of deep learning features was not investigated in our research, as we focused solely on the feasibility of the texture features from semiquantitative kinetic parameter maps. In addition, MRI scans in our study included more phases for post-contrast acquisitions, based on which of the seven kinetic maps were obtained. To improve the performance of the prediction model, further work should be conducted to develop the classification models by combining texture and deep learning features from kinetic maps of DCE-MRI for preoperative prediction of HER2 status in patients with breast cancer.
This study had several limitations. First, the sample size in this study was relatively small and this may have impeded the generalizability of the findings. Second, our results may not be applicable in all other institutions as our study was performed in a single institution with uniform MRI parameters. Additional studies are needed to increase cohort size and consider various conditions. Third, only the texture features extracted from semiquantitative kinetic parameter maps were used to discriminate different HER2 status in breast cancer in this study. Other clinical and histopathological characteristics such as volume, tumor location, lymph-vascular invasion, and diffusion-weighted imaging radiomics may also be good signatures for distinguishing positive and negative status of hormone receptors in breast cancer (25,34). Combining texture features and clinical characteristics in models may improve prediction performance. Fourth, the slices with maximum tumor diameter were selected and utilized in our study. Texture analysis was performed on two-dimensional images, and the representation of the entire volume of the tumor may have been limited compared with threedimensional analysis. Finally, Piantadosi et al. (35) reported a U-shaped deep convolutional neural network that exploited the well-known 3TP approach for the automatic lesion segmentation task, which showed a good result in breast DCE-MRI segmentation. However, our study used the semi-automatic segmentation based on Otsu's algorithm, which was timeconsuming. Thus, in our future research, automatic lesion segmentation will be performed using the U-shaped deep convolutional neural network.

CONCLUSION
In conclusion, our results indicated that texture features derived from kinetic parameter maps, calculated based on breast DCE-MRI, have the potential to be used as imaging biomarkers for distinguishing HER2-positive and HER2-negative breast cancer.
Further studies with larger sample sizes are necessary to verify the results of this study.

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 author.

ETHICS STATEMENT
Written informed consent was not obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.