A Machine-Learning Model for the Prediction of Triple-Negative Breast Cancer Based on Multiparameter MRI

Introduction

Breast cancer is the most common cancer among women¹ and is a group of heterogeneous diseases with various clinical behaviours and responses to treatment. Triple-negative breast cancers are characterized by a lack of expression of the estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2).² BRCA1 and BRCA2 (BRCA1/2) are the most important TNBC susceptibility genes, but only BRCA1 has been identified as a triple-negative susceptibility gene.³ Furthermore, compared with other breast cancer subtypes, pathogenic variations of PALB2, BRIP1, BARD1, and RAD51C are considered to be enriched in TNBC.³ Epidemiological data show that TNBC mainly occurs in young premenopausal women under 40 years of age, accounting for about 15–20% of all breast cancer patients.^4,5 Women who have TNBC cannot benefit from endocrine therapy or trastuzumab and have the worst prognosis.⁶ Thus, the accuracy of the identification of TNBC is essential.

Core needle biopsy (CNB) is the gold standard for identifying TNBC but has limitations. Several studies had indicated significant variability in the consistency between CNB samples and the assessment of tumor excision.^7–10 The coincidence rates of ER, PR and HER2 states ranged from 61% to 99%.^11,12 Intratumor heterogeneity could explain these results because of undersampling of the most informative areas by CNB.¹²

Dynamic contrast enhancement magnetic resonance imaging (DCE-MRI) has been extensively used in the diagnosis of breast cancer because its sensitivity is higher than other auxiliary tests. It can show malignant lesions that cannot be visualized by mammography or ultrasound.¹³ Previous studies mostly used functional and morphological parameters such as dynamic contrast enhancement patterns or apparent diffusion coefficient (ADC) to identify TNBC but rarely used semi-quantitative parameters to predict TNBC.^14–18 Angiogenesis plays a crucial role in the development and progression of breast cancer, and its pathological quantification often relies on indirect measures, such as alternative markers of microvascular density (MVD). Angiogenesis could be displayed by a breast MRI, enhancement features such as the enhancement ration phases or time to peak (TTP) of the breast MRI were associated with microvascular density of breast cancer. An adjacent vessel sign (AVS) was first defined by Carriero et al¹⁹ as a vessel (either artery or vein) in contact with a lesion or entering it. We used maximum intensity projection (MIP) series constructed from early postcontrast images in DCE-MRI.

In recent years, due to the development of statistical theory and the advancement of scientific and technological innovations, machine learning (ML) has attracted the attention and recognition of clinicians. ML methods have been applied to prediction of TNBC neoadjuvant chemotherapy response,²⁰ survival prediction,²¹ differentiation of TNBC and non-TNBC at the genetic level,²² and prediction of breast cancer outcomes in tumor tissue images.²³ However, few studies have analyzed the association between multiparametric MRI (semi-quantitative parameters and the number and diameter of the vessels) with TNBC. The purpose of this study was to develop a model based on multiparametric MRI for predicting TNBC.

Materials and Methods

Patients

Medical Ethics Committee of Tianjin Cancer Hospital approved this study (Approved No. of ethic committee: Ek2021235) and required waiving the requirement for informed patient consent. Because this was a retrospective study that analyzed pre-existing and unalterable historical data, without any intervention or experimental manipulation of the patients, and the use of these data would not directly affect the patients’ current health or living conditions. A total of 1353 consecutive breast cancer patients who underwent DCE-MRI without additional treatment before surgery from January 2019 to December 2020 were retrospectively analyzed. Surgeries were performed on all of the lesions following examination. The exclusion criteria were as follows: (1) patients with occult breast cancer and lymphoma; (2) patients treated with therapy before surgery; (3) poor image quality, incomplete image artifact or post-processing failure; (4) insufficient clinical-pathological information; (5) those who underwent vacuum-assisted or excisional biopsy for diagnosis (Figure 1).

