DCE-MRI Texture Features for Early Prediction of Breast Cancer Therapy Response.

This study investigates the effectiveness of hundreds of texture features extracted from voxel-based dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) parametric maps for early prediction of breast cancer response to neoadjuvant chemotherapy (NAC). In total, 38 patients with breast cancer underwent DCE-MRI before (baseline) and after the first of the 6-8 NAC cycles. Quantitative pharmacokinetic (PK) parameters and semiquantitative metrics were estimated from DCE-MRI time-course data. The residual cancer burden (RCB) index value was computed based on pathological analysis of surgical specimens after NAC completion. In total, 1043 texture features were extracted from each of the 13 parametric maps of quantitative PK or semiquantitative metric, and their capabilities for early prediction of RCB were examined by correlating feature changes between the 2 MRI studies with RCB. There were 1069 pairs of feature-map combinations that showed effectiveness for response prediction with 4 correlation coefficients >0.7. The 3-dimensional gray-level cooccurrence matrix was the most effective feature extraction method for therapy response prediction, and, in general, the statistical features describing texture heterogeneity were the most effective features. Quantitative PK parameters, particularly those estimated with the shutter-speed model, were more likely to generate effective features for prediction response compared with the semiquantitative metrics. The best feature-map pair could predict pathologic complete response with 100% sensitivity and 100% specificity using our cohort. In conclusion, breast tumor heterogeneity in microvasculature as measured by texture features of voxel-based DCE-MRI parametric maps could be a useful biomarker for early prediction of NAC response.

Introduction

Neoadjuvant (preoperative) chemotherapy (NAC) was introduced in the 1970s, and over the past 2 decades, it has been established as a standard of care for patients with locally advanced breast cancer (LABC) for both initially operable and inoperable tumors (1–3). Compared with adjuvant (postoperative) chemotherapy, NAC has been shown to increase the breast-conserving surgery rate. Furthermore, the pathologic complete response (pCR) to NAC or minimal post-NAC residual disease has been found to be clearly correlated with disease-free and overall survival rates (1, 4–9). However, patients undergoing NAC do not always achieve pCR, and the pCR rate is reported to vary in the range of 6%–45% depending on breast cancer subtypes and treatment regimens (10–13). In the emerging era of precision medicine, early prediction of NAC response may allow rapid, personalized treatment regimen alterations for nonresponding patients with breast cancer, and spare them from potential short- and long-term toxicities associated with ineffective therapies. In addition, accurate evaluation of residual disease after NAC is vital for surgical decision-making and could result in surgical treatment plans that are more tailored to individual patients.

As a noninvasive imaging method for in vivo measurement of tissue microvascular perfusion and permeability, dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is increasingly used in clinical trial and research settings to assess breast cancer response to NAC (14). Changes in tumor size observed on the basis of DCE-MRI images are routinely used in clinical trials to assess tumor response to treatment. However, size change in response to therapy is often found to manifest later compared with changes in underlying tumor functions (15–18), such as vascularization and vascular permeability, cellularity, and metabolism. There is substantial literature showing that semiquantitative analysis (19–24) or quantitative pharmacokinetic (PK) analysis (25–35) of DCE-MRI time-course data can provide better early prediction of breast cancer pathologic response to NAC than image tumor size changes after 1–2 NAC cycles. We have recently shown (25, 36) that changes in tumor mean PK parameters after only 1 NAC cycle are good early predictors of pCR versus non-pCR.

However, the approaches to compute tumor mean DCE-MRI metrics, whether semiquantitative or quantitative, cannot capture the spatial heterogeneity of the tumor functions as measured by DCE-MRI, and therefore, valuable information could be missed during therapeutic response evaluation.

Texture analysis of semiquantitative or quantitative PK metrics from breast DCE-MRI has been shown to be effective in applications such as automatic lesion segmentation (37–39) and cancer diagnosis (40–46). Recently, Teruel et al. (47) presented a detailed analysis of 16 textural statistical features from T1-weighted DCE-MRI images that are capable of predicting early tumor response to NAC. However, this study used a 2-dimensional (2D) statistical texture description without taking advantage of the 3-dimensional (3D) information provided by the T1-weighted DCE-MRI. Moreover, these features were computed in a specific gray-level condition that reduced the number of potentially useful features. Golden et al. (48) used similar DCE-MRI texture features to predict therapeutic responses in terms of pCR, residual lymph node metastases, and residual tumor with lymph node metastases in patients with triple-negative breast cancer.

Here, we conducted a thorough analysis of hundreds of 3D statistical features extracted from parametric maps of semiquantitative and quantitative PK metrics that were obtained from a DCE-MRI study (25) of breast cancer response to NAC. We report our preliminary findings on the effectiveness of these features for early prediction of NAC response.

Materials and Methods

Patient Cohort and Study Schema

In this institutional review board-approved and HIPAA (Health Insurance Portability and Accountability Act of 1996)-compliant study, 38 women (age range, 27–79 years) diagnosed with grades 2–3 invasive breast tumors and scheduled to undergo NAC consented to participate in a longitudinal research MRI study that includes DCE-MRI. In total, 31 of the 38 women were treated with standard-of-care therapy regimens that included 4 cycles of doxorubicin–cyclophosphamide administration every 2 weeks followed by 4 cycles of a taxane every 2 weeks, or 6 cycles of the combination of all 3 drugs every 3 weeks. The other 7 patients were enrolled in the NAC ISPY-2 trial,[1] where patients were randomized to receive standard-of-care regimen or standard-of-care regimen plus experimental drugs. The ISPY-2 standard-of-care regimen started with a taxane administration followed by doxorubicin–cyclophosphamide administration. If used, the experimental drug was usually added to the taxane. In total, 4 of the 7 patients were placed in the treatment arm with experimental drugs—3 received neratinib, a tyrosine kinase inhibitor, and 1 received ganitumab, a human monoclonal antibody against type 1 insulin-like growth factor receptor. The targeted agent, trastuzumab, was added to the NAC regimen for tumors with positive HER2 (human epidermal growth factor receptor 2) receptor status (n = 23).

MRI examinations for this research study were performed before NAC (visit 1, V1), after 1 cycle of NAC (V2), at midpoint of NAC (V3, usually after 3 or 4 NAC cycles, or before change of NAC agents), and after NAC completion but before surgery (V4). For the studies, the MRI scan was undertaken at least 7 days after the administration of the previous NAC cycle to allow time for drug effects. This paper reports results for early prediction of NAC response, and thus, only the data from the V1 and V2 studies were used for texture analysis.

DCE-MRI Data Acquisition and Analysis

All breast MRI studies were performed using a 3T Siemens Tim Trio system with a body coil and a 4-channel bilateral phased-array breast coil as the radiofrequency transmitter and receiver, respectively. During each MRI session, following pilot scans and precontrast agent (CA) axial T2-weighted MRI with fat saturation and axial T1-weighted MRI without fat saturation, axial bilateral DCE-MRI images with fat saturation and full breast coverage were acquired with a 3D gradient echo-based TWIST (Time-resolved angiography WIth Stochastic Trajectories) sequence, which uses the strategy of k-space undersampling during acquisition and data sharing during reconstruction (49). DCE-MRI spatial resolution = 1.0 × 1.0 × 1.4 mm3 and temporal resolution = 14–20 seconds. Details of the acquisition parameters are described in Tudorica et al.'s study (25).

Breast tumor regions of interest (ROIs) were drawn by 2 experienced breast radiologists on postcontrast (∼90–120 seconds after the gadolinium CA injection) multisection DCE images covering the entire contrast-enhanced tumor. For voxels within the tumor ROI, the DCE-MRI time-course data were separately fitted with the following 2 PK models: 1 with a 1-compartment-2-parameter standard Tofts model (SM) (50, 51) and another with a 2-compartment-3-parameter shutter-speed model (SSM) (52). Both model fittings returned parameters K (rate constant for CA plasma-to-interstitium transfer) and v (volume fraction of extravascular and extracellular space), whereas the SSM also returned an additional parameter τ (mean intracellular water molecule lifetime), which is used to account for the effect of cross cell membrane water exchange kinetics (in the extravascular space) in the SSM (52). Details of PK data analysis and mathematical formulations for the SM and SSM are described in Tudorica et al.'s study (25). CA intravasation rate constant, k, was calculated as . The dK parameter, defined as , provided a measure of the exchange effect on K estimation (53) and was also calculated. Voxel-based multisection parametric maps of the PK parameters (within the tumor ROI) were generated following the SM and SSM fittings.

In addition, the following 5 semiquantitative metrics were quantified from the DCE-MRI signal time-course data, and voxel-based parametric maps were generated (41, 54, 55): where t0 and t are the first and last time points, respectively.