Figure 1 Study flow chart. Abbreviation: DCE-MRI, dynamic contrast enhancement magnetic resonance imaging.

Histopathological Analysis

The histopathological and IHC results were retrieved from the Hospital Information System (HIS). Stained nuclei indicated negative ER or PR expression in <1% of cancer cells on ten high-power fields. In HER2 immunohistochemical staining, a score of 0 or 1+ was negative, 3+ was positive. Tumors scored at 2+ were equivocal and sent for fluorescence in-situ hybridisation to test. If the ratio of the HER2 gene signal to the chromosome 17 signal was less than 2, it was classified as negative.¹⁴ TNBC was defined as ER negative, PR negative and HER2 negative, or ER negative, PR negative, HER2 positive and FISH negative. Breast cancer was classified into two types based on the expression of ER/PR/HER2 in immunohistochemistry.

Imaging Protocol

All MRI examinations were performed on a 3.0 T scanner using a dedicated 8-channel phased-array breast coil (Discovery MR750, GE Healthcare) and a 1.5 T scanner using a dedicated 4-channel phased-array breast coil (Signa Infinity Excite II, GE Healthcare). The routine scans including T1-weighted (T1W) sequence, a fat-saturated T2-weighted (T2W) sequence and a diffusion-weighted sequence were performed using axial planes. The same sequence parameters were: field of view, 30cm×30cm; slice thickness, 5mm; slice gap, 0.5mm. On T1W and T2W images: matrix size, 384×224; number of excitations (NEX), 2 and flip angle, 90°. Repetition time/echo time (TR/TE) of T1W, 700ms/10ms; TR/TE of T2W, 4500ms/85ms.

The different parameters were: diffusion-weighted imaging (DWI) with readout-segmented echo-planar technique: TR/TE, 6300/64ms; b-values, 0, 500 and 1000 mm²/s; matrix, 128×128. The volume imaging for breast assessment (VIBRANT) bilateral breast imaging technique: field of view, 26cm×26cm; slice thickness, 1.8mm; flip angle, 15°; matrix size, 256×128; NEX, 1 and TR/TE, 6.1ms /2.9ms.

Six dynamic sequence analyses were performed before and after contrast injection. After unenhanced acquisition, 0.2mL/kg body weight of the contrast agent (Gd-DTPA) was injected at the rate of 2.0mL/s and then followed by an equal volume of saline solution. The acquisition was repeated five times, and each phase took 60s–90s.

Image Analysis

According to American College of Radiology BI-RADS MR lexicon, magnetic resonance images were independently reviewed by two radiologists (YXC, YBL, with 3 and 8 years of experience in breast MRI, respectively); in particular, they proceeded with the method of visual inspection of MR dynamic contrast-enhanced MIP of the subtraction between unenhanced images in the dynamic sequence and the first series of contrast-enhanced images. MIP reconstruction was performed on subtracted images. Both radiologists, blinded to clinicopathological findings, compared their conclusions and resolved inconsistencies by consensus. They assessed MR imaging findings including margins, enhancement, necrosis, edema, and ADC values. For multifocal or multicentric disease, the largest lesion was measured. The distance from the lesion to the nipple was measured on the reconstructed multi-planar reformation (MPR) by turning the rotation line to align the lesion with the nipple in the sagittal position and measuring the distance from the frontal edge of the lesion to the nipple in the transverse axis position.

Coronal, sagittal and transverse MIPs were prepared from subtracted MR images and 3D rotation using Functool II software on an Advantage Workstation (AW 4.7, GE Healthcare) to assess the AVS of breast cancer. The presence of vessels either entering the enhancing lesion or in contact with the lesion edge on MIP images was considered as the presence of AVS.^24–27 Count the number of these vessels, then find the largest vessel and measure its diameter.

DCE-MRI used semi-quantitative enhancement kinetic features, which had been hypothesized as the reflection of tumor angiogenesis to describe vascular changes after intravenous injection of a contrast agent.^28,29 These quantitative enhancement kinetic features were based on a region of interest (ROI) plotted on the brightest enhancement region to avoid bleeding and necrosis.

ROIs plotted to be as large as possible were placed inside the tumor on the first post-contrast image of DCE-MRI scans and ADC maps, avoiding areas of hemorrhage and necrosis. When measuring the gland enhancement ratio, the ROIs were drawn in the mirror image situation of the contralateral breast. If the patients underwent the radical mastectomy and had only one side of the breast, the ROIs plotted in the situation of the most apparent glandular enhancement in the ipsilateral breast. To assess each phase signal intensity increase, we calculated the enhancement for every post-contrast image: percentage enhancement at each phase(%)=[(post-contrast signal intensity−pre-contrast signal intensity)/pre-contrast signal intensity]×100%.³⁰ TTP, considered as the time (s) reached to achieve the maximum enhancement signal, was also assessed. Tumor size (the longest diameter of the tumor), the distance from the lesion to the nipple and the gland enhancement ratio were measured in the first phase of DCE-MRI. The gland enhancement ratio was calculated using the above formula. The maximum, minimum, mean, standard deviation and relative deviation of ADC were recorded on ADC maps.

XGBoost Model and Boruta Algorithm

The XGBoost algorithm optimizes the loss function using Newton’s method, achieved through a second-order Taylor expansion of the loss function combined with regularization terms to enhance traditional gradient boosting frameworks. During training, the objective function comprises two components: the loss function and the regularization term, which govern predictive accuracy and model complexity, respectively. Furthermore, the algorithm incorporates a “feature subsampling” strategy, where a random subset of features is selected for training each decision tree (analogous to the approach employed in Random Forest). This technique enhances the model’s generalization capability by promoting diversity among base learners, thereby effectively mitigating overfitting risks. The operational workflow of the XGBoost algorithm can be formalized as follows: feature vector with the corresponding (output) category yi³¹:

Feature selection: Boruta is a feature selection algorithm based on Random Forest. The Random Forest classifier is relatively fast, typically operable without parameter tuning, and provides numerical estimates of feature importance. The importance metric for a feature quantifies the loss in classification accuracy caused by the random permutation of its values across instances. A Z-score is calculated by dividing the mean loss by its standard deviation, serving as the importance measure. A shadow feature was obtained from the initial features by value resetting in a dataset copy. Shadow features were used as reference standards for selecting initial features. This method compared the Z-score between shadow features and each feature and evaluated them using a Random Forest classifier. A feature was chosen if its Z-score exceeded the highest Z-score of the shadow features.

Statistical Analysis

Continuous variables were expressed as medians with interquartile ranges (IQRs) and categorical variables as frequencies and percentages. The Mann–Whitney U-test compared continuous variables, while the Fisher exact test or Pearson χ2 test compared categorical variables. Statistical analysis was performed using R software (version 4.2.3). Univariate logistic regression identified variables associated with TNBC, and those with a p-value of < 0.1 were further analyzed using the Boruta method. An XGBoost model was built based on features selected by Boruta. Model performance for predicting TNBC was evaluated using the AUC with 95% confidence intervals (CI), and AUCs were compared using the DeLong test.

Results

Study Population

A total of 1376 lesions were detected in 1353 patients, including 23 patients with double breast cancer. Postoperative pathology confirmed 272 TNBC lesions and 1104 non-TNBC lesions. The age of patients ranged from 23 to 81 years old, and the median age was 48. The lesions were randomly divided into a primary cohort and a validation cohort at a ratio of 7:3. The clinical and MRI features of the primary cohort (963 lesions; median [interquartile range {IQR}] age, 48 [42–56] years) and validation cohort (413 lesions; median [interquartile range {IQR}] age, 47 [41–56] years) are shown in Tables 1 and 2 respectively. No significant differences were observed between the primary cohort and validation cohort. There was no difference between ADC values. (Supplementary Table 1).