Initial enhancement (WashIn) describes the initial signal increase from the precontrast value to the postcontrast value, as defined by the following equation: where s0 is the average of the precontrast baseline signal intensity and s is the postcontrast signal intensity ∼90 seconds after completion of CA injection.

Signal enhancement ratio (SER) characterizes the postcontrast signal curve shape, as defined by the following equation: where s is the last signal value.

After initial enhancement (WashOut), the signal curve behavior is described after the initial phase of contrast enhancement, as follows:

Wash-in slope (SlopeIn) is a measure of the CA uptake rate, and is calculated using the following equation: where t is the time when s was measured.

Area under the curve (AUC), the signal integral value, is defined using the following equation:

Texture Feature Extraction

This section presents the mathematical definition of the various texture features under study. We extracted several sets of features. The first 3 sets are based on the gray-level cooccurrence matrix (GLCM) (54, 56–60), the gray-level run length matrix (RLM) (59–63), and the size zone matrix (SZM) (60, 64, 65), and all belong to the family of statistical matrices. Once these matrices are constructed, texture features [such as Haralick features (56) and moments] can be derived. The fourth set of features is based on a local binary pattern (LBP) representation (66, 67), and the fifth set is based on the morphological operation called pattern spectrum (PS) (68, 69). In addition to these advanced features, we also include classic texture features that are based on moments of the image intensity pattern (70).

Statistical matrices have been extensively used in texture characterization, the best known of which is the GLCM, which leads to the definition of Haralick's features (56). As a second-order statistic, and for an image f, the represents the joint probability that a pixel with gray-level i occurs jointly with another pixel having a gray-level j, for a given spatial offset Δ between the pair of pixels. For an offset , the GLCM is defined as follows:

By design, the GLCM is dependent on the offset and is therefore not rotation-invariant (Figure 1A). When using 8-connexity, this is addressed by computing the GLCM in 4 directions, corresponding to offsets , and , and the average matrix over all offsets can be used (57–59). In this study, we used the GLCM 3D formulation with 26-connexity, computed with 26 different offsets, and then we averaged the resulting matrices into a single matrix. On the basis of this matrix, we derived the following 16 second-order statistical features [known as the Haralick features (56)] for each image: average, contrast, correlation, energy (or second angular momentum), entropy, homogeneity, dissimilarity, inertia, variance, inverse difference momentum, sums average, sums variance, sources entropy, differences variances, differences entropy, and maximum probability. The detailed mathematical definition of these Haralick features is listed in the Supplemental Appendix.

Figure 1.

Example of statistical matrices construction. A 4-Gy levels image texture of size 4 × 4 (A). The gray-level cooccurrence matrix (GLCM) with Δ = (1, 0) (B). The run length matrix (RLM) with θ = 0° (C). The size zone matrix (SZM) using 8-connexity (D).

In addition to GLCM, another classical technique is a statistical matrix called the gray-level RLM (61), which has been extensively used for texture classification (60, 62, 63). Although the GLCM represents second-order statistical features, the RLM extracts higher-order statistical features. The matrix element counts the total number of runs in the image f with the gray-level g and run length l (ie, collinear pixels with the same intensity in the direction θ), as shown in Figure 1C. This method is particularly effective for periodic textures, and it complements the information provided by the GLCM. From the RLM, we can extract features as the moments of order varying from −2 to 2. Similar to the GLCM, this matrix requires computations using various directions to be rotation-invariant. In our study, we used the 3D RLM formulation computed in 26 directions.

Recently, Thibault et al. (60, 64, 65) introduced the gray-level SZM as an alternative to the joint RLM formulation. In this statistical matrix-based method, the value of the matrix element is equal to the number of zones in the image f with size s and gray level g (Figure 1D). The resulting matrix has a fixed number of rows equal to t (the total number of gray levels), and a dynamic number of columns, determined by the largest zone size and size quantization. The structure of SZM reflects the image texture—the more homogeneous the texture (large, flat zones with similar gray levels), the wider and flatter the matrix. From this statistical matrix representation, we can calculate all the second-order moments as compact texture features (62), plus 2 more features, which are specific weighted variances (64).

Unlike GLCM and RLM, which depend on the offset Δ and the orientation θ, respectively, the SZM matrix is invariant with respect to rotation and translation. In general, the RLM and GLCM are appropriate for periodic textures, whereas the SZM is more suitable for heterogeneous nonperiodic textures.

There are several SZM variants (65, 71). For example, the multiple gray-level SZMs incorporate the gray-level quantization. For an 8-bit image, it is computed from 8 SZMs for 8 different gray-levels (). The resulting matrices are combined by a weighted average using the following equation: where, is SZM calculated from quantized in N gray levels.

Remark.

By design, all these statistical matrices are sensitive to image acquisition noise. To improve their robustness, the number of gray levels can be reduced before a matrix is constructed. Different methods exist to reduce the number of gray levels to N possible values, for example, using a monotonically decreasing function or a cumulated histogram.

In addition to these statistical matrices-based features, we also extract a set of features known as the LBP (66). LBP is a simple yet very efficient texture operator, which labels the pixels of an image by thresholding the neighborhood of each pixel and produces a binary number. The LBP is widely used in pattern recognition because of its simplicity and computational speed and its efficacy in describing the local spatial structure of an image. The LBP provides a unique score for each pixel, and the scores' distribution represents the texture features. Practically, for an image f and a given pixel p with 8 ordered neighbors p,, with v = 1 if p ≥ p, or 0 otherwise.

One of the most important properties of the LBP operator is its robustness to monotonic gray-scale changes caused, for example, by illumination variations. In Ojala et al.'s study (67), a rotation-invariant formulation is presented, making this feature even more versatile. However, no proper definition of LBP in 3D exists. The most effective solution remains to compute the 2D LBP score on the 3 plans (XY, XZ, and YZ) for each pixel, and then to use all these values to fill the same histogram.

The next texture-characterization technique we included in our study is the PS (68, 69). The PS features describe the shape and size of structures in an n-dimensional signal. Measurement of the PS is based on morphological operations, which use a variety of structuring elements to filter a signal at multiple spatial scales. The PS is the combination of a granulometry, and its dual operation is the antigranulometry. A granulometry is the distribution study of all the object sizes present in an image. Formally, a granulometry is a family of morphological openings that depend on a positive parameter n, which expresses a size factor for a fixed structuring element. The granulometric analysis of an image f with respect to Γ involves evaluating each opening of size n with a measurement . The PS curve PS of f with respect to Γ and Φ [the antigranulometry, , a family of closings] is defined by the following normalized mapping:

The PS value for each size n is a probability density function (ie, a histogram), and it corresponds to a structure measurement: a peak in PS at a given scale n indicates the presence of many image structures of this scale or size. PS size distribution is a powerful gray-level and rotational-invariant texture descriptor (72).

In addition to these advanced texture features, we also added a set of classical moment-based texture features that include volume, average, and standard deviation. The statistical matrices were computed using different gray-level reduction values and algorithms. Counting all the texture features (GLCM, RLM, SZM, LBP, PS, and moments features), the total number of extracted features for each parametric map was 1043, distributed as shown in Table 1.

Table 1.

Texture Features' Distribution Among Different Feature-Extraction Techniques

Technique	Number of Features	%
Moments	6	0.58
GLCM	187	17.93
RLM	301	28.86
SZM	352	33.75
LBP	184	17.64
PS	13	1.24
Total	1043	100

Abbreviations: GLCM, gray-level cooccurrence matrix; RLM, run length matrix; SZM size zone matrix; LBP, local binary pattern; PS, pattern spectrum.

All the features presented in this paper are deterministic mathematical functions. They do not use any random parameter and, consequently, are reproducible.