Table 1 The Clinicopathologic of the Study Population

Table 2 The MRI Features of the Study Population

Identification of Predictive Factors

Tables 3 and 4 summarize respectively the clinicopathological and MRI characteristics of TNBC versus non-TNBC in the primary cohort. Detailed information of ADC values is present in Supplementary Table 2. Univariate analysis revealed that TNBC had higher clinical N stage (p < 0.001) and histological grade (p < 0.001) compared to non-TNBC. MRI findings showed TNBC more frequently exhibited rim enhancement (p < 0.001), edema (p < 0.001), necrosis (p < 0.001), high T2 signal (p = 0.04), and unifocal tumors (p = 0.04). Lesions in TNBC had larger maximum diameters (p < 0.001) and were further from the nipple (p = 0.03). Semi-quantitative parameters indicated lower enhancement ratios (p < 0.05) and TTP (p < 0.001) for TNBC, while gland enhancement ratio (p = 0.001) and wash-in slope (p = 0.02) were higher than non-TNBC. Additionally, the number and diameter of the vessels on MIP maps (p < 0.001) were higher and larger in TNBC compared to non-TNBC (Table 5).

Table 3 Comparison of the Clinicopathological Between TNBC and Non-TNBC in the Primary Cohort

Table 4 Comparison of the MRI Characteristics Between TNBC and Non-TNBC in the Primary Cohort

Table 5 Univariate Logistic Regression Analysis to Determine the Factors Related to TNBC

Development and Validation of the Predictive Model

The Boruta algorithm was used for feature selection to identify significant parameters for distinguishing TNBC from non-TNBC. The selection criterion required that a feature’s Z-score surpass the peak Z-score of the shadow features. (Figure 2a) The Boruta R package (accessed February 18, 2022)³² was employed, running the algorithm 100 times with 250 iterations each. (Figure 2b) Parameters selected in 90% or more of the runs were deemed significant, resulting in 13 critical parameters.

Figure 2 Boruta result plot. Feature selection analyzed by Boruta algorithm. Green and red boxplots represent confirmed and rejected attributes respectively. Yellow boxplot represent tentative attribute. Blue boxplots correspond to maximum, average, and minimal Z score of a shadow attribute (a). Z score evolution during Boruta run. Green lines correspond to confirmed attributes, yellow to tentative attribute, red to rejected ones and blue to respectively minimal, average and maximal shadow attribute importance. The vertical axis represents the importance of features, and the horizontal axis represents the number of runs of the classifier (b). Abbreviation: TTP, time to peak.

The XGBoost model was developed using critical variables identified by Boruta in the primary cohort (Figure 3). SHAP, a unified framework for model interpretation, assigned significance values to features, indicating their impact on model decisions. A higher SHAP value for a given feature in the model means a more significant impact on the decision of the model. The prediction model calculated the probability of TNBC based on necrosis, edema, lesion diameter, enhancement ratios, TTP, gland enhancement ratio, wash-in slope, and the vessel characteristics. Typical MRI examples of patients with TNBC and non-TNBC and their models are shown in Figure 4. Model performance was assessed using the AUC of the ROC curve, as well as accuracy, sensitivity, and specificity. The AUCs were 0.795 (95% CI: 0.758–0.832) for the primary cohort and 0.705 (95% CI: 0.640–0.770) for the validation cohort, demonstrating good predictive ability for TNBC (Figure 5).

Figure 3 SHAP model construction. The SHAP importance score of each feature in the XGBoost classifier. The darker the color (black), the greater the contribution to the prediction result (a). The relationship between feature value and the predicted probability is visually reflected by the color of the scatter plot of the SHAP summary plot, including positive and negative predictive effects. The color from light (yellow) to dark (black) indicates that the feature value is from high to low. SHAP value greater than zero, the greater the possibility of TNBC (b). Abbreviation: TTP, time to peak.