The PS and LBP are powerful texture-characterization techniques. They describe a texture by providing distributions of patterns or sizes; therefore, their effectiveness comes from the combination of all the provided feature values. However, our small cohort size prevents simultaneous use of multiple features. In our study, we analyzed each feature individually (see section on Features' Evaluation below), by using one feature at a time for the subsequent RCB prediction.

Pathological Analysis

The status of pathologic response (to NAC) for each breast tumor was determined by pathological analysis of the post-NAC resection specimen. The following pathology parameters were measured from the resection specimen under light microscopy: cross-sectional tumor size in 2D, estimated invasive tumor cellular density, number of involved lymph nodes, and the greatest tumor dimension in the largest involved node. These measures were used to compute the RCB index value using an equation published by Symmans et al. (6). pCR is defined as the absence of residual invasive tumor with RCB = 0. A pathologic nonresponse (pNR) is defined as tumor cell density in a resection specimen that is either equal to or greater than the tumor cell density in a core biopsy specimen. Pathologic partial response is defined as findings intermediate between pCR and pNR. Non-pCR includes both pathologic partial response and pNR with RCB > 0.

Features' Evaluation

This paper aims to determine the capability of each map–feature pair to early predict the RCB index value for patients with LABC undergoing NAC. To do so, we have extracted a total of 1043 texture features in 3D from each of the 13 parametric maps, as presented in the previous section. They include Haralick features (from GLCM), RLM, SZM, LBP, PS, and basic moments, computed with various parameters and with gray-level reduction values and techniques. The capability of each feature is then evaluated individually using the following pipeline (see Algorithm 1):

For a given feature f and a given parametric map p, we extracted the feature for each patient at visit V1 and V2.

We computed the feature gradient, defined as the difference of feature values between V1 and V2.

Gradient outliers were then removed to improve the algorithm robustness.

The remaining gradients[2] were then centered (with an average of 0), whitened (with a normalized standard deviation of 1), and used as inputs in a ridge regression model (73) coupled with a leave-one-out cross validation to predict the corresponding normalized (to increase the computational precision) pathologically measured RCB index values.

The ridge regression was preferred over the classical linear regression because it adds a penalty term to the coefficients as a way of regularization, leading to reduced risk of overfitting and improved generalization capability. These penalty advantages are particularly beneficial because of the small number of patients available in our cohort. Moreover, by coupling the ridge regression with the leave-one-out cross validation protocol, the risk of overfitting is further reduced.

The predicted RCB index values using the feature were finally compared with the pathologically measured RCB values, and 4 correlations (using different paradigms) were computed for feature evaluation, namely, the Pearson product moment (linear), the Spearman rank-order (rank), the Kendall tau (rank), and the Goodman–Kruskal gamma (rank).

Results

According to the pathology analysis of the resection specimens, 9 patients were pCRs with RCB = 0, whereas the other 29 patients were non-pCRs, with RCB index values ranging from 0.43 to 4.1. These RCB index values were the inputs to the pipeline described by Algorithm 1 to evaluate the predictive ability for each of the 1043 3D texture features extracted from a total of 13 DCE-MRI parametric maps (both quantitative PK parameters and semiquantitative metrics). Consequently, a total of 13 559 map–feature pairs were investigated. To be considered as an effective predictor, a map–feature pair must have all 4 correlation coefficients >0.7. This “effectiveness” threshold value was empirically determined by the visual analysis of “true RCB versus predicted RCB curves”.

Among the 13 559 map–feature pairs investigated, 1069 (7.9%) were found to meet the effectiveness conditions.

Texture Features for Prediction of Therapy Response

The feature distribution among the 1069 feature–map pairs with all 4 correlation coefficients >0.7 is presented in Table 2. It is worth noting that the texture-extraction techniques (GLCM, RLM, SZM, etc.) provide an imbalanced number of features (Table 1). To fairly compare the effectiveness of these techniques, the distribution of the 1069 feature–map pairs was weighted according to the original distribution shown in Table 1—the higher the number of features provided by a technique, the lower the weight. For example, in Table 2, among the 1069 feature–map pairs, only 6 use moments, or 0.57%. The moments were also the technique generating the lowest number of features (6 features, representing 0.58% of all the generated features in Table 1). After correction through weighting, the 6 feature–map pairs using moments now represent 7.69% of the best pairs (Table 2). Because of this correction, the GLCM is the most frequently used feature-extraction technique, whereas the other techniques have comparable percent distribution.

Table 2.

Distribution of the Best Feature–Map Pairs With all 4 Correlations >0.7

Technique	Number of Features	%	Weighted %
Moments	6	0.57	7.69
GLCM	304	28.44	12.5
RLM	287	26.84	7.33
SZM	306	28.62	6.69
LBP	156	14.59	6.52
PS	10	0.94	5.92
Total	1069	100	N/A

Note: Features distribution among the 1069 best feature/map pairs, the matching percentage, and the weighted percentage. The weight computation is detailed in the section on “Texture Features for Prediction of Therapy Response”.

Abbreviations: GLCM, gray-level cooccurrence matrix; RLM, run length matrix; SZM size zone matrix; LBP, local binary pattern; PS, pattern spectrum.

MRI Metrics for Prediction of Therapy Response

Figure 2 shows how often quantitative PK parameters and semiquantitative metrics provided good texture features for early prediction of therapy response under the effectiveness condition of all 4 correlation coefficients >0.7. It appears that the SSM parametric maps (K, dK, τ, k, and v) are more likely to provide a good predictive feature than the SM PK parameters or the semiquantitative metrics.

Figure 2.

Percentage distribution of quantitative pharmacokinetic (PK) parameters and semiquantitative metrics among the 1069 best map–feature pairs (with all 4 correlations >0.7) for early prediction of therapy response.

After increasing the effectiveness condition threshold from 0.7 to 0.875, we find 9 feature–map pairs with all 4 correlation coefficients >0.875, as shown in Table 3. The GLCM (coupled with Haralick features) and the SSM PK parameters are the most frequently used features that are effective for prediction of NAC response, accounting for 6 out of the 9 best pairs. This is consistent with the results presented in Table 2 and Figure 2. Figure 3 shows examples of the correlations between the normalized pathologically measured RCB index values and the predicted RCB index values using + RLM + Gray-level nonuniformity feature (Figure 3A) and + Haralick + Contrast feature (Figure 3B) as predictive features, respectively. These correlations are nonlinear, even though the predictive models were built using the linear model of ridge regression. It is not mathematically possible to model a nonlinear problem using a linear approach, which sets a mathematical limit to the performance prediction. It also explains why the predicted RCB index values are accurately ranked, but the range of predicted values is smaller than that of the pathologically measured RCB index values. Nonetheless, even with a ridge regression, we can observe in Figure 3 that the predictive models can separate the patients with pCR (green dots; n = 9) from those with non-pCR (red dots; n = 29) with high accuracy. At 100% sensitivity (ie, correctly classifying all pCRs), the specificities are 100% and 96.7% for + RLM + Gray-level nonuniformity feature and + Haralick + Contrast feature, respectively. Figure 4 shows the map at visits V1 and V2 for a tumor with pCR (Figure 4A) and a tumor with non-pCR (Figure 4B). The Haralick contrast feature value increases by ∼450% for the tumor with pCR and by ∼30% for the tumor with non-pCR.

Table 3.

The 9 Best Map–Feature Pairs With all 4 Correlations >0.875

Map	Technique	Gray Level	Feature
Ktrans(SM)	Haralick	64	Entropy differences
Ktrans(SM)	RLM	N/A	Gray-level nonuniformity
Ktrans(SM)	RLM	N/A	Long-run emphasis
dKtrans	Haralick	128	Contrast
dKtrans	Haralick	128	Variance differences
dKtrans	Haralick	128	Inertia
τi	Haralick	8	Mean
υe(SSM)	Haralick	256	Contrast
υe(SSM)	Haralick	256	Inertia

Figure 3.

Fitted curves between the normalized pathologically measured residual cancer burden (RCB) index value and the predicted RCB: feature (A) and feature (B). Patients with pathologic complete response (pCR) and non-pCR are represented with green and red dots, respectively.

Figure 4.

Example that shows the changes in maps change examples for the following 2 tumors: one tumor with pCR at V1 (A) and V2 (B), and the other tumor with non-pCR at V1 (C) and V2 (D). The Haralick contrast feature value (Figure 3B) increased by ∼450% for the tumor with pCR and increased by ∼30% for the tumor with non-pCR.

Discussion

This preliminary study with a 38-patient cohort shows that changes in breast tumor heterogeneity as measured by changes in texture features of DCE-MRI voxel-based parametric maps can be useful markers for early prediction of breast cancer response to NAC, discriminating pCR versus non-pCR, as well as predicting low versus high post-NAC RCB index value. Although it clearly needs to be validated with larger patient populations, this noninvasive 3D imaging feature-extraction approach has the potential to become an important clinical tool in the emerging era of precision medicine to identify, in the early stages of treatment, nonresponding patients for alternative personalized therapy regimens, and stratify patients for better surgical decision-making, and after surgery care planning based on accurate prediction of RCB. For example, using the RLM + Gray-level nonuniformity feature of the map (Figure 3A), the 29 non-pCRs can be discriminated from the 9 pCRs with 100% sensitivity and specificity, and consequently, an AUC of receiver operating characteristic equal to 1, after only 1 out of 6–8 NAC cycles. In comparison, the mean parametric value previously suggested as an imaging marker for therapy response prediction (25) never reaches the effectiveness condition (all 4 correlations >0.7). When this feature is used to classify pCR versus non-pCR, the best AUC of the receiver operating characteristic is 0.91, when the feature is extracted from . These results suggest that, if texture feature analysis of DCE-MRI PK parameters had been used in clinical care, these 29 patients could potentially have been spared from the morbidity caused by continuous treatment with ineffective and toxic chemotherapy agents and could have been treated with different personalized therapy regimens or stratified for novel therapy trials. Furthermore, accurate prediction of RCB index value at the early stage of NAC course using DCE-MRI textures could potentially allow drug dose escalation to prevent heavy burden of residual disease and possibly enable longer survival.

One important observation in this study is the complete absence of the tumor size among the best features for early response prediction. This confirms the results presented in Martincich et al.'s study (24), suggesting that the tumor does not really shrink after a single NAC cycle, even for pCRs. However, the gray-level variance feature is present among the best features, indicating that changes in DCE-MRI texture heterogeneity are potentially important markers for early RCB prediction. Our findings suggest that one of the initial effects of NAC (before shrinkage) is changing the spatial heterogeneity of the tumor microenvironment, which, in this study, includes perfusion and permeability as measured by DCE-MRI. The importance of assessing changes in tumor heterogeneity for evaluation of therapy response is confirmed by the prevalence of predictive features from SZM and RLM, which are powerful techniques for characterizing the texture homogeneity/heterogeneity. Moreover, this finding is further strengthened by the high frequency of Haralick features such as entropy, variance, homogeneity, and contrast, which use different paradigms to describe texture homogeneity, as illustrated in Figure 4.

This study shows that texture features of DCE-MRI quantitative PK parameters are likely to be more useful than those of the semiquantitative metrics for early prediction of breast cancer NAC response, showing the benefit of performing PK modeling of the DCE time-course data. In addition, we found that the good predictive texture features are more likely to come from SSM than from SM parameters, although we observed comparable predictive capabilities between percent changes of tumor mean SM and SSM parameters after 1 NAC cycle in a separate study (25). This is possibly because of the fact that compared with SSM, the SM has a tendency to underestimate PK parameters (K and v) in malignant breast tissue, but not in benign or normal breast tissue (53). The SSM PK parametric map of the same breast tumor is expected to provide a larger dynamic range of the voxel parameter values than its SM counterpart, and thus, provide a more robust characterization of heterogeneity and its changes induced by therapy.

The results presented in this study were obtained using a single 3D texture feature extracted from a single parametric map, when previous work used a set of 2D features (47, 48) for early prediction response. Features in 3D offer more robust and efficient characterizations of early breast tumor changes than classical 2D features. This study, conducted on 1043 features and 13 maps, provides better characterization and evaluation of feature and map capabilities compared with the study that used 2D features. If a larger cohort of patients becomes available in the future, such information will allow us to build more complex models combining multiple features and maps for a better early prediction of breast cancer response to NAC.

This study has several limitations. First, the sample size of the study is small, and thus, it is important to validate the initial findings from this study with a larger patient cohort in the future. Second, mean tumor DCE-MRI parameter values were used in this study to assess breast tumor response to preoperative therapy. The tumor heterogeneity that is reflected in the imaging metrics, for example, in the K maps, was not captured in the mean DCE-MRI parameter values. The potential integration of mean values and texture features of DCE-MRI metrics may further improve the effectiveness of quantitative DCE-MRI for assessment of therapy response. The third limitation is the use of the ridge regression model (a linear technique) to predict a nonlinear phenomenon. A nonlinear regression technique trained with >1 feature from different maps, and on a larger cohort of patients, should increase the early response prediction. The last limitation is that the statistical matrices and the PS are not robust to resolution variations, and consequently, their efficacy may be reduced if the resolution is degraded. However, resolution degradation does not automatically lead to performance degradation. Further study is needed to elucidate the relationship between the 2.

In conclusion, we have investigated the capabilities of thousands of 3D texture features derived from different quantitative DCE-MRI maps for early prediction of NAC, on a cohort of 38 patients with LABC. Tumor texture heterogeneity changes captured by 3D statistical features measured using quantitative DCE-MRI parameters such as are promising markers for early prediction of pathologic response outcome.

Algorithm 1.

Pipeline of Feature Evaluation for RCB Index Value Prediction