Figure 4 Pretreatment MR Images of non-TNBC patients and TNBC patients. (a–f) Axial fat-saturated T2W images, first sagittal post-contrast MR images acquired approximately 60s-90s after the contrast medium injection, kinetics, axial volume imaging for breast assessment, MR dynamic contrast-enhanced maximum intensity projection of non-TNBC and measurement of vessels. (g–l) Axial fat-saturated T2W images, first sagittal post-contrast MR images acquired approximately 60s-90s after the contrast medium injection, kinetics, axial volume imaging for breast assessment, MR dynamic contrast-enhanced maximum intensity projection of TNBC and measurement of vessels. (m) The waterfall chart of the XGBoost model predicted the possibility of non-TNBC in this case. (n) The waterfall chart of the XGBoost model (showing log odds ratio) predicted the possibility of TNBC in this case. Negative values indicate the probability that a lesion is a TNBC is less than 0.5. The features that are the yellow bars contribute to a prediction that a lesion is TNBC, while the features that are the purple bars contribute to a prediction that a lesion is non-TNBC. Features are ranked from top to bottom in terms of importance. Abbreviation: TTP, time to peak.

Figure 5 The ROC curves of the primary cohort’s model and the validation cohort’s model. The AUC is 0.795 for the primary cohort (a) and 0.705 for the validation cohort (b). Abbreviation: AUC, area under curve.

Discussion

The early identification of breast cancer subtypes is crucial for personalized treatment. It is one of the leading causes of cancer-related deaths.³³ This study utilized non-invasive MRI to differentiate between TNBC and non-TNBC and developed a prediction model based on multi-parametric MRI features. The XGBoost model has the characteristics of being fast and accurate. Compared with traditional machine learning algorithms, XGBoost is more original and better.³⁴ In this study, necrosis, edema, maximum lesion diameter, enhancement ratio in each stage, peak time, gland enhancement ratio, washout slope, number and diameter of blood were independently associated with TNBC. The model combining these predictors demonstrated good predictive performance, with AUCs of 0.795 (95% CI: 0.758–0.832) and 0.705 (95% CI: 0.640–0.770) in the primary and validation cohorts. This model demonstrates significant clinical value by enabling early TNBC identification on preoperative MRI, thus guiding personalized therapeutic strategies.

The vessel formation is crucial for the growth of malignant tumors including breast cancer.³⁵ TNBC exhibits the highest MVD among the breast cancer subtypes, partly due to elevated vascular growth factor expression.^36–39 Tumors with high MVD may benefit from antiangiogenic therapy.³⁵ TNBC is often associated with central necrosis and fibrosis, with rim enhancement as a characteristic feature.^14,15,40 Our study found rim enhancement in 45.31% of TNBC cases, aligning with the 41% to 80% range reported in the literature.¹⁵ Rim enhancement, a marker of high angiogenesis and tumor aggression, is linked to the overexpression of vascular endothelial growth factor and intratumoral hypoxia.^16,41 Reduced vessel numbers in central necrotic or fibrotic areas also contribute to rim enhancement.^28,42–44 Therefore, we suggest that rim enhancement in TNBC is associated with vessels.

While some studies have compared vessel differences between benign and malignant lesions,^45,46 there is limited focus on vessel differences across breast cancer subtypes. Chen et al⁴⁷ examined the role of vessel presence in distinguishing TNBC from non-TNBC using MRI. Our study included the number and diameter of the vessels in MIP. We found that the number of the vessels in TNBC was higher and the diameter of the vessels was larger than non-TNBC.

In this study, kinetic semi-quantitative parameters were crucial for identifying TNBC. Our data indicated that enhancement ratios at each stage, TTP, and wash-in slope were independently correlated with TNBC. Specifically, the enhancement ratios for TNBC were approximately 10% lower than that of non-TNBC. A recent study¹⁷ found that the initial enhancement ratio, indicative of early contrast wash-in, was elevated in HER2+ and TNBC compared to ER+/PR+ breast cancers. Previous research focused on early enhancement ratio and peak enhancement ratio with inconsistent results.^14,18 And Kim et al⁴⁸ reported that TTP had no concern with molecular typing. Our findings may differ due to a larger sample size. The enhancement ratio may be linked to tumor angiogenesis, with higher peak enhancement ratio associated with increased blood flow.^49,50 We speculated that TNBC’s high peripheral and low central vascular density contributes to higher TTP and wash-in slope and lower enhancement ratios at each stage.

High signal intensity on T2-weighted images was associated with necrosis within the tumor.⁴⁰ Similar to previous reports,^13,15,40,51 TNBCs were most frequently observed with high intratumoral signal intensity on T2-weighted images in our study.

Background parenchymal enhancement (BPE) in the breast, seen on MRI post-contrast injection, was influenced by endogenous estrogens, fibroglandular tissue, and menstrual phase. BPE was higher in women with more glandular tissue, extensive blood vessels and hormone sensitivity.⁵² While some studies^53,54 found higher BPE in basal-like and TNBC subtypes, others⁵⁵ reported lower BPE in TNBC compared to luminal A and B subtypes. However, Ha et al⁵⁶ demonstrated that there is no statistical correlation between BPE and molecular subtypes. This study agreed with Dilorenzo et al⁵⁴ in finding higher BPE in TNBC compared to non-TNBC, suggesting that BPE could be a valuable predictive factor for TNBC diagnosis despite variability in measurement standards. In this article, BPE was represented by gland enhancement ratio.

TNBC exhibits distinct MRI characteristics compared to other subtypes, including rim enhancement, unifocal lesions, higher histological grade, and high T2 signal intensity.^14,15,57 In this article, we found that the histological grade of TNBC had higher than that of non-TNBC, but it was not included in the model because the histological grade of some lesions was unknown. Edema, a marker of aggressiveness, is more common in TNBC.^58,59 Our study found that TNBC had a higher probability of edema. Chen et al⁵⁸ reported non-luminal breast cancer, including TNBC, showed more focal edema and higher breast edema scores (BES) than luminal breast cancer, with both measures correlating with tumor size and axillary lymph node metastasis. However, the link between tumor subtype and lymph node status is inconsistent. Some studies reported no significant differences in HER-2, ER, and PR positivity between ALNM and NALNM groups,^60,61 while others noted higher lymph node metastasis rates in TNBC.^62,63 This study suggested TNBC’s higher clinical N stage might be due to specific gene expression changes associated with survival and spread.⁶⁴

We observed that the distance from the TNBC lesion to the nipple was farther away. Kim et al⁶⁵ found that TNBC was located significantly closer to the chest wall. The results may account for our observation.

DWI is an MRI technique based on Brownian water motion, with the ADC reflecting proton diffusion. A recent meta-analysis indicated that ADC cannot predict breast cancer subtypes⁶⁶ and that some angiogenesis markers are not correlated with ADC values.⁶⁷ This study found no statistically significant ADC differences between TNBC and non-TNBC, suggesting that histopathological differences among subtypes may be too subtle for DWI to predict reliably.⁶⁶

There are some limitations to our study. Firstly, it is a monocentric and retrospective analysis; therefore, results need to be further validated in the prospective studies. Secondly, branching vessels may cause counting problems. In this study, we treated a Y-shaped vessel as a dual vessel if it satisfied the counting criteria. Finally, we manually mapped the ROI on the largest slice of the tumor and the statistics for the selected area might differ from the measurements for the entire tumor.

In conclusion, the model combining MRI parameters and machine learning method could predict the probability of TNBC before surgery and help implement targeted treatment for patients as early as possible